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

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

54 views (last 30 days)

Show older comments

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.

##### 6 Comments

William Rose
on 3 Oct 2022

Your equaiton 2:

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

and here

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.

### Answers (2)

William Rose
on 3 Oct 2022

Edited: William Rose
on 4 Oct 2022

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

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 Comments

RITIKA Jaiswal
on 4 Oct 2022

Edited: RITIKA Jaiswal
on 4 Oct 2022

##### 8 Comments

William Rose
on 4 Oct 2022

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.

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!