- 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?
40 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1143640/image.jpeg)
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 comentarios
William Rose
el 3 de Oct. de 2022
Your equaiton 2: ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1144185/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1144185/image.png)
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.
Respuestas (2)
William Rose
el 3 de Oct. de 2022
Editada: William Rose
el 4 de Oct. de 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: ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1144330/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1144330/image.png)
Discretized version of equation 2 :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1144335/image.jpeg)
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.
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1144555/image.png)
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 comentarios
RITIKA Jaiswal
el 4 de Oct. de 2022
Editada: RITIKA Jaiswal
el 4 de Oct. de 2022
8 comentarios
William Rose
el 4 de Oct. de 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1145570/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1145575/image.png)
Ver también
Categorías
Más información sobre Eigenvalue Problems en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!