کد متلب روش های عددی، مثال های کتاب سادیکو، فصل 3، مثال 3-9

% RJA - 2/23/07, revised 4/14/07

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MATLAB CODE FOR EXAMPLE 3.9:

% AN AXISYMMETRIC PROBLEM OF AN EARTHED CYLINDER PARTIALLY FILLED WITH 

% CHARGED LIQUID SOLVED USING FINITE DIFFERENCE SCHEME

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A=1; % Radius of cylindrical tank (meters)

B=1; % Height of liquid in tank

C=1; % Height of gas in tank

ER1=1;ER2=2;EO=8.854E-12; %dielectric parameters

H = 0.05;   %spacial step size

NA = A/H;   %number of points along A

NB = B/H;   %number of points along B                    

NC = C/H;   %number of points along C                    

NBC = NB + NC;  %%number of points along B & C

NMAX = 500;     %number of iterations

RHOV = 1E-5;                       

G = -RHOV/(ER2*EO);                

GH2 = G*H^2;                       

%INITIALIZE - THIS ALSO TAKES CARE OF DIRICHLET CONDITIONS

V = zeros(NA+1,NBC+1);   %V(radius, height) 

%NOW, APPLY FINITE DIFFERENCE SCHEME

for N = 1:NMAX            

    for I = 2:NA %step through radius       

        FM = (2*(I-1) - 1)/(2*(I-1)); %note that indicie starts at 1 instead of 0

        FP = (2*(I-1) + 1)/(2*(I-1));

        

        %step through liquid (z+ direction)        

        for J = 2:NB      

            V(I,J) = 0.25*( V(I,J-1) + V(I,J+1). 

                + FM*V(I-1,J) + FP*V(I+1,J) - GH2 );

        end                                                                     

        

        %step through gas (z+ direction)

        for J = NB+2:NBC                             

            V(I,J) = 0.25*( V(I,J-1) + V(I,J+1).

                + FM*V(I-1,J) + FP*V(I+1,J) );

        end                                                   

        %ALONG THE GASS-LIQUID INTERFACE

        V(I,NB+1) = 0.5*( V(I,NB+2)*ER1/(ER1+ER2)+.   

             V(I,NB)*ER2/(ER1+ER2) ).

             + 0.25*(FM*V(I-1,NB+1)+FP*V(I+1,NB+1) );

    end                                                       

    %IMPOSE NEUMANN CONDITION ALONG THE Z-AXIS

    for J = 2:NBC                                           

        V(1,J) = ( 4.0*V(2,J) + V(1,J-1)+V(1,J+1))/6.0;

    end                                                                

end                                                                    

% OUTPUT THE POTENTIAL ALONG RHO = 0.5, 0 < Z < 1.0

NR5 = 0.5/H;                                                          

                                                                           

radius = 0:H:A;

height = 0:H:(B+C);                                                                           

                                                                           

figure(1), 

    subplot(211),

        plot(height,V(NR5+1,(0:NBC)+1)/1000.0)

            ylabel('V(kV)')

            title('Potential distribution in the tank of Fig. 3.38 along \rho=0.5m')

    subplot(212)                                                                            

        plot(radius ,V((0:NA)+1,NB+1)/1000.0)

            title('Potential distribution in the tank of Fig. 3.38 along gas-liquid interface')

            ylabel('V(kV)')  

                  

    figure(2),surf(height,radius,V/1000)

        ylabel('radius in tank')

        xlabel('height in tank')

        zlabel('Potential (kV)')