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

Question

RITIKA Jaiswal 2022-10-3

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1816490-how-to-solve-systems-of-non-linear-partial-differential-equations-using-finite-difference-method

评论： RITIKA Jaiswal 2023-4-9

1-s2.0-S0098135417300960-main.pdf

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

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

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

RITIKA Jaiswal 2022-10-3

编辑：RITIKA Jaiswal 2022-10-3

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

William Rose 2022-10-3

@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 2022-10-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1816490-how-to-solve-systems-of-non-linear-partial-differential-equations-using-finite-difference-method#answer_1065910

编辑：William Rose 2022-10-4

在 MATLAB Online 中打开

@RITIKA Jaiswal,

Edited:

Corrected error in the code in the equation for Tcool.
Changed vz=1 to vz=0.5, to match value given by @RITIKA Jaiswal in a comment.
Changed Nr to 26 and Nt to 6251, to match the ratios and , given in a comment.
Fixed error in formula for dTdr which caused run time error at line T(Nr,j,k)=...
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

%% 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.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2022-10-3

在 MATLAB Online 中打开

@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

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

William Rose 2022-10-3

@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.

Torsten 2022-10-3

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

RITIKA Jaiswal 2022-10-4

编辑：RITIKA Jaiswal 2022-10-4

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 2022-10-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1816490-how-to-solve-systems-of-non-linear-partial-differential-equations-using-finite-difference-method#answer_1066415

编辑：RITIKA Jaiswal 2022-10-4

@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?

8 个评论
显示 6更早的评论隐藏 6更早的评论

William Rose 2022-10-4

@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 2023-4-9

1-s2.0-S0098135417300960-main (2).pdf

@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.

请先登录，再进行评论。