How to Solve system of nonlinear PDE
조회 수: 8 (최근 30일)
이전 댓글 표시
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
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);
채택된 답변
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
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개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!