%=========================================================================%
% HIE2DCYL.M %
% %
% This is the code for calculating solid angle C, surface pressure ps, %
% and field pressure pf coming out from a pulsating cylinder with radius %
% of r and normal velocity vn using the solution of Helmholtz Integral %
% Equation. For verifying the calculation result, it will also presents %
% the analytical solution. In case of non-uniqueness problem, r = 1 m %
% and f = 131 Hz, up to 10 CHIEF points can be implemented. %
% %
% The code is written in a very simple and traditional way so that any %
% beginner in MATLAB programming would be able to understand it easily. %
% %
% Developed by Agustinus Oey <oeyaugust@gmail.com> %
% As a short project of Acoustics Laboratory Summer Camp 2006, Week 5-6 %
% Center of Noise and Vibration Control (NoViC) %
% Department of Mechanical Engineering %
% Korea Advanced Institute of Science and Technology (KAIST) %
% Daejeon, South Korea %
% %
% Ref., T.W. Wu, "Boundary Element Acoustics: Fundamentals and Computer %
% Codes", WIT Press, Computational Mechanics, 2001" %
%=========================================================================%
% preparing the environment
close all;
clear;
clc;
% defining medium and wave parameter
rho0=1.21;
cair=343;
% asking problem parameter
r=0;
while r<0.5 | r>2,
r = input('Cylinder radius [0.5 ... 2 m]? ');
end
nnode=0;
while nnode<16 | nnode>128,
nnode = floor(input('Number of node [16 ... 128]? '));
end
freq=0;
while freq<20 | freq>500,
freq = input('Pulsating frequency [25 ... 250 Hz]? ');
end
vr=10;
while vr<-2 | vr>2,
vr = input('Surface velocity vr [-2 ... 2 m/s]? ');
end
disp('-');
vn=-ones(1,nnode)'*vr;
% calculating angular frequency and wave number
omega=2*pi*freq;
kwave=omega/cair;
% calculating angle spacing inside the cylinder
dteta=360/nnode;
teta=-360:dteta:0-dteta;
teta=-teta/90*1.5708;
% calculating coordinate of nodes
xnode=r*cos(teta);
ynode=r*sin(teta);
clear teta dteta
% defining basic linier shape function
tshp1=[0.5 -0.5];
tshp2=[0.5 0.5];
% calculating x and y as functions of zeta
xzeta=zeros(2, nnode);
yzeta=zeros(2, nnode);
for k=1:nnode-1,
xzeta(:,k)=(xnode(k)*tshp1+xnode(k+1)*tshp2)';
yzeta(:,k)=(ynode(k)*tshp1+ynode(k+1)*tshp2)';
end
k=nnode;
xzeta(:,k)=(xnode(k)*tshp1+xnode(1)*tshp2)';
yzeta(:,k)=(ynode(k)*tshp1+ynode(1)*tshp2)';
% calculating dN1/dzeta and dN2/dzeta
dtshp1=-0.5;
dtshp2= 0.5;
dxzta=zeros(1, nnode);
dyzta=zeros(1, nnode);
for k=1:nnode-1,
dxzta(k)=xnode(k)*dtshp1+xnode(k+1)*dtshp2;
dyzta(k)=ynode(k)*dtshp1+ynode(k+1)*dtshp2;
end
k=nnode;
dxzta(k)=xnode(k)*dtshp1+xnode(1)*dtshp2;
dyzta(k)=ynode(k)*dtshp1+ynode(1)*dtshp2;
% calculating Jacobian and normal vectors
Jacb=sqrt(dxzta.^2+dyzta.^2);
nvect=[dyzta./Jacb; -dxzta./Jacb];
clear dtshp1 dtshp2 dxzta dyzta
% defining spacing and weighting of 6 points Gaussian Quadrature Int.
ngaus=6;
wzeta=[-0.93246951 -0.66120939 -0.23861918 ...
0.23861918 0.66120939 0.93246951];
wgaus=[0.17132449 0.36076157 0.46791393 ...
0.46791393 0.36076157 0.17132449];
% calculating [C], [H*], and [G]
czero=eye(nnode);
hzero=zeros(nnode, nnode);
gzero=zeros(nnode, nnode);
for k=1:nnode,
xp=xnode(k);
yp=ynode(k);
% calculate Ck, hkj, and gkj
for j=1:nnode,
for m=1:ngaus,
xpoin=[1 wzeta(m)]*xzeta(:,j);
ypoin=[1 wzeta(m)]*yzeta(:,j);
rdist=sqrt((xpoin-xp)^2+(ypoin-yp)^2);
dpsdr=-1/(2*pi*rdist);
drdn=[xpoin-xp ypoin-yp]*nvect(:,j)/rdist;
% arranging Ck into [C]
czero(k,k)=czero(k,k)-dpsdr*drdn*Jacb(j)*wgaus(m);
psi=-i*besselh(0,2,kwave*rdist)/4;
dpsidn=i*kwave*besselh(1,2,kwave*rdist)*drdn/4;
n1=tshp1*[1; wzeta(m)];
n2=tshp2*[1; wzeta(m)];
% arranging hkj into [H*] and gkj into [G]
hzero(k,j)=hzero(k,j)+dpsidn*n1*Jacb(j)*wgaus(m);
gzero(k,j)=gzero(k,j)-i*rho0*omega*psi*n1*Jacb(j)*wgaus(m);
if j==nnode
hzero(k,1)=hzero(k,1)+ ...
dpsidn*n2*Jacb(1)*wgaus(m);
gzero(k,1)=gzero(k,1)- ...
i*rho0*omega*psi*n2*Jacb(1)*wgaus(m);
else
hzero(k,j+1)=hzero(k,j+1)+ ...
dpsidn*n2*Jacb(j+1)*wgaus(m);
gzero(k,j+1)=gzero(k,j+1)- ...
i*rho0*omega*psi*n2*Jacb(j+1)*wgaus(m);
end
end
end
end
% calculating condition number Cn
Cn=max(svd(czero+hzero))/min(svd(czero+hzero));
% offering CHIEF point if needed
nchif=0;
teks=['Condition number is ' num2str(Cn)];
disp(teks);
clear teks
k = floor(input('Add CHIEF points [0 or 1]? '));
disp('-');
if k==1,
% ploting solid angle
while nchif<1 | nchif>10,
nchif = floor(input('Number of CHIEF points [1 .. 10]? '));
end
disp('-');
xchif=0:r/nchif:r-r/nchif;
for k=1:nchif,
xchif(k)=xchif(k)*(-1)^(k-1);
end
ychif=zeros(1,nchif);
% calculating [Hc*] and [Gc] from CHIEF points
hchif=zeros(nchif, nnode);
gchif=zeros(nchif, nnode);
for k=1:nchif,
xp=xchif(k);
yp=ychif(k);
% calculating hc and gc from CHIEF points
for j=1:nnode,
for m=1:ngaus,
xpoin=[1 wzeta(m)]*xzeta(:,j);
ypoin=[1 wzeta(m)]*yzeta(:,j);
rdist=sqrt((xpoin-xp)^2+(ypoin-yp)^2);
drdn=[xpoin-xp ypoin-yp]*nvect(:,j)/rdist;
psi=i*besselh(0,2,kwave*rdist)/4;
dpsidn=i*kwave*besselh(1,2,kwave*rdist)*drdn/4;
n1=tshp1*[1; wzeta(m)];
n2=tshp2*[1; wzeta(m)];
% arranging hc into [Hc*] and gc into [Gc]
hchif(k,j)=hchif(k,j)+dpsidn*n1*Jacb(j)*wgaus(m);
gchif(k,j)=gchif(k,j)+i*rho0*omega*psi*n1*Jacb(j)*wgaus(m);
if j==nnode
hchif(k,1)=hchif(k,1) + ...
dpsidn*n2*Jacb(1)*wgaus(m);
gchif(k,1)=gchif(k,1) + ...
i*rho0*omega*psi*n2*Jacb(1)*wgaus(m);
else
hchif(k,j+1)=hchif(k,j+1) + ...
dpsidn*n2*Jacb(j+1)*wgaus(m);
gchif(k,j+1)=gchif(k,j+1) + ...
i*rho0*omega*psi*n2*Jacb(j+1)*wgaus(m);
end
end
end
end
Hchif = [hzero; hchif];
Gchif = [gzero; gchif];
Cchif = [czero; zeros(nchif, nnode)];
Cn=max(svd(Cchif+Hchif))/min(svd(Cchif+Hchif));
teks=['New condition number is ' num2str(Cn)];
disp(teks);
clear teks
% calculating surface pressure HIE solution ps with CHIEF points
ps=pinv(Cchif+Hchif)*(Gchif*vn);
else
% calculating surface pressure HIE solution ps w/o CHIEF points
ps=inv(czero+hzero)*(gzero*vn);
end
% calculating surface pressure analytical solution psa
psa=i*rho0*cair*vr*besselh(0,2,kwave*1)/besselh(1,2,kwave*1);
% ploting elements, nodes, normal vector, and CHIEF points
k=0;
k = floor(input('Ploting the geometry [0 or 1]? '));
disp('-');
if k==1,
%ploting the elements, nodes, and CHIEF points
if nchif>0 & nchif<11,
plot([xnode xnode(1)], [ynode ynode(1)], 'r', ...
xnode, ynode, '*', ...
xchif, ychif, 'x')
else
plot([xnode xnode(1)], [ynode ynode(1)], 'r', ...
xnode, ynode, '*')
end
%lining normal vectors
for k=1:nnode-1,
xelm(k)=0.5*(xnode(k)+xnode(k+1));
yelm(k)=0.5*(ynode(k)+ynode(k+1));
end
k=nnode;
xelm(k)=0.5*(xnode(k)+xnode(1));
yelm(k)=0.5*(ynode(k)+ynode(1));
for k=1:nnode,
line([xelm(k) xelm(k)+0.05*nvect(1,k)], ...
[yelm(k) yelm(k)+0.05*nvect(2,k)])
end
clear xelm yelm
title('Geometry of Cylinder')
xlabel('x axis (m)')
ylabel('y axis (m)')
axis square
disp('Legend: asterix - node, green x - CHIEF point');
disp(' red line - surface element, blue line - normal vector');
disp('-');
end
% ploting solid angles and surface pressure, if needed
k=0;
k = floor(input('Ploting Cp, ps, and pf [0 or 1]? '));
if k==1,
% ploting solid angle
figure
subplot(3,1,1)
k=1:nnode;
plot(k, diag(czero), k, 0.5*ones(1,nnode),'r')
axis([1 nnode 0.4950 1.01*max(diag(czero))])
legend('HIE', 'Analytic', 4)
title('Solid Angle C(p) and Surface Pressure ps')
ylabel('C(p)')
% ploting real surface pressure
subplot(3,1,2)
plot(k, real(ps), k, real(psa)*ones(1,nnode),'r')
legend('HIE', 'Analytic', 4)
ylabel('real ps (Pa)')
xaxis=[real(psa) min(real(ps)) max(real(ps))];
axis([1 nnode 0.999*min(xaxis) 1.001*max(xaxis)])
% ploting real surface pressure
subplot(3,1,3)
plot(k, imag(ps), k, imag(psa)*ones(1,nnode),'r')
legend('HIE', 'Analytic', 4)
xlabel('Number of Node (n)')
ylabel('imaginer ps (Pa)')
xaxis=[imag(psa) min(imag(ps)) max(imag(ps))];
axis([1 nnode 0.999*min(xaxis) 1.001*max(xaxis)])
clear xaxis
end
% noticing for next process
disp('-');
disp('Calculating field pressure pf.');
disp('Press any key to continue.');
pause
% defining field points
xfild=r:0.25:20;
yfild=0*xfild;
% calculating [H*] and [G] from field points
npoin=length(xfild);
hfild=zeros(npoin, nnode);
gfild=zeros(npoin, nnode);
for k=1:npoin,
xp=xfild(k);
yp=yfild(k);
% calculating hc and gc from field points
for j=1:nnode,
for m=1:ngaus,
xpoin=[1 wzeta(m)]*xzeta(:,j);
ypoin=[1 wzeta(m)]*yzeta(:,j);
rdist=sqrt((xpoin-xp)^2+(ypoin-yp)^2);
drdn=[xpoin-xp ypoin-yp]*nvect(:,j)/rdist;
psi=i*besselh(0,2,kwave*rdist)/4;
dpsidn=i*kwave*besselh(1,2,kwave*rdist)*drdn/4;
n1=tshp1*[1; wzeta(m)];
n2=tshp2*[1; wzeta(m)];
% arranging hc into [Hc*] and gc into [Gc]
hfild(k,j)=hfild(k,j)+dpsidn*n1*Jacb(j)*wgaus(m);
gfild(k,j)=gfild(k,j)+i*rho0*omega*psi*n1*Jacb(j)*wgaus(m);
if j==nnode
hfild(k,1)=hfild(k,1) + ...
dpsidn*n2*Jacb(1)*wgaus(m);
gfild(k,1)=gfild(k,1) + ...
i*rho0*omega*psi*n2*Jacb(1)*wgaus(m);
else
hfild(k,j+1)=hfild(k,j+1) + ...
dpsidn*n2*Jacb(j+1)*wgaus(m);
gfild(k,j+1)=gfild(k,j+1) + ...
i*rho0*omega*psi*n2*Jacb(j+1)*wgaus(m);
end
end
end
end
% calculating field pressure HIE solution pf
pf=gfild*vn-hfild*ps;
% calculating surface pressure analytical solution psa
pfa=i*rho0*cair*vr*besselh(0,2,kwave*xfild)/besselh(1,2,kwave*r);
% ploting pf and pfa
figure;
plot(xfild, real(pfa), 'b', xfild, real(pf), ...
'b*', xfild, imag(pfa), 'r', xfild, imag(pf), 'r*')
legend('re(pf), analytic', 're(pf), HIE', 'im(pf), analytic', 'im(pf), HIE')
title('Field Pressure (pf) in Acoustic Domain')
xlabel('distance from the surface (m)')
ylabel('real and imaginer pressure (Pa)')
axis tight
disp('-');
disp('End of calculation, thank you.');
