chengxu.txt

常数元BEM求解平行双线电容和特性阻抗的电磁计算

clear;

format long;
a=1;
d=100;
e0=8.854e-12;
u=4*pi*10^(-7);
k=1;
for N=1:200;
%N=20;
    A=[ones(N,1);-1*ones(N,1)];
    
    x=[-0.86113,-0.33998,0.33998,0.86113];
    alfa=[0.34785,0.65214,0.65214,0.34785];
    
        
for i=1:2*N
    for j=1:2*N
        if (i<=N&&j<=N)||(i>N&&i<=2*N&&j>N&&j<=2*N)          %同一个圆
            if i==j             %i==j
                   H(i,j)=1/2;
               for m=1:4              %计算Gii
                  g(m)= -1/4/pi*log(a*a*((1-cos(pi/N))^2 +(x(m)*sin(pi/N))^2));
               end
                   G(i,j)=alfa(1)*g(1)+alfa(2)*g(2)+alfa(3)*g(3)+alfa(4)*g(4); 

              else                  %i~=j           
                 for m=1:4              %计算Gij,Hij
                    g(m)= -1/2/N*log(2*a*abs(sin(x(m)*pi/2/N-pi/N*(i-j))))*sqrt(abs(a^2+d^2/4-a*d*cos(x(m)*pi/N+pi*(2*j-2)/N)));
                    h(m)=-1/(4*a*N)*sin(pi*(j-i)/N-pi*x(m)/2/N)/sin(pi*(j-i)/N+pi*x(m)/2/N)*sqrt(abs(a^2+d^2/4-a*d*cos(x(m)*pi/N+pi*(2*j-2)/N)));
                 end
                   G(i,j)=alfa(1)*g(1)+alfa(2)*g(2)+alfa(3)*g(3)+alfa(4)*g(4);
                   H(i,j)=alfa(1)*h(1)+alfa(2)*h(2)+alfa(3)*h(3)+alfa(4)*h(4);
            end
        end
    end
end


for i=1:2*N
    for j=1:2*N
        if(i<=N&&j>N&&j<=2*N)||(i>N&&i<=2*N&&j<=N)   %%%%%%%不在同一圆上
            
            for m=1:4              %计算Gij
            g(m)= -1/2/N*log(sqrt(4*a*a*cos(pi*x(m)/2/N+pi*(i+j-2)/N)+d*d-2*a*d*(cos(pi*x(m)/N+pi*(2*j-2)/N)+cos(pi*(2*i-2)/N))))*sqrt(a*a+d*d/4-a*d*cos(pi*x(m)/N+pi*(2*j-2)/N));
            s(m)=(a*cos(2*pi*(i-j)/N)+a*cos(pi*x(m)/N+4*pi*(j-1)/N-d*cos(2*pi*(j-1)/N)));
            z(m)=(2*a*a*(1-cos(pi*x(m)/N+2*pi*(j-i)/N))-2*a*d*(cos(pi*x(m)/N+2*pi*(j-1)/N)+cos(2*pi*(i-1)/N))+d*d);
            h(m)=1/(2*N)*s(m)/z(m)*sqrt(a*a+d*d/4-a*d*cos(pi*x(m)/N+pi*(2*j-2)/N));
            end
              G(i,j)=alfa(1)*g(1)+alfa(2)*g(2)+alfa(3)*g(3)+alfa(4)*g(4); 
              H(i,j)=alfa(1)*h(1)+alfa(2)*h(2)+alfa(3)*h(3)+alfa(4)*h(4);
        end
    end
 end
 
B=-2*pi*e0*a*inv(G)*H*A;
b=0;
for n=1:N
        b=b+B(n);
end
%C=b/N
%Z=sqrt(u*e0)/C
%c=pi*e0/log((d+sqrt(d*d-4*a*a))/2/a)   %电容解析解
%z=sqrt(u*e0)/c
 
C(N)=b/N;
Z(N)=sqrt(u*e0)/C(N);
c(k)=pi*e0/log((d+sqrt(d*d-4*a*a))/2/a);    %电容解析解
z(k)=sqrt(u*e0)/c(k);
 
k=k+1;
end
N=1:200;
figure(1)
plot(N,C,'b',N,c,'r','lineWidth',2)
title('单位长度电容值随分段数N变化的曲线图')
xlabel('分段数N');
ylabel('电容值,单位(F)')
grid on
legend('边界元法','解析解')


figure(2)
plot(N,Z,'b',N,z,'r','lineWidth',2)
title('单位长度阻抗值随分段数N变化的曲线图')
xlabel('分段数N');
ylabel('阻抗值,单位(Ω)')
grid on
legend('边界元法','解析解')