How to Solve system of nonlinear PDE

3 次查看（过去 30 天）
Mohammed Hamza
Mohammed Hamza 2022-4-2
编辑： Torsten 2022-4-3
Hello,
I'm trying to solve this system of non-linear equations for a while. Unfortunatly it seems that the code doesn't work as requested. The code attached below is used to model a PFR system. Would be really thankful if anyone could help :)
The equations are attached above.
clc
clear all
close all
Tin = 445; % Feed concentration
L = 17; % Reactor length
t0 = 0; % Initial Time
tf = 120; % Final time
nt = 100; % Number of time steps
t = linspace(t0, tf, nt); % Time vector
time = t;
n = 10; % Number of axial steps
z = linspace(0,L,n); % Axial vector
n = numel(z); % Size of mesh grid
T0 = zeros(n,1); % t = 0, T = 0
T0(1) = Tin; % t = 0, T = 440 for z = 1
Ya0 = zeros(n,1); % t = 0, Ya = 0 for all z,
Ya0(1) = 0.5; % t = 0, Ya = 0.4 for z = 1
Yb0 = zeros(n,1); % t = 0, Yb = 0 for all z,
Yb0(1) = 0.5; % t = 0, Yb = 0.3 for z = 1
y0 = [ T0 ; Ya0 ; Yb0 ]; % All intial values enter to the same vector
% Appends conditions together
[t, y] = ode45(@(t,y) f(t,y,z,n,time,nt),[t0 tf],y0);
% t is the time
% y is T,Ya,Yb
% z is the axial vector
% n is the number of axial steps
% time is
% Plotting
figure; plot(z, y(end,n+1:2*n));
title('Ya at final time & z=1');
xlabel('distance')
ylabel('Ya')
figure; plot(z, y(end,2*n+1:3*n));
title('Yb at final time at final time & z=1');
xlabel('distance')
ylabel('Yb')
figure; plot(z, y(end,1:n));
title('T at final time at final time & z=1' );
xlabel('distance')
ylabel('Temperature')
function dydt=f(t,y,z,n,time,nt)
% Constant Parameters
D_p = 0.003;
mu = 0.18*(10^-4);
epsilon = 0.4;
alpha = 0.19038;
rho_cat = 2000;
lambda = 23237;
E = 69710;
R_rate = 8.314; %kJ/kmol.K
R_conc = 0.08314; % m^3.bar/kmol.K
MA = 15;
MB = 20;
c = 2;
D_r = 1.71; %Diameter of reactor
L_r = 17; %Length of reactor
c_cat = 0.5;
Area = pi*(D_r^2)/4;
Pressure = 50 ;
% Initiallizing the derivatives
dTdt = zeros(n,1);
dYadt = zeros(n,1);
dYbdt = zeros(n,1);
dTdz = zeros(n,1);
dYadz = zeros(n,1);
dYbdz = zeros(n,1);
dt = zeros(n,1);
% Extracting the initial values
T = y(1:n);
Ya = y(n+1:2*n);
Yb = y(2*n+1:3*n);
% Defining the axial change
%for i=2:n-1
for i=2:n
dTdz(i)= (T(i)-T(i-1))/(z(i)-z(i-1));
dYadz(i)= (Ya(i)-Ya(i-1))/(z(i)-z(i-1));
dYbdz(i)= (Yb(i)-Yb(i-1))/(z(i)-z(i-1));
end
% Calculated Parameters
round_n = nnz(dTdt)+1;
rho_molar = Pressure/(T(round_n)*1.01325*0.082057);
MM_avg = Ya(round_n)*MA+Yb(round_n)*MB;
rho = rho_molar*MM_avg;
V_rate = Pressure*MM_avg/rho;
v = V_rate/Area;
%Re = D_p*v*rho/mu;
%f = (1-epsilon)*(1.75+150*(1-epsilon)/Re)/(epsilon^3);
r_c = alpha*rho_cat*exp(-E/(R_rate* T(round_n)))*Ya(round_n)*Yb(round_n)*(Pressure)^2;
C = Pressure/(R_conc*T(round_n));
% Calculating outputs
for i=2:n
dTdt(i) = (-c*rho*v*(dTdz(i))-lambda*r_c)/(rho_cat*c_cat);
% T(i) = dTdt(i)*dt(i)+T(i-1);
dYadt(i) = (-v*dYadz(i)-(1-Ya(round_n))*r_c/C)/epsilon;
% Ya(i) = dYadt(i)*dt(i)+Ya(i-1);
dYbdt(i) = (-v*dYbdz(i)-(1-Yb(round_n))*r_c/C)/epsilon;
% Yb(i) = dYbdt(i)*dt(i)+Yb(i-1);
end
% Send the latest outputs back
dydt=[dTdt;dYadt;dYbdt];
end
  15 个评论
显示 13更早的评论
Torsten
Torsten 2022-4-2
编辑：Torsten 2022-4-2
You set an initial condition for T to 0 K over the complete z-range except for T(1) which is 445 K (same for the molar fractions).
T0 = zeros(n,1); % t = 0, T = 0
T0(1) = Tin; % t = 0, T = 440 for z = 1
Ya0 = zeros(n,1); % t = 0, Ya = 0 for all z,
Ya0(1) = 0.5; % t = 0, Ya = 0.4 for z = 1
Yb0 = zeros(n,1); % t = 0, Yb = 0 for all z,
Yb0(1) = 0.5; % t = 0, Yb = 0.3 for z = 1
y0 = [ T0 ; Ya0 ; Yb0 ]; % All intial values enter to the same vector
% Appends conditions together
[t, y] = ode45(@(t,y) f(t,y,z,n,time,nt),[t0 tf],y0);
Mohammed Hamza
Mohammed Hamza 2022-4-2
But that doesn't seem to change a lot. Still the numbers barely changes, and doesn't make sense. So, I think this won't solve the big problem yet.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-4-2
编辑：Torsten 2022-4-3
Looks better, doesn't it ?
Tin = 445; % Feed concentration
L = 17; % Reactor length
t0 = 0; % Initial Time
tf = 120; % Final time
nt = 100; % Number of time steps
t = linspace(t0, tf, nt); % Time vector
time = t;
n = 100; % Number of axial steps
z = linspace(0,L,n); % Axial vector
n = numel(z); % Size of mesh grid
T0 = 293.15*ones(n,1); % t = 0, T = 0
T0(1) = Tin; % t = 0, T = 440 for z = 1
Ya0 = 0.01*ones(n,1); % t = 0, Ya = 0 for all z,
Ya0(1) = 0.5; % t = 0, Ya = 0.4 for z = 1
Yb0 = 0.01*ones(n,1); % t = 0, Yb = 0 for all z,
Yb0(1) = 0.5; % t = 0, Yb = 0.3 for z = 1
y0 = [ T0 ; Ya0 ; Yb0 ]; % All intial values enter to the same vector
% Appends conditions together
[t, y] = ode15s(@(t,y) f(t,y,z,n,time,nt),[t0 tf],y0);
% t is the time
% y is T,Ya,Yb
% z is the axial vector
% n is the number of axial steps
% time is
% Plotting
figure; plot(z,[y(20,n+1:2*n);y(40,n+1:2*n);y(60,n+1:2*n);y(80,n+1:2*n);y(end,n+1:2*n)]);
title('Ya at final time & z=1');
xlabel('distance')
ylabel('Ya')
figure; plot(z,[y(20,2*n+1:3*n);y(40,2*n+1:3*n);y(60,2*n+1:3*n);y(80,2*n+1:3*n);y(end,2*n+1:3*n)]);
title('Yb at final time at final time & z=1');
xlabel('distance')
ylabel('Yb')
figure; plot(z,[y(20,1:n);y(40,1:n);y(60,1:n);y(80,1:n);y(end,1:n)]);
title('T at final time at final time & z=1' );
xlabel('distance')
ylabel('Temperature')
function dydt=f(t,y,z,n,time,nt)
% Constant Parameters
D_p = 0.003;
mu = 0.18*(10^-4);
epsilon = 0.4;
alpha = 0.19038;
rho_cat = 2000;
lambda = 23237;
E = 69710;
R_rate = 8.314; %kJ/kmol.K
R_conc = 0.08314; % m^3.bar/kmol.K
MA = 15;
MB = 20;
c = 2;
D_r = 1.71; %Diameter of reactor
L_r = 17; %Length of reactor
c_cat = 0.5;
Area = pi*(D_r^2)/4;
Pressure = 50 ;
% Initiallizing the derivatives
dTdt = zeros(n,1);
dYadt = zeros(n,1);
dYbdt = zeros(n,1);
dTdz = zeros(n,1);
dYadz = zeros(n,1);
dYbdz = zeros(n,1);
dt = zeros(n,1);
% Extracting the initial values
T = y(1:n);
Ya = y(n+1:2*n);
Yb = y(2*n+1:3*n);
% Defining the axial change
%for i=2:n-1
for i=2:n
dTdz(i)= (T(i)-T(i-1))/(z(i)-z(i-1));
dYadz(i)= (Ya(i)-Ya(i-1))/(z(i)-z(i-1));
dYbdz(i)= (Yb(i)-Yb(i-1))/(z(i)-z(i-1));
end
% Calculating outputs
for i=2:n
% Calculated Parameters
rho_molar = Pressure/(T(i)*1.01325*0.082057);
MM_avg = Ya(i)*MA+Yb(i)*MB;
rho = rho_molar*MM_avg;
V_rate = Pressure*MM_avg/rho;
v = V_rate/Area;
%Re = D_p*v*rho/mu;
%f = (1-epsilon)*(1.75+150*(1-epsilon)/Re)/(epsilon^3);
r_c = alpha*rho_cat*exp(-E/(R_rate* T(i)))*Ya(i)*Yb(i)*(Pressure)^2;
C = Pressure/(R_conc*T(i));
dTdt(i) = (-c*rho*v*dTdz(i)-lambda*r_c)/(rho_cat*c_cat);
% T(i) = dTdt(i)*dt(i)+T(i-1);
dYadt(i) = (-v*dYadz(i)-(1-Ya(i))*r_c/C)/epsilon;
% Ya(i) = dYadt(i)*dt(i)+Ya(i-1);
dYbdt(i) = (-v*dYbdz(i)-(1-Yb(i))*r_c/C)/epsilon;
% Yb(i) = dYbdt(i)*dt(i)+Yb(i-1);
end
% Send the latest outputs back
dydt=[dTdt;dYadt;dYbdt];
end
  1 个评论
Torsten
Torsten 2022-4-3
Yes, I had put your script into a function and forgot to remove the "end" when copying it back.
Thanks for pointing it out.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by