用于二维和三维声学设计灵敏度分析的奇异边界法（Matlab代码实现）

%=========================================================================
% Essay topic:
%              Singular boundary method for 2D and 3D acoustic design 
%                             sensitivity analysis
%
%                                                      Date:    2022/05/30     
%=========================================================================
% benchmark example:
%                    3D acoustic sensitivity analysis for 
%                          Oscillating sphere
%                    ---Neumann boundary condition
%                    ---Acoustic sensitivity with respect 
%                            to the wave number k
%  %% Attention:
%               The results of the paper are the results obtained by 
%                    running the code in the software Matlab2019a.
%               ---Computations are implemented on a laptop with 
%                  an AMD Ryzen 7 4800H CPU at 2.90 GHz 
%                  on a 64-bit Windows 10 with 16G RAM.
%
%=========================================================================
clear all;
clc;
format long
%% Parameter settings
t0 = clock;              % Time Start 
n = 100;                 % Total number of soource points
m = 100;                 % The number of equal parts for design variable
a = 0.1;                 % Radius of pulsating cylinder
SS = 4*pi*a^2;
Lj = SS/n;               % The influence ranges of source point.
rho = 1.2;               % the air density
c = 341;                 % The speed of sound in air
v0=1;                    % the radial velocity
k = linspace(0.01,6,m);  % The value range of the wave number
dk = (k(end)-k(1))/(m-1);
xi = 5 ; yi = 0; zi = 0; % test point
i = sqrt(-1); 
%% Analytical solution of the acoustic sensitivity with respect to the wave number 
pt_exact = @(ro,k,a) i*rho*c*v0*a^2*exp(i*k.*(ro-a))./(ro.*(1-i*k*a).^2).*((1-i*k*a).^2+i*k*ro.*(1-i*k*a)+i*k*a);
%% Arrangement of Boundary Points
[nodes] = Creat_Sources(n);
xp = a.*nodes(:, 1);
yp = a.*nodes(:, 2);
zp = a.*nodes(:, 3);
n_x = -xp/a; n_y = -yp/a; n_z = -zp/a;
r=sqrt((xp-xp').^2+(yp-yp').^2+(zp-zp').^2);
%% Unknown coefficient solving process 
alpha = zeros(n,m);
for j=1:m
    lamda=i/(k(j)+1);
    b = i*rho*k(j)*c*v0*ones(n,1) ; % boundary condition (Nuemann)
    [~,F,~,H] = GFEH(k(j),r,xp,yp,zp,xp',yp',zp',n_x,n_y,n_z,n_x',n_y',n_z');
    A=F+lamda*H;
    [~,H0] = FH0(r,xp,yp,zp,xp',yp',zp',n_x,n_y,n_z,n_x',n_y',n_z');
    [~,E0] = GE0(r,xp,yp,zp,xp',yp',zp',n_x',n_y',n_z');
    E0(logical(eye(n)))=0; 
    Asum1 = sum(E0,2);
    H0(logical(eye(n)))=0; 
    Asum2 = sum(H0,2);
    H(logical(eye(n)))=0; 
    Asum3 = sum(H,2);
    %Source intensity factors (SIFs) 
    uii=1/(4*pi)*(i*k(j)+pi^4./(25*sqrt(Lj))+(log(pi))^2/SS);
    qii=1/Lj-Asum1;
    qii_BM=-Asum2+k(j)^2/2*uii;
    A(logical(eye(n))) = qii+lamda*qii_BM;  % Source intensity factors (SIFs) 
    alpha(:,j)=A\b;  
end
%% Central difference method
alpha_k(:,1) = (alpha(:,2)-alpha(:,1))/dk; %(-alpha(:,3)+4*alpha(:,2)-3*alpha(:,1))/(2*dk);
alpha_k(:,2:m-1) = (alpha(:,3:m)-alpha(:,1:m-2))/(2*dk);
alpha_k(:,m) = (alpha(:,m)-alpha(:,m-1))/dk;  %(3*alpha(:,m)-4*alpha(:,m-1)+alpha(:,m-2))/(2*dk);
%% Numerical solution and exact solution
Pe=zeros(m,1);
for j = 1:m
    lamda=i/(k(j)+1);
    lamda_dk=-i/(k(j)+1).^2;
    r=sqrt((xi-xp').^2+(yi-yp').^2+(zi-zp').^2);
    [G,E] = GE(k(j),r,xi,yi,zi,xp',yp',zp',n_x',n_y',n_z');
    [G_dk,~,E_dk,~] = GFEHk(k(j),r,xi,yi,zi,xp',yp',zp',n_x,n_y,n_z,n_x',n_y',n_z');
    Pe(j)=(G_dk+lamda_dk*E+lamda*E_dk)*alpha(:,j)+(G+lamda*E)*alpha_k(:,j);  %Numerical solution
end
ri=sqrt(xi.^2+yi.^2+zi.^2);
Exact = pt_exact(ri,k',a);  % Exact solution
%% Global erroe
Ge1 = sqrt(sum( (real(Pe -Exact)).^2 ))/sqrt(sum((real(Exact)).^2));
Ge2 = sqrt(sum( (imag(Pe -Exact)).^2 ))/sqrt(sum((imag(Exact)).^2));
fprintf('Numerical solutions(real part): %12.8e\n',real(Pe)) 
fprintf('Exact solutions(real part): %12.8e\n',real(Exact)) 
fprintf('Numerical solutions(imaginary part): %12.8e\n',imag(Pe)) 
fprintf('Exact solutions(imaginary part): %12.8e\n',imag(Exact)) 
fprintf('Global erroe(real part): %10.4e\n',Ge1) 
fprintf('Global erroe(imaginary part): %10.4e\n',Ge2) 
%% Comparision of numerical and analytical solutions by figures
figure(1)
plot(k,real(Exact),'bo--',k,real(Pe),'r*');
legend('Exact','BM-SBM');
xlabel('k'); 
ylabel('\partial\it{p}/\partial\it{k} (real part)');
figure(2)
plot(k,imag(Exact),'bo--',k,imag(Pe),'r*');
legend('Exact','BM-SBM');
xlabel('k');
ylabel('\partial\it{p}/\partial\it{k} (imaginary part)');
%% 
function [F0,H0] = FH0(r,x1,x2,x3,s1,s2,s3,n_x1,n_x2,n_x3,n_s1,n_s2,n_s3)
% G0 =1/(4*pi*r);
G0_dx1 = -1/(4*pi)*(x1-s1)./r.^3;
G0_dx2 = -1/(4*pi)*(x2-s2)./r.^3;
G0_dx3 = -1/(4*pi)*(x3-s3)./r.^3;
G0_ds1_dx1 = (r.^2-3*(x1-s1).^2)./(4*pi*r.^5);
G0_ds2_dx2 = (r.^2-3*(x2-s2).^2)./(4*pi*r.^5);
G0_ds3_dx3 = (r.^2-3*(x3-s3).^2)./(4*pi*r.^5);
G0_ds2_dx1 = -3*(x1-s1).*(x2-s2)./(4*pi*r.^5);
G0_ds1_dx2 = G0_ds2_dx1;
G0_ds3_dx1 = -3*(x1-s1).*(x3-s3)./(4*pi*r.^5);
G0_ds1_dx3 = G0_ds3_dx1;
G0_ds2_dx3 = -3*(x2-s2).*(x3-s3)./(4*pi*r.^5);
G0_ds3_dx2 = G0_ds2_dx3;
G0_dns_dx1 = G0_ds1_dx1.*n_s1 + G0_ds2_dx1.*n_s2 + G0_ds3_dx1.*n_s3;
G0_dns_dx2 = G0_ds1_dx2.*n_s1 + G0_ds2_dx2.*n_s2 + G0_ds3_dx2.*n_s3;
G0_dns_dx3 = G0_ds1_dx3.*n_s1 + G0_ds2_dx3.*n_s2 + G0_ds3_dx3.*n_s3;
F0 = G0_dx1.*n_x1+G0_dx2.*n_x2+G0_dx3.*n_x3;
H0 = G0_dns_dx1.*n_x1+G0_dns_dx2.*n_x2+G0_dns_dx3.*n_x3;
end
%% 
function [G,E] = GE(k,r,x1,x2,x3,s1,s2,s3,n_s1,n_s2,n_s3)
i = sqrt(-1);
G = exp(i*k*r)./(4*pi.*r);
G_ds1 = -exp(i*k*r).*(x1-s1).*(i*k*r-1)./(4*pi.*r.^3);
G_ds2 = -exp(i*k*r).*(x2-s2).*(i*k*r-1)./(4*pi.*r.^3);
G_ds3 = -exp(i*k*r).*(x3-s3).*(i*k*r-1)./(4*pi.*r.^3);
E = G_ds1.*n_s1+G_ds2.*n_s2+G_ds3.*n_s3;   %G_ns
end
%%  
function [G0,E0] = GE0(r,x1,x2,x3,s1,s2,s3,n_s1,n_s2,n_s3)
G0 =1./(4*pi*r);
G0_ds1 = 1/(4*pi)*(x1-s1)./r.^3;
G0_ds2 = 1/(4*pi)*(x2-s2)./r.^3;
G0_ds3 = 1/(4*pi)*(x3-s3)./r.^3;
E0 = G0_ds1.*n_s1+G0_ds2.*n_s2+G0_ds3.*n_s3;  %G0_ns
end
%%
function [G,F,E,H] = GFEH(k,r,x1,x2,x3,s1,s2,s3,n_x1,n_x2,n_x3,n_s1,n_s2,n_s3)
i = sqrt(-1);
G_ds1 = -exp(i*k*r).*(x1-s1).*(i*k*r-1)./(4*pi.*r.^3);
G_ds2 = -exp(i*k*r).*(x2-s2).*(i*k*r-1)./(4*pi.*r.^3);
G_ds3 = -exp(i*k*r).*(x3-s3).*(i*k*r-1)./(4*pi.*r.^3);
G_dx1 = -G_ds1;  
G_dx2 = -G_ds2;   
G_dx3 = -G_ds3;   
G_ds1_dx1 = -exp(i*k*r)./(4*pi.*r.^5).*((i*k*r-1).*((i*k*r-3).*(x1-s1).^2+r.^2)+i*k*r.*(x1-s1).^2);
G_ds2_dx2 = -exp(i*k*r)./(4*pi.*r.^5).*((i*k*r-1).*((i*k*r-3).*(x2-s2).^2+r.^2)+i*k*r.*(x2-s2).^2);
G_ds3_dx3 = -exp(i*k*r)./(4*pi.*r.^5).*((i*k*r-1).*((i*k*r-3).*(x3-s3).^2+r.^2)+i*k*r.*(x3-s3).^2);
G_ds2_dx1 = -exp(i*k*r).*(x1-s1).*(x2-s2)./(4*pi.*r.^5).*((i*k*r-1).*(i*k*r-3)+i*k*r);
G_ds1_dx2 = G_ds2_dx1;
G_ds3_dx1 = -exp(i*k*r).*(x1-s1).*(x3-s3)./(4*pi.*r.^5).*((i*k*r-1).*(i*k*r-3)+i*k*r);
G_ds1_dx3 = G_ds3_dx1;
G_ds2_dx3 = -exp(i*k*r).*(x2-s2).*(x3-s3)./(4*pi.*r.^5).*((i*k*r-1).*(i*k*r-3)+i*k*r);
G_ds3_dx2 = G_ds2_dx3;
G_dns_dx1 = G_ds1_dx1.*n_s1 + G_ds2_dx1.*n_s2 + G_ds3_dx1.*n_s3;
G_dns_dx2 = G_ds1_dx2.*n_s1 + G_ds2_dx2.*n_s2 + G_ds3_dx2.*n_s3;
G_dns_dx3 = G_ds1_dx3.*n_s1 + G_ds2_dx3.*n_s2 + G_ds3_dx3.*n_s3;
G = exp(i*k*r)./(4*pi.*r);
F = G_dx1.*n_x1+G_dx2.*n_x2+G_dx3.*n_x3;
E = G_ds1.*n_s1+G_ds2.*n_s2+G_ds3.*n_s3;
H = G_dns_dx1.*n_x1+G_dns_dx2.*n_x2+G_dns_dx3.*n_x3;
end
%%  
function [G_dk,F_dk,E_dk,H_dk] = GFEHk(k,r,x1,x2,x3,s1,s2,s3,n_x1,n_x2,n_x3,n_s1,n_s2,n_s3)
i = sqrt(-1);
G_dx1_dk = -k*(x1-s1).*exp(i*k*r)./(4*pi*r);
G_dx2_dk = -k*(x2-s2).*exp(i*k*r)./(4*pi*r);
G_dx3_dk = -k*(x3-s3).*exp(i*k*r)./(4*pi*r);
G_ds1_dk = -G_dx1_dk;  
G_ds2_dk = -G_dx2_dk;  
G_ds3_dk = -G_dx3_dk;  
G_ds1_dx1_dk = k.*exp(i*k*r)./(4*pi*r.^3).*((i*k*r-1).*(x1-s1).^2+r.^2);
G_ds2_dx2_dk = k.*exp(i*k*r)./(4*pi*r.^3).*((i*k*r-1).*(x2-s2).^2+r.^2);
G_ds3_dx3_dk = k.*exp(i*k*r)./(4*pi*r.^3).*((i*k*r-1).*(x3-s3).^2+r.^2);
G_ds1_dx2_dk = k.*exp(i*k*r)./(4*pi*r.^3).*(i*k*r-1).*(x1-s1).*(x2-s2);
G_ds1_dx3_dk = k.*exp(i*k*r)./(4*pi*r.^3).*(i*k*r-1).*(x1-s1).*(x3-s3);
G_ds2_dx3_dk = k.*exp(i*k*r)./(4*pi*r.^3).*(i*k*r-1).*(x2-s2).*(x3-s3);
G_ds2_dx1_dk = G_ds1_dx2_dk;
G_ds3_dx1_dk = G_ds1_dx3_dk;
G_ds3_dx2_dk = G_ds2_dx3_dk;
G_dns_dx1_dk = G_ds1_dx1_dk.*n_s1+G_ds2_dx1_dk.*n_s2+G_ds3_dx1_dk.*n_s3;
G_dns_dx2_dk = G_ds1_dx2_dk.*n_s1+G_ds2_dx2_dk.*n_s2+G_ds3_dx2_dk.*n_s3;
G_dns_dx3_dk = G_ds1_dx3_dk.*n_s1+G_ds2_dx3_dk.*n_s2+G_ds3_dx3_dk.*n_s3;
G_dk = i*r.*exp(i*k*r)./(4*pi*r);
F_dk = G_dx1_dk.*n_x1+G_dx2_dk.*n_x2+G_dx3_dk.*n_x3;
E_dk = G_ds1_dk.*n_s1+G_ds2_dk.*n_s2+G_ds3_dk.*n_s3;
H_dk = G_dns_dx1_dk.*n_x1+G_dns_dx2_dk.*n_x2+G_dns_dx3_dk.*n_x3;
end
%%  
function [rn, vn]=countnext(r,v,G)     
num=size(r,1);
dd=zeros(3,num,num); 
dd(1, :, :) = r(:, 1) - r(:, 1)';
dd(2, :, :) = r(:, 2) - r(:, 2)';
dd(3, :, :) = r(:, 3) - r(:, 3)';
L=sqrt(sum(dd.^2,1));     
L(L < G) = G;              
F=sum(dd./repmat(L.^3,[3 1 1]),3)';    
Fr=r.*repmat(dot(F,r,2),[1 3]);              
Fv=F-Fr; 
rn=r+v;  
rn=rn./repmat(sqrt(sum(rn.^2,2)),[1 3]);
vn=v+G*Fv;   
end
%%  
function [rn] = Creat_Sources(N)
a = linspace(0, 2*pi, N+1); a = a(1:N)';
b = linspace(-1, 1, N+1); b = b(1:N)';
b=asin(b);
r0=[cos(a).*cos(b), sin(a).*cos(b), sin(b)];
v0=zeros(size(r0));
G=1e-2;      
for ii=1:600      
    [rn,vn]=countnext(r0,v0,G);
    r0=rn; v0=vn;
end
end