差分法による電位分布計算

denkai.m

N = 300; %r方向の格子数
M = 300; %z方向の格子数
RMAX = 30; %r座標の最大値
ZMAX = 30; %z座標の最大値
Conv = 1e-6; %収束判定誤差
dr = RMAX / N; %r方向の刻み幅
dz = ZMAX / M; %z方向の刻み幅
phi = zeros(N,M); %0埋め
flag = -10*ones(N,M); %非境界部分の-10埋め

for i = 1:N
    for j = 1:M
        %格子点(i,j)の座標(r,z)
        r = i*dr;
        z = j*dz;
        %ディリクレ条件(1kV): 高圧電極
        if ((r<1 && z > 16) || r^2+(z-14)^2<9)
            flag(i,j)=-3;
            phi(i,j)=1.0;
        %ディリクレ条件(0kV): 低圧電極、接地面、r遠方境界
        elseif ((r<1 && z < 5) || r^2+(z-7)^2<9 || j==1 || i==N)
            flag(i,j)=-3;
        %z方向ノイマン条件
        elseif (j == M)
            flag(i,j)=-1;
        %r方向ノイマン条件
        elseif (i == 1)
            flag(i,j)=-2;
        end
    end
end

MaxPhi = 1.0;
MaxErr = 1e20;

loop = 0;
while MaxErr > Conv
    Prev_phi=phi;
    for i = 1:N
        for j = 1:M 
            r = i*dr; %格子点(i,j)の座標(r,z)
            %非境界のphi
            if (flag(i,j)==-10)
                phi(i,j)=0.25*(Prev_phi(i+1,j)+Prev_phi(i-1,j)+Prev_phi(i,j+1)+Prev_phi(i,j-1)+0.5*dr*(Prev_phi(i+1,j)-Prev_phi(i-1,j))/r);
            %z方向ノイマン境界
            elseif (flag(i,j)==-1) 
                phi(i,j)=0.25*(Prev_phi(i+1,j)+Prev_phi(i-1,j)+2*Prev_phi(i,j-1)+0.5*dr*(Prev_phi(i+1,j)-Prev_phi(i-1,j))/r);
            %r方向ノイマン条件
            elseif (flag(i,j)==-2)
                phi(i,j)=0.25*(2*Prev_phi(i+1,j)+Prev_phi(i,j+1)+Prev_phi(i,j-1));
            end
        end
    end
    %前回の行列との誤差。MaxPhiで規格化
    Err = (abs(phi-Prev_phi))/MaxPhi;
    MaxErr=max(max(Err));%行列の最大の要素
    loop = loop + 1;
    %1000ループごとにMaxMaxErr表示
    if rem(loop,1000)==0
       disp(loop);
       disp(MaxErr);
    end
end

rdata=linspace(dr,dr*N,N);
zdata=linspace(dz,dz*M,M);
createfigure(rdata,zdata,phi');

function createfigure(xdata1,ydata1,zdata1)
    figure1 = figure; % figure を作成
    axes1 = axes('Parent',figure1);% axes を作成 
    hold(axes1,'on');
    % contour を作成
    [c1,h1] = contour(xdata1,ydata1,zdata1,'TextStep',0.1,'LevelStep',0.01);
    clabel(c1,h1);
    ylabel({'z方向'}); % ylabel を作成
    xlabel({'r方向'}); % xlabel を作成
    box(axes1,'on');
    axis(axes1,'tight');
    % 残りの座標軸プロパティの設定
    set(axes1,'BoxStyle','full','FontName','Times New Roman','FontSize',24,'Layer','top','XGrid','on','XMinorGrid','on','XTick',[5 10 15 20 25 30],'YGrid','on','YMinorGrid','on','YTick',[5 10 15 20 25 30]);
    colorbar(axes1); % colorbar を作成 
end