How to Solve system of nonlinear PDE

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

[Torsten]

I don't understand all the places in your code where you index arrays with "round_n".

Why don't you insert the actual values for the variable in gridpoint i ?

[Mohammed Hamza]

round_n refers to the round that ode45 is working on.

[Mohammed Hamza]

Also, if you can notice those values which depends on round_n, they are not under any type of loops other than ode45.

[Mohammed Hamza]

Also those values are dependent on distance/time so I can't keep them constant.

[Torsten]

Take away all these round_n's. They don't make sense. Replace the values by their actual values in the respective grid points.

[Mohammed Hamza]

I tried doing that, but the results looks unrealistic. Like how the amount of a material becomes in negative?

[Torsten]

Then your equations must be wrong.

For example, in the material balances, the sink term must be proportional to the material that is present in the grid point, but in your case, it's proportional to 1-(the material present in the grid point).

This will lead to negative concentrations for A and B.

[Mohammed Hamza]

All equations were taken from Chemical Reactor Design and Control. By William L. Luyben. It is hard to belive that the equations are wrong. Can you emphisis on the point of grid you have mentioned? Or can you copy the code edit the wrong protion, run it & then send it back?

[Torsten]

Do you have a copy of the equations from the book you could include ?

And why do you start with T = 0 K ?

[Mohammed Hamza]

Yes they are attached above in the question. I will reattach them again now.

[Mohammed Hamza]

Dear, I haven't started with T = 0 K. I have re-assigned a value for T @z=0 t=0 to be 445 so that as we have initial value & as we move the ode45 generate the next numerical values.

[Torsten]

T0 = zeros(n,1);
y0 = [ T0 ; Ya0 ; Yb0 ];

This will assign T=0 K over the complete z-range except T(1) which is set to 445 K.

It will take a while until T=445 reaches the end of the z-interval by convection ...

[Mohammed Hamza]

I still didn't get it. Can you elaborte ?

[Torsten]

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]

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.

Answer 1

[Torsten]

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

[Torsten]

Yes, I had put your script into a function and forgot to remove the "end" when copying it back.

Thanks for pointing it out.