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

65 views (last 30 days)
RITIKA Jaiswal on 3 Oct 2022
Commented: RITIKA Jaiswal on 9 Apr 2023
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.
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.

William Rose on 3 Oct 2022
Edited: William Rose on 4 Oct 2022
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)=...
@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
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');
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.
RITIKA Jaiswal on 4 Oct 2022
Edited: RITIKA Jaiswal on 4 Oct 2022
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.

RITIKA Jaiswal on 4 Oct 2022
Edited: RITIKA Jaiswal on 4 Oct 2022 @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?
RITIKA Jaiswal on 9 Apr 2023
@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.

### Categories

Find more on Boundary Conditions in Help Center and File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!