一、简介

基于matlab二维边界单元法计算腐蚀电位

二、源代码

clear;

fid=fopen('input.dat','r');
indat=fscanf(fid,'%g%g%d%g',[4,inf]);
indat=indat';
xb=indat(:,1);
yb=indat(:,2);
bt=indat(:,3);
bv=indat(:,4);
n=length(xb)-1;

for i=1:n
    xm(i)=0.5d0*(xb(i)+xb(i+1));
    ym(i)=0.5d0*(yb(i)+yb(i+1));
    lm(i)=sqrt((xb(i+1)-xb(i))^2d0+(yb(i+1)-yb(i))^2d0);
    nx(i)=(yb(i+1)-yb(i))/lm(i);
    ny(i)=(xb(i)-xb(i+1))/lm(i);
end

for m=1:n
    b(m)=0d0;
    for k=1:n
        if(k==m)
            G=0.0;
            F=lm(k)/(2.0*pi)*(log(lm(k)/2.0)-1.0);
            del=1.0;
        else
            [F,G]=findfg(xm(m),ym(m),xb(k),yb(k),nx(k),ny(k),lm(k));
            del=0.0;
        end
        if(bt(k)==0)
            A(m,k)=-F;
            b(m)=b(m)+bv(k)*(-G+0.5d0*del);
        else
            A(m,k)=G-0.5d0*del;
            b(m)=b(m)+bv(k)*F;
        end
    end
end
z=A\b';

for m=1:n
   u(m)=(1-bt(m))*bv(m)+bt(m)*z(m);
   q(m)=(1-bt(m))*z(m)+bt(m)*bv(m);
end

for j=1:9
   y(j)=0.1*j;
   for i=1:9
      x(i)=0.1*i;
      s(j,i)=0d0;
      for k=1:n
         [F,G]=findfg(x(i),y(j),xb(k),yb(k),nx(k),ny(k),lm(k));
         s(j,i)=s(j,i)+u(k)*G-q(k)*F;
      end
   end
end

三、运行结果

四、备注

版本：2014a