How to solve systems of non linear partial differential equations using finite difference method?

I am trying to solve Sets of pdes in order to get discretize it.Using finite difference method such that the resulting ODEs approximate the essential dynamic information of the system.

I am not sure whether the code is correct. I have used first order forward difference and 2nd order centered difference.

i am unclable to solv equation (2).Please guide.

clear all;

close all;

M=1000;

c=0.25;%lets dt/dr^2 =c

a=0.02;%lets dt/dr=a

r=0.01;

v=0.5;

for i =2:25

for j =2:25

p(i,j)=200;

end

end

dt=0.001;

dr=0.25;

for t=1:M

for i=2:25

for j =2:24

pp(i,j)=p(i,j)+(0.5*a)*(-v+(1/r))*(p(i,j+1)-p(i,j-1))+c*(p(i,j+1)-2*p(i,j)+p(i,j-1));

end

end

n=25;

pp(i,1:n)=500; %lets assume

pp(n,1:n)=500;

pp(1:n,1)=500;

pp(1:n,n)=500;

p=pp;

t=t;

end

figure

contourf(p,25,'linecolor','non')

this code is running i have also attached paper here for reference.

[Torsten]

pp(i,j)=p(i,j)+(0.5*a)*(-v+(1/r))*(p(i,j+1)-p(i,j-1))+c*(p(i,j+1)-2*p(i,j)+p(i,j-1));

1. r depends on j !
2. (-v+(1/r))*(p(i,j+1)-p(i,j-1)) is wrong: v refers to dp/dz, not dp/dr !

[RITIKA Jaiswal]

yeah some parts i have assumed constants as i dont know how to write that

[RITIKA Jaiswal]

I have used finite difference method in order to discretize the equation (Dr,a,r,Vz) here i am struggling. how to use the boundary condition?

[Torsten]

I suggest

http://www.scholarpedia.org/article/Method_of_lines

for an introduction to the discretization of convection-diffusion equations in two spatial dimensions.

If you need some further literature, don't hesitate to ask.

[RITIKA Jaiswal]

i used finite difference method to discretize it. and final equation will be like this

[William Rose]

@RITIKA Jaiswal,

Your equaiton 2:

This is not really an answer but a suggesiton: see here

https://www.mathworks.com/matlabcentral/answers/1775970-1d-heat-conduction-using-eulers-explicit-discretisation

and here

https://www.mathworks.com/matlabcentral/answers/1641800-explicit-heat-conduction-in-3d

for examples of solving a partial differential equation for heat.

The two examples above are considerably simpler than the example you posted. In your problem:

• There are 7 variables to solve for: 6 gases plus temperature.
• The 6 PDEs for gases are relatively sraightforward. Each gas partial differential equaiton is independent of the other gases and they are all independent of temperature. Therefore they can be solved independently, before solving for temperature.
• The temperature boundary conditions are given for the temperature gradient at r=0 and r=R.
• The temperature PDE involves the total gass density (sum of 6 gases) at each point and time. Therefore you solve for the gas densities first. Then you compute the 3D array of total gas density at every point and every time. Then you work on the temperature PDE.
• The arrays are 3D. The dimensions are r=radius, z=axial position, t=time.

Answer 1

[William Rose]

@RITIKA Jaiswal,

Edited:

1. Corrected error in the code in the equation for Tcool.
2. Changed vz=1 to vz=0.5, to match value given by @RITIKA Jaiswal in a comment.
3. Changed Nr to 26 and Nt to 6251, to match the ratios and , given in a comment.
4. Fixed error in formula for dTdr which caused run time error at line T(Nr,j,k)=...
5. Added plots.

@Torsten gives excellent suggesitons.

Let us define 3D arrays for rho=total gas density and T=temperature.

Equation 2 in the paper:

Discretized version of equation 2 :

i=radius, j=z (axial coordinate), k=time. The right hand side involves only T(..,..,k-1), and never T(..,..,k). That means you can step through time, computing T(i,j,k) using only the T values for previous times. In the code below, I compute T at all radii except the centerline and the outside edge (r=R). Then I use the boundary conditions to compute T at r=0 and r=R.

Nr=26; Nz=101; Nt=6251;                         %number of values for radius, z, time
R=2; L=8; tmax=10;                              %max values along each axis
dr=R/(Nr-1); dz=L/(Nz-1); dt=tmax/(Nt-1);       %deltas
r=(0:Nr-1)*dr; z=(0:Nz-1)*dz; t=(0:Nt-1)*dt;    %vectors for r, z, t
rho=zeros(Nr,Nz,Nt);  T=zeros(Nr,Nz,Nt);        %allocate arrays
rhoInlet=[1,1,1,1,1,1];                         %inlet density for each gas
Lr=1; cp=1; vz=0.5; kw=1;                       %constants
Tin=400; Tcoolin=273; Tcoolout=350;             %constants
k1=10;                                          %constant for boundary condition on T at r=R
%Compute rho1, rho2, rho3, rho4, rho5, rho6.
rho1=rhoInlet(1)*ones(Nr,Nz,Nt);                %actual calculation will be complicated
rho2=rhoInlet(2)*ones(Nr,Nz,Nt);
%etc
rho=rho1+rho2;                                  %actual rho=rho1+...+rho6;
%After we compute rho, we can compute T.
%First, compute Tcool(z)
Tcool=Tcoolin*(1.-z/L)+Tcoolout*z/L;
%Now compute T(i,j,k)
T(:,:,1)=Tin*ones(Nr,Nz);       %temperature distribution at t=0 - is correct?
T(:,1,:)=Tin*ones(Nr,1,Nt);     %temperature distribution at z=0
for k=2:Nt
    for j=2:Nz
        for i=2:Nr-1
            T(i,j,k)=T(i,j,k-1)+dt*...
                      ( (-vz/dz)*(T(i,j,k-1)-T(i,j-1,k-1))...
                       +(Lr/(rho(i,j,k)*cp))*( (1/dr^2)*(T(i+1,j,k-1)-2*T(i,j,k-1)+T(i-1,j,k-1)) ...
                                              +(1/(r(i)*dr))*(T(i,j,k-1)-T(i-1,j,k-1))));
        end
        T(1,j,k)=T(2,j,k);                          %boundary condition at r=0: dT/dr=0
        dTdr=-k1+(kw/Lr)*(Tcool(j)-T(Nr,j,k-1));    %boundary condition at r=R
        T(Nr,j,k)=T(Nr-1,j,k)+dTdr*dr;              %compute T(r=R)
    end
end
%% Plot results
idxR=[1,6,11,16,21,26];         %radial indices for plotting
idxZ=[1,26,51,76,89,101];       %axial indices for plotting
idxT=[100,626,1251,2501,6251];  %time indices for plotting
Tlimits=[min(min(min(T))),max(max(max(T)))];
%Plot 1
subplot(121);
for k=1:length(idxT)
    plot(z,T(idxR(4),:,idxT(k)))
    legstr1{k}=sprintf('T=%.2f',t(idxT(k)));
    hold on
end
xline(z(idxZ(5)),'--r');
ylabel('Temperature'), xlabel('Z=Axial Position');
legend(legstr1); grid on
ylim(Tlimits);
titlestr=sprintf('Temp. vs Z, r=%.1f',r(idxR(4))); 
title(titlestr)
%Plot 2
subplot(122);
for k=1:length(idxT)
    plot(r,T(:,idxZ(5),idxT(k)))
    legstr2{k}=sprintf('T=%.2f',t(idxT(k)));
    hold on
end
xline(r(idxR(4)),'--r');
ylabel('Temperature'), xlabel('Radius');
legend(legstr2); grid on
ylim(Tlimits);
titlestr=sprintf('Temp. vs R, z=%.1f',z(idxZ(5))); 
title(titlestr)

The code runs without error, which is a good start. I do not promise that the formulas above are correct. I think they are. If it were my project, I would do a lot of error checking. Check the parentheses carefully, in case I made a typographical error.

You will need to put in correct values for all the constants. .

You will need to calculate the correct value for k1 using the first term on the right side of equation 4 in the paper.

I have made Tcool() a function of z only: Tcool(z)=linear interpolation from inlet to outlet, as suggested in the paper. The paper says Tcool() is a function t as well as z, but the portion of the paper shown above does not give a formula for Tcool,in(t) or Tcool,out(t), so I assume constant values for Tcool,in and Tcool,out.

rho1 ..., rho6, should be computed using the respective PDEs, before computing temperature. I made them constant 3D arrays, for simplicity. After you compute rho1, ..., rho6, then compute rho=sum(rho1,...,rho6).

I have assumed T(z=0)=Tin, in accordance with equation 4 in the paper. I have also assumed T(t=0)=Tin. The portion of the paper shown above does not specifiy the initial temparature distribution.

[Torsten]

@William Rose

Better use

Nt = 20001 

instead of

Nt = 201

and

dTdr=-k1+(kw/Lr)*(Tcool(j)-T(Nr,j,k-1));   %boundary condition at r=R

instead of

dTdr=-k1+(kw/Lr)*(Tcool-T(Nr,j,k-1));   %boundary condition at r=R

to see reasonable results in a final plot command for the axial temperature profile:

plot(z,[T(10,:,10000);T(10,:,15000);T(10,:,20000)])

Nr=21; Nz=161; Nt=20001;                        %number of values for radius, z, time
R=1; L=16; tmax=10;                             %max values along each axis
dr=R/(Nr-1); dz=L/(Nz-1); dt=tmax/(Nt-1);       %deltas
r=(0:Nr-1)*dr; z=(0:Nz-1)*dz; t=(0:Nt-1)*dt;    %vectors for r, z, t
rho=zeros(Nr,Nz,Nt);  T=zeros(Nr,Nz,Nt);        %allocate arrays
rhoInlet=[1,1,1,1,1,1];                         %inlet density for each gas
Lr=1; cp=1; vz=1; kw=1;                         %constants
Tin=400; Tcoolin=273; Tcoolout=350;             %constants
k1=10;                                          %constant for boundary condition on T at r=R
%Compute rho1, rho2, rho3, rho4, rho5, rho6.
rho1=rhoInlet(1)*ones(Nr,Nz,Nt);                %actual calculation will be complicated
rho2=rhoInlet(2)*ones(Nr,Nz,Nt);
%etc
rho=rho1+rho2;                                  %actual rho=rho1+...+rho6;
%After we compute rho, we can compute T.
%First, compute Tcool(z)
Tcool=Tcoolin*(1.-z)/L+Tcoolout*z/L;
%Now compute T(i,j,k)
T(:,:,1)=Tin*ones(Nr,Nz);       %temperature distribution at t=0 - is correct?
T(:,1,:)=Tin*ones(Nr,1,Nt);     %temperature distribution at z=0
for k=2:Nt
    for j=2:Nz
        for i=2:Nr-1
            T(i,j,k)=T(i,j,k-1)+dt*...
                      ( (-vz/dz)*(T(i,j,k-1)-T(i,j-1,k-1))...
                       +(Lr/(rho(i,j,k)*cp))*( (1/dr^2)*(T(i+1,j,k-1)-2*T(i,j,k-1)+T(i-1,j,k-1)) ...
                                              +(1/(r(i)*dr))*(T(i,j,k-1)-T(i-1,j,k-1))));
        end
        T(1,j,k)=T(2,j,k);                      %boundary condition at r=0: dT/dr=0
        dTdr=-k1+(kw/Lr)*(Tcool(j)-T(Nr,j,k-1));   %boundary condition at r=R
        T(Nr,j,k)=T(Nr-1,j,k)+dTdr*dr;          %compute T(r=R)
    end
end
plot(z,[T(10,:,10000);T(10,:,15000);T(10,:,20000)])

[William Rose]

@RITIKA Jaiswal,

I assumed that you knew how to solve the PDEs for rho1 to rho6. Maybe I was wrong about that. The solution I gave above for T(r,z,t) should give you guidance.

is, in a way, simpler than T(r,z,t), since , whereas .

The outer loop should be time. Compute rho(r,z,t=2) (i.e. rho at all r and z values), based on the time 1 values. Then rho(r,z,t=3), using rho(r,z,t=2). And so on.

Boundary conditrions:

For rho and T, the boundary conditions at r=0 are and . This makes sense since we are using cylindrical coordinates. The profiles of rho and T are expected to smooth across the centerline, which means they must be flat at r=0. To implement this numerically: At each time step, compute rho and T up to but not including the centerline value (using the values at the previous time step, including the value at the centerline in the previuos time step, which you will already know). Then make the centerline value equal to the value immediately adjacent to the centerline. This will satisfy the boundary condition. Likewise for T.

For T, the boundary condition at r=R is , where k1 i a somewhat complicated value that appears on the right hand side of equation 4 in the paper. To implement this numerically: At each time step, compute T up to but not including r=R, using the values at the previous time step, which you already know. Then compute using the formula above, where, for T in the formula, you use the value for T(r=R) from the previous time step. Then compute T(r=R) for the current time step by T(R)=T(R-dr)+dr*dT/dr. In other words, use the the current T value at r=R-dr, plus dT/dr, which you know from the boundary condition equation, to estimate T at r=R.

For rho1,...,rho6, the boundary condition at r=R is , where k3 is a somewhat complicated value that appears on the right hand side of equation 3 in the paper. To implement this numerically: At each time step, compute rho up to but not including r=R, using the values at the previous time step, which you already know. Then compute . Then compute rho(r=R) for the current time step by . In other words, use the the current value for at r=R-dr, plus , which you know from the boundary condition equation, to estimate at r=R.

[William Rose]

@Torsten, thank you. I have made the sggestions you suggested in the code I posted, and some other fixes, as noted now in my edited answer above.

[Torsten]

This was a lot of effort you invested. My appreciation !

[RITIKA Jaiswal]

Thankyou @Torsten @William Rose for putting so much effort and dedication.I really appricate to both of your hardwork.I am going through the solution and suggestions which you have given and trying to solve each step .I will keep updating about my approach and steps.

Answer 2

[RITIKA Jaiswal]

@Torsten @William RoseI read the paper found that we dont need to discretize in time domain .To get the in the form of equation 5 we need to discretize in space domain.In the paper within one finite volume (FV), all physicalcoefficients (e.g., diffusion coefficient Dr,˛, thermal conductiv-ity r) are assumed to be constant. How to get in the form of equation 5?

[Torsten]

Did you read my link:

http://www.scholarpedia.org/article/Method_of_lines

?

If you had, you would know.

Tdot(i,j) = ( (-vz/dz)*(T(i,j)-T(i,j-1))...
             +(Lr/(rho(i,j)*cp))*( (1/dr^2)*(T(i+1,j)-2*T(i,j)+T(i-1,j)) ...
             +(1/(r(i)*dr))*(T(i,j)-T(i-1,j))));

is the system of ordinary differential equations you have to solve (using ODE15S, e.g.).

[RITIKA Jaiswal]

yeah read it trying to implement on equation 1 first then 2nd.

[RITIKA Jaiswal]

It should be like this right? in case of space variable r, there will be change in j right?Please correct me if i am wrong

[William Rose]

@RITIKA Jaiswal,

I see that the full paper is here:

https://pure.mpg.de/rest/items/item_2417133_9/component/file_2582706/content

Do you aim to reproduce the work in that paper? In equation 5, x and xdot are vectors of length 4375. A1, A2 are matrices of size 4375x4375. B1 is 4375x175, and B2 is 4375x2. I expect that A1 and A2 will be mainly tridiagonal, i.e. they will have non-zero elements along the main diagonal and the diagonal lines just above and below the main diagonal. This is because each volume element's time evolution is determined by itself and its immediate neighbors only. Even with this large reduction in the number of non-zero matrix elements, it is still non-trivial.

[RITIKA Jaiswal]

yeah i am trying to solve these pdes equation given in the https://pure.mpg.de/rest/items/item_2417133_9/component/file_2582706/content

my aim is to first use the finite difference method to discretize the equations and then solve it to get in the form equation 5.

instead of finite volume method i am using finite difference method

[RITIKA Jaiswal]

@william if you know do help.

[William Rose]

@RITIKA Jaiswal,

It is not obvious to me exactly what the formulas are for the elements of A1, A2, B1, and B2, in equation 5 of Bremer et al (2017). That would take further thought. I would have to read and deeply understand the Bremer paper, which also appears in published version here. I do not have the time to do so.

The code I provided demonstrates how you can solve the PDE for temperature using an explicit approach. With a few modifications, it can also be used to solve for the six gas concentrations. The method of lines approach, recommended by @Torsten, uses an implicit approach. Advantages of the implicit approach are stability and, related to that, computational efficiency: it will not diverge, the way the explicit method can. The disadvantage of the MOL is the extra calculations you must do to derive the implicit equations. If you choose a small enough value for , the explicit approach will work, but such a small value of can lead to long computation times.

[RITIKA Jaiswal]

@Torsten Hi, Suppose I have ode form on the LHS side and it is in the form of equation 5. I have derived the equation 5 I have used finite volume method. Earlier i used finite difference method but i did mistake because in LHS side i was doing discretisation .Now i have converted Both two equations 1 and 2 to get equation 5. Now in order to solve ode of size 4375 I am unable to solve it and code it on the matlab.