Help with Solving PDE System with DAE in MATLAB

I hope you're doing well. My name is William. I’m currently working on a problem and was hoping you might be able to help. I have a system of PDEs, but one of the variables does not include a time derivative (e.g., no ∂u/∂t), even though the variable depends on both time and space. I'm not sure how to approach solving this in MATLAB. I’ve been able to implement it in gPROMS, but translating it into MATLAB has been challenging.

I tried to follow your DAE method, but my variables depend on both time and space—not just time—so I’m unsure how to proceed. If you could offer any guidance or point me in the right direction, I would greatly appreciate it.

Thank you for your time and support.

Best regards,

William

[Torsten]

We need to know your PDE system together with the DAE and the solution method you try to implement. And the MATLAB code so far - if you already started.

[Torsten]

Define a mass matrix M for your resulting DAE system.

It then takes the form

M*dy/dt = f(t,y).

Set the rows in M to 0 that correspond to the algebraic variables (velocity) and define the right-hand side as "velocity - value".

For an example, see

Solve Robertson Problem as Semi-Explicit Differential Algebraic Equations (DAEs)

under

https://uk.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html

I wonder what the algebraic equation for the velocity is and if it's not possible to derive its values directly from the other variables without defining it as an extra solution variable.

Don't you have a pdf file with a mathematical description of your model ? It's laborious to deduce everything just from your coding.

Answer 1

[Torsten]

You have discretized your equations in space. Thus the y variables are only functions of time for fixed z-coordinates.

Here is an example where y1 and y2 have time derivative:

dy1/dt = z + 2 , y1(z,0) = 0
dy2/dt = z - 4 , y2(z,0) = 0

and the third solution variable y3 is given by an algebraic equation

y3 = z + t.

tspan = 0:0.1:1;
z = (0:0.1:1).';
n = numel(z);
fun = @(t,y) [z+2;z-4;z+t-y(2*n+1:3*n)];
y0 = [zeros(n,1);zeros(n,1);zeros(n,1)];
M = [eye(2*n),zeros(2*n,n);zeros(n,2*n),zeros(n)];
options = odeset('Mass',M);
[T,Y] = ode15s(fun,tspan,y0,options);
figure(1)
plot(T,Y(:,1:n))   % should be y1 = t*(z+2)

figure(2)
plot(T,Y(:,n+1:2*n))  % should be y2 = t*(z-4)

figure(3)
plot(T,Y(:,2*n+1:3*n))  % should be y3 = z + t

[William]

Thank you very much for sharing the code. I’ll go through it in MATLAB to understand it in more detail.

I have a question regarding the initial conditions. In your example, you’ve provided initial conditions for y3. However, in my case, I don’t have initial conditions for u. In the gPROMS code, there’s no initial condition specified for u—the value I used in MATLAB was something I assigned manually.

From what I understand, u might be calculated from P and y, but I’m still quite confused about how that works. I’d really appreciate your guidance on this.

Thanks again!

[Torsten]

If you don't supply correct initial values for the algebraic variables, MATLAB tries to determine them consistently. Sometimes this works, sometimes not. In the above code, I changed the y0 vector from

y0 = [zeros(n,1);zeros(n,1);z];

to

y0 = [zeros(n,1);zeros(n,1);zeros(n,1)];

and it still works.

Did your code work in gProms ? Or why did you switch to MATLAB ?

Why don't you solve the Ergun equation for u in function f (maybe iteratively) ? This way, you don't have u as solution variable, but only as a deduced quantity.

[William]

The original code runs well in gPROMS. However, my current workplace does not have a gPROMS license—we primarily use MATLAB—so I’m working on converting the code to MATLAB.

I’ve also tried implementing the model in Aspen Custom Modeler (ACM). While it runs, the performance is significantly slower: gPROMS solves everything in under 5 seconds, whereas ACM takes 1–2 minutes. When I increase the number of discretization points (N) and reduce the timestep, gPROMS still finishes in under a minute, but ACM slows down drastically, taking 5–10 minutes.

Given these limitations, I’m now trying to implement the model in MATLAB to see if it can improve the speed and meet the requirements of my workplace.

I attempted to solve for u = f(P1, ...) using the Ergun equation, but the code performs poorly—it runs extremely slowly and fails to converge, even in the first few time steps (0.0xx seconds).

Thank you again for taking the time to support me. I really appreciate your help!

[Torsten]

 0== -dP1dz(z)-150 * mu .* u(z) .* ((1-void)^2) / (4 * 2.46e-6 * void^3)- 1.75 * (1-void) .* u(z) .* (abs(u(z))) .* Rho_gas(z) ./ (2 * 1.57e-3 * void^3) ;

In one spatial dimension, it's almost impossible that u changes sign over the length of the reactor. Thus u*abs(u) is either u^2 or -u^2 throughout depending on whether the flow is from the left or from the right. So it's a quadratic equation in u that has to be solved. Nothing complicated in my opinion.

[William]

Thank you so much for your continued support.

I actually tried using the method you suggested, but unfortunately the results still seem incorrect. Specifically, after solving, the values of y2 remain nearly flat around 0.14, whereas in both gPROMS and ACM, y2 gradually increases over time.

I’ve attached the code below for your reference. I spent quite a bit of time reviewing it, but I couldn’t identify any issues.

In addition, I included comparison plots of q from both MATLAB and gPROMS. In gPROMS, the two q curves intersect around 80 seconds, while in MATLAB, they cross much earlier—around 30 seconds. You can also observe that the rate of change for q2 is significantly different between the two.

I also compared the P values: the results from MATLAB are not consistent, whereas the values from gPROMS and ACM align closely. In MATLAB, P is much higher, which leads to higher concentrations and a faster increase in q—which doesn’t match the behavior in the other platforms.

I’m not sure where the issue lies in the MATLAB code, and I’d really appreciate it if you could help me take a look.

clc; clear; close all;

% Initial variables

L=1;

R_bed=1.1e-2;

k11=-9.471;

k12=9.158;

k21=0.05369;

k22=0.01731;

k31=6.273e-6;

k32=1.167e-6;

k41=2595.6;

k42=3879.16;

k51=1;

k52=1;

k61=0;

k62=0;

H_Ads1=21581.1;

H_Ads2=32253.3;

De1=9.646e-8;

De2=8.109e-9;

Dp=3.14e-3;

R_p=1.57e-3;

m1=0.614786238793;

m2=0.051682579415;

y_feed1=0.85;

y_feed2=0.15;

R=8.3145;

Rho_p=1.16e3;

void=0.315;

tot_void=0.76;

Cp_gas=0.977e3;

Rho_b=0.795e3;

Cp_s=9.24e2;

ke=0.09;

Ui=60;

mu=1.75e-5;

Patm=1e5;

BedArea=3.14*R_bed^2;

Rho_wall=7800;

Cp_wall=504;

A_wall=4.24e-4;

R_bed_outside=R_bed+5e-3;

h_outside_air=200;

gamma=1.4;

Pfeed=10000;

Ppurge=100000;

T_ads=298;

T_des=373;

u_feed=0.20;

LLength=1.0;

PressTime=100;

%=========== Meshing ===========

tf = 100;

N=101;

Nt=101;

z = linspace(0, L, N);

dz = z(2) - z(1);

t = linspace(0, tf, Nt);

%======Initial Conditions========

ICT=T_des*ones(1,N);

ICP=Ppurge*ones(1,N);

ICT_Wall=T_des*ones(1,N);

ICq1=zeros(1,N);

ICq2=zeros(1,N);

ICy1= ones(1,N);

ICy2= zeros(1,N);

ICy1(1)=y_feed1;ICy2(1)=y_feed2;

IC=[ICT,ICP,ICT_Wall,ICq1,ICq2,ICy1,ICy2];

%options= odeset('RelTol',1e-12,'AbsTol',1e-12);

[t, y]=ode15s(@f,t,IC,[],void, A_wall, Cp_gas, Rho_p, Rho_b, Rho_wall, H_Ads1, H_Ads2, h_outside_air, Ui, mu,T_ads,T_des,PressTime, Pfeed,R, R_bed, R_bed_outside, tot_void, Cp_s, Cp_wall,dz,k11,k12, k21,k22, k31,k32, k41,k42, k51,k52, k61, k62, u_feed, y_feed1, y_feed2, m1,m2, R_p, N,LLength);

y1=y(:,1+5*N:6*N);y2=y(:,1+6*N:7*N);y1(:,1)=y_feed1;y2(:,1)=y_feed2; y1(:,end)=y1(:,end-1);y2(:,end)=y2(:,end-1);

q1=y(:,1+3*N:4*N);q2=y(:,1+4*N:5*N);

T=y(:,1:N);

T(:,end)=T(:,end-1);

P=y(:,N+1:2*N); P(:,end)=Pfeed;P(:,1)=P(:,2);

figure (1)

plot(t,q1(:,end),'b','linewidth',2);

hold on;grid on;

plot(t,q2(:,end),'r','linewidth',2)

xlabel('Time t (s)'); ylabel('Value of q');

function dydt=f(t,y,void, A_wall, Cp_gas, Rho_p, Rho_b, Rho_wall, H_Ads1, H_Ads2, h_outside_air, Ui, mu,T_ads,T_des,PressTime, Pfeed,R, R_bed, R_bed_outside, tot_void, Cp_s, Cp_wall,dz,k11,k12, k21,k22, k31,k32, k41,k42, k51,k52, k61, k62, u_feed, y_feed1, y_feed2, m1,m2, R_p,N ,LLength)

dydt=zeros(length(y),1);

dTdt=zeros(N,1);

dPdt=zeros(N,1);

dT_Walldt=zeros(N,1);

dq1dt=zeros(N,1);

dq2dt=zeros(N,1);

dy1dt = zeros(N,1);

dy2dt = zeros(N,1);

T=y(1:N);

P=y(N+1:2*N);

T_Wall=y(2*N+1:3*N);

q1=y(1+3*N:4*N);

q2=y(1+4*N:5*N);

y1=y(1+5*N:6*N);

y2=y(1+6*N:7*N);

for z = 2:(N) %*BFDM1 Method*%

Px(z)=(P(z)-P(z-1))/(dz);

Tx(z)=(T(z)-T(z-1))/(dz);

y1x(z)=(y1(z)-y1(z-1))/(dz);

y2x(z)=(y2(z)-y2(z-1))/(dz);

end

%Px(1)=(-P(1)+P(2))/(dz);

%Tx(1)=(-T(1)+T(2))/(dz);

%y1x(1)=(-y1(1)+y1(2))/(dz);

%y2x(1)=(-y2(1)+y2(2))/(dz);

%/*Boundary Condition*/%

dTdt(1)=(T_ads - T_des)/PressTime;

T(N)=T(N-1);

P(N)=Pfeed;

P(1)=P(2);

y1(1)=y_feed1;

y1(N)=y1(N-1);

y2(1)=y_feed2;

y2(N)=y2(N-1);

%/*Governent Equation*/

for z = 1:N

P_i1(z)=max(min(y1(z),1),1e-6)*(abs(P(z))+1e-20)/1e5;

P_i2(z)=max(min(y2(z),1),1e-6)*(abs(P(z))+1e-20)/1e5;

Rho_gas(z)=P(z)*(max(min((y1(z)),1),1e-6)*0.028+max(min(y2(z),1),1e-6)*0.0440)/R/(T(z));

end

for z=1:N

P1(z)=P(z);

end

for z = 2:(N-1) %*FFDM1 Method*%

P1x(z)=(-P(z)+P(z+1))/(dz);

end

P1x(N)=(P(N)-Px(N-1))/(dz);

%/*Ergun Equation*/

for z=1:N

u(1)=u_feed;

for z = 2:N-1

AA(z)=1.75*Rho_gas(z)*(1-void)./(2*1.57e-3*void^3);

BB=(150*mu*(1-void)^2)/(4*2.46e-6*void^3);

CC(z)= P1x(z)/LLength;

u(z)=(-BB + sqrt(BB^2-4.*AA(z).*CC(z)))./(2.*AA(z));

ux(z)=(u(z)-u(z-1))/(dz);

end

u(N)=u(N-1);

end

%===============Langmuir Equations======%

for z = 1: N

qm1(z)=k11+k21*T(z);

B1(z)=k31*exp((abs(k41)+1e-20)/(T(z)));

n1(z)=k51+k61/(T(z));

qm2(z)=k12+k22*T(z);

B2(z)=k32*exp((abs(k42)+1e-20)/(T(z)));

n2(z)=k52+k62/(T(z));

q_star1(z)=(qm1(z)*B1(z)*(abs(P_i1(z))+1e-20)^n1(z))/(1+(B1(z)*(abs(P_i1(z))+1e-20)^n1(z)+B2(z)*(abs(P_i2(z))+1e-20)^n2(z)));

q_star2(z)=(qm2(z)*B2(z)*(abs(P_i2(z))+1e-20)^n2(z))/(1+(B1(z)*(abs(P_i1(z))+1e-20)^n1(z)+B2(z)*(abs(P_i2(z))+1e-20)^n2(z)));

end

for z=1:N

%/*Linear Driving Force*/

dq1dt(z)=m1*(q_star1(z)-q1(z));

dq2dt(z)=m2*(q_star2(z)-q2(z));

end

%===========Component Mass Balance================%

for z = 2:N-1

%dy1dt(z) = - (max(u(z),1e-6)*y1x(z)/LLength/void+(1-void)/void*(R.*T(z)./P(z))*Rho_p*(dq1dt(z)-max(min(y1(z),1),1e-6)*(dq1dt(z)+dq2dt(z))));

%dy2dt(z) = - (max(u(z),1e-6)*y2x(z)/LLength/void+(1-void)/void*(R.*T(z)./P(z))*Rho_p*(dq2dt(z)-max(min(y2(z),1),1e-6)*(dq1dt(z))+dq2dt(z)));

dy1dt(z) =- (max(u(z),1e-6)*y1x(z)/LLength/void+(1-void)/void*(R*T(z)/P(z))*Rho_p*(dq1dt(z)-max(min(y1(z),1),1e-6)*(dq1dt(z)+dq2dt(z))));

dy2dt(z) =- (max(u(z),1e-6)*y2x(z)/LLength/void+(1-void)/void*(R*T(z)/P(z))*Rho_p*(dq2dt(z)-max(min(y2(z),1),1e-6)*(dq1dt(z)+dq2dt(z))));

end

%==========Overal Mass Balance============%

for z=2:N-1

%dPdt(z) =P(z)*(-ux(z)/LLength -u(z)/P(z)*Px(z)/LLength/void+u(z)/(T(z))*Tx(z)/LLength/void+1/(T(z))*dTdt(z)-Rho_p*(R*T(z)/P(z))*(1-void)/void*(dq1dt(z)+dq2dt(z)));

dPdt(z)= P(z)*(-ux(z)/LLength -u(z)/P(z)*Px(z)/LLength/void+u(z)/(T(z))*Tx(z)/LLength/void+1/(T(z))*dTdt(z)-Rho_p*(R*T(z)/P(z))*(1-void)/void*(dq1dt(z)+dq2dt(z)));

end

for z =1 : N

T_out(z)=T_ads;

end

for z = 2:N-1

dTdt(z)=-(Rho_gas(z)*Cp_gas*u(z)*Tx(z)/LLength-Rho_b*(H_Ads1*dq1dt(z)+H_Ads2*dq2dt(z))+2*Ui/R_bed*(T(z)-T_Wall(z)))/(tot_void*Rho_gas(z)*Cp_gas+Rho_b*Cp_s);

end

for z = 1 : N

dT_Walldt(z)=(2*3.14159*R_bed*Ui*(T(z)-T_Wall(z))-2*3.14159*R_bed_outside*h_outside_air*(T_Wall(z)-T_out(z)))/Rho_wall*Cp_wall*A_wall;

end

disp(t);

dydt = [dTdt;dPdt;dT_Walldt;dq1dt;dq2dt;dy1dt;dy2dt];

end

[Torsten]

One problem is that you included the boundary conditions after you computed the gradients. This should be done in advance.

And if you use SI units: are you sure that velocities in the order of 1 m/s are normal for your process ?

And you should remove the outer loop when you compute u:

%/*Ergun Equation*/
u(1)=u_feed;
for z = 2:N-1
    AA=1.75*Rho_gas(z)*(1-void)./(2*1.57e-3*void^3);
    BB=(150*mu*(1-void)^2)/(4*2.46e-6*void^3);
    CC= P1x(z)/LLength;
    u(z)=(-BB + sqrt(BB^2-4*AA*CC))/(2*AA);
    ux(z)=(u(z)-u(z-1))/(dz);
end
u(N)=u(N-1);

instead of

%/*Ergun Equation*/
for z=1:N
    u(1)=u_feed;
for z = 2:N-1
    AA(z)=1.75*Rho_gas(z)*(1-void)./(2*1.57e-3*void^3);
    BB=(150*mu*(1-void)^2)/(4*2.46e-6*void^3);
    CC(z)= P1x(z)/LLength;
    u(z)=(-BB + sqrt(BB^2-4.*AA(z).*CC(z)))./(2.*AA(z));
    ux(z)=(u(z)-u(z-1))/(dz);
end
    u(N)=u(N-1);
end

And in the former code, you didn't divide by LLength (which doesn't matter in this case because it's 1, but ...)

[William]

Hi,

Thank you very much for your comment—I appreciate your input.

I gave it a try, but unfortunately, there's no change in the results. When I place the boundary conditions before computing the gradient, MATLAB throws an error. It recognizes u(1) = u_feed, but not u(N) = u(N-1), so I have to apply the boundary condition exactly at that point.

I also removed the outer loop as you suggested, but the results remain unchanged.

Is there anything else you think might be causing the issue?

Just to clarify, the unit of u is m/s, while the gPROMS code doesn’t specify units—so I’ve assumed they are consistent.

Thanks again for your continued support.

Best regards,

William

[Torsten]

I used

%/*Boundary Condition*/%
dTdt(1)=(T_ads - T_des)/PressTime;
T(N)=T(N-1);
P(N)=Pfeed;
P(1)=P(2);
y1(1)=y_feed1;
y1(N)=y1(N-1);
y2(1)=y_feed2;
y2(N)=y2(N-1);

for z = 2:N %*BFDM1 Method*%
    Px(z)=(P(z)-P(z-1))/(dz);
    Tx(z)=(T(z)-T(z-1))/(dz);
    y1x(z)=(y1(z)-y1(z-1))/(dz);
    y2x(z)=(y2(z)-y2(z-1))/(dz);    
end

instead of

for z = 2:(N) %*BFDM1 Method*%
    Px(z)=(P(z)-P(z-1))/(dz);
    Tx(z)=(T(z)-T(z-1))/(dz);
    y1x(z)=(y1(z)-y1(z-1))/(dz);
    y2x(z)=(y2(z)-y2(z-1))/(dz);
    
    
    end
    %Px(1)=(-P(1)+P(2))/(dz);
    %Tx(1)=(-T(1)+T(2))/(dz);
    %y1x(1)=(-y1(1)+y1(2))/(dz);
    %y2x(1)=(-y2(1)+y2(2))/(dz);
    
    
   
%/*Boundary Condition*/%
dTdt(1)=(T_ads - T_des)/PressTime;
T(N)=T(N-1);
P(N)=Pfeed;
P(1)=P(2);
y1(1)=y_feed1;
y1(N)=y1(N-1);
y2(1)=y_feed2;
y2(N)=y2(N-1);

and got the gProms result (at least for the graph you included).

I'm really surprised about your setting P(1) = P(2) because without a pressure gradient, there is no flow.

[Torsten]

You get values from the integrator in the first and last point when you write the values in local variables:

    T=y(1:Nz);
    P=y(Nz+1:2*Nz);
    T_Wall=y(2*Nz+1:3*Nz);
    q1=y(3*Nz+1:4*Nz);
    q2=y(1+4*Nz:5*Nz);
    yy1=y(1+5*Nz:6*Nz);
    yy2=y(1+6*Nz:7*Nz);    

These values have to be overwritten by the boundary conditions before computing the gradients because they have no physical meaning.

I don't know how gProms handles this conflict between the values stored in the solution vector y for the unknowns and the values you explicitly define in the boundary points. MATLAB at least is a sequential programming language: things that should be used in line x have to be supplied somewhere before line x.