Hi:
The integral to calculate acoustic field given below is a 2d convolution,
though it has non-zeros only at transducer surface. So one should be
able to use 2d fft to fastly calculate this convolution integral. But I
tried
the following MatLab script to firstly use quadrature method to get that
convolution and then use 2d fft to calculate that convolution. In the last
step, I plot results and compare them. But results are totally off.
What's wrong?
The following is the test using MatLab.
CI_D = sqrt(-1.0E0);
F = 0.5e6; % wave frequency
C = 1.5e5; % wave speed
k0 = 2.0E0*pi*F/C; % constant wave number
U0 = 5; % constant vibration velocity at transducer surface
% The circular transducer locates at z = 0 and with radius of R = 1.
R = 1.0E0; % transducer is a disk with radius R in xy-plane
% Computation domain is x from -2 to 2, y from -2 to 2, and z from 1 to 2.
% set computational domain nodes in x-dir
NoX = 32; xa = -2; xb = 2; NoMidX = ( NoX+mod(NoX,2) )/2;
LenX = xb-xa; dx = LenX/NoX; dfX = NoX/LenX;
x = linspace(xa, xb-dx, NoX);
% set computational domain nodes in y-dir
NoY = 32; ya = -2; yb = 2; NoMidY = ( NoY+mod(NoY,2) )/2;
LenY = yb-ya; dy = LenY/NoY; dfY = NoY/LenY;
y = linspace(ya, yb-dy, NoY);
% set computational domain nodes in z-dir
NoZ = 11; za = 1; zb = 2;
LenZ = zb-za; dz = LenZ/NoZ;
z = linspace(za, zb, NoZ);
% Method-1
% solve convolution using quad rule
% set discrete nodes that represent circular transducer disk
ni = 32; nj = 32;
xt = zeros(ni,nj); yt = xt;
TWOPI_D = 2.0E0*pi;
rho = linspace(0, ni-1, ni);
rho = rho*R/( ni-1 );
drho = abs( rho(1)-rho(2) );
psi = linspace(0, nj-1, nj);
psi = psi*TWOPI_D/( nj-1 );
dpsi = abs( psi(1)-psi(2) );
tmp = ones(ni,1);
for j = 1:1:nj
xt(:,j) = rho.*cos(psi(j))*tmp;
end % j loop
for j = 1:1:nj
yt(:,j) = rho.*sin(psi(j))*tmp;
end % j loop
dfactor = drho*dpsi;
dA = dfactor*ones(ni,nj);
for i = 1:1:ni
dA(i,:) = rho(i)*dA(i,:);
end % i loop
% calculate convolution integral using quad rule
P = ( 1.0E0+CI_D )*zeros(ni,nj);
for k = 1:1:NoZ
vec_z = z(k)*ones(ni,nj);
for j = 1:1:NoY
vec_y = y(j)-yt;
for i = 1:1:NoX
vec_x = x(i)-xt;
r = sqrt( vec_x.*vec_x+vec_y.*vec_y+vec_z.*vec_z );
P(i,j,k) = sum( sum( U0.*dA.*exp( -CI_D*k0*r )./r ) );
end % i loop
end % j loop
end % k loop
% calculate power field
I = real( P.*conj(P) );
% Method-2
% solve convolution using fft
U4Fft = zeros(NoX,NoY); Ker4Fft = ( 1.0E0+CI_D )*U4Fft;
P4Fft = ( 1.0E0+CI_D )*U4Fft;
for k = 1:1:NoZ
for j = 1:1:NoY
for i = 1:1:NoX
tmp = x(i)*x(i)+y(j)*y(j);
rr = sqrt(tmp);
r = sqrt( tmp+z(k)*z(k) );
% assign vibration velocity and wave kernel
% ONLY at circular transducer disk
if abs(rr) < R | abs( rr-R ) <= eps
U4Fft(i,j) = U0;
Ker4Fft(i,j) = exp( -CI_D*k0*r )/r;
else
U4Fft(i,j) = 0.0;
Ker4Fft(i,j) = 0.0;
end % if abs(r) <= R
end % i loop
end % j loop
% 2D (xy) FFT for vibration velocity and for wave kernel
UFft = dx*dy*fft2(U4Fft);
KerFft = dx*dy*fft2(Ker4Fft);
% 2D inverse FFT for product of the above
P4Fft(:,:,k) = dfX*dfY*fftshift(real(ifft2( UFft.*KerFft )));
end % k loop
% calculate the power filed
I4Fft = real( P4Fft.*conj(P4Fft) );
% plot results
figure
subplot(211)
plot(z, squeeze(I(NoMidX,NoMidY,:)), z, squeeze(I(NoMidX,NoMidY,:)), '.')
title('Quad method')
subplot(212)
plot(z, squeeze(I4Fft(NoMidX,NoMidY,:)), z, squeeze(I4Fft(NoMidX,NoMidY,:)),
'.')
title('Fft method')
.
|
