How to numerically solve a system of coupled partial differential and algebraic equations?

I have a system of coupled partial differential and algebraic equations.

Two 1-D parabolic pdes coupled (function of x and time) with two algebraic equations. What would be the way to solve this sytem?

[Torsten]

Without knowing your equations together with initial and boundary conditions, nothing useful can be said.

[Arpita]

Hi Torsten,

I really appreciate you getting back to me.

There is a possibility of singularity in the matrix. My equations include:

Solve for C1, C2, C3, φ.

χ1 -χ7 are known constants.

Initial Conditions:

Boundary conditions

[Torsten]

This will be hard to solve, and there is no MATLAB code that could help you.

The usual procedure is to discretize the spatial derivatives in equations (1) and (2) and solve the resulting system of differential-algebraic equations using ODE15S.

But you'll have a lot of trouble with it, that's for sure.

[Arpita]

Would that something similar to what's done in this example:

https://www.mathworks.com/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines

[Torsten]

Yes, the method-of-lines means to use difference quotients as approximations for the spatial derivatives (d/dx, d^2/dx^2) and transform the partial differential-algebraic system into an ordinary differential-algebraic system.

This can then be solved using a MATLAB solver like ODE15S.

[Arpita]

I convert my pdes into algebraic equations (finite difference method) as also done in the example. Do I just add the algebraic equations into the for loop?

[Torsten]

Yes, but you must tell the solver that these are algebraic and not differential equations.

This is done by zero rows within the mass matrix M which is part of the differential-algebraic system in the form of

M*y' = f(t,y)

[Arpita]

Okay, got it. Thank you!

Let me try that and get back to you.

[Arpita]

In solving these coupled PDEs and algebraic equations I am facing the following error:

"Warning: Failure at t=1.048722e-05. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.725809e-20) at time t."

I tried changing the solver to ode23t to skirt around stiffness issues, if any, but that doesn't solve the problem.

Setting 'Mass Singular" to 'yes" does not help either.

Could someone help with this problem?

[Torsten]

This means that ODE15S doesn't like your DAE system.

The reasons can be manifold: programming errors, difficulty of the system to be integrated, ...

So no advice can be given here.

[Arpita]

I don't think I have any programming errors, but could you possibly help me?

If the issue is the difficulty of the system to be integrated, should I move to a different software?

clear
clc
%% Parameters %%
% Experimental conditions %
p.Q1=37.33/60;% 4.4664/60; % flowrate to the stack (mL/s)
p.I=-1258/1120*10^(-3); % applied current density (A/cm2)
p.T1=500;% round(41/Q1); % charging time(s)
p.C0=0.0342; % initial concentration (mmol/ml)

% Electrode and cell physical properties %
p.Mass=10; %electrode weight (g)
p.N=1; % number of electrode pairs
p.psp=0.708; % spacer porosity
p.pma=0.4; % macroprosity
p.pmi=0.28;  % microporosity
p.Lelec=0.028; % each electrode thickness (cm)
p.Lsp=0.01;  % spacer thickness (cm)
p.a=1120; % each electrode area (cm^2)
p.hrt1=p.Lsp*p.a*p.N/p.Q1; % hydraulic retention time (s) in the spacer
p.tlag=round(1*p.hrt1);
p.alpha=50;
p.Rext=200;% resistance (om*cm2)

% Electrochemical properties %
p.De=1.68*10^(-5);  % effective diffusion coefficient (cm^2/s)
p.R=6.1;% specific electrode resistance (om*mmol/cm), from zhao, WR, 'optimation' 2013
p.u=0; % chemical attraction term (kT) 
p.Vt=0.0256; % thermal voltage (V)
p.Cst1=72; % stern capacitance (F/mL) value from JE.Dykstra W 2016.
p.F=96.5; % farady constant (C/mmol)
p.i=p.I/p.F; % convert (A/cm2) to (mmol/s/cm2)
p.Vapplied = (1.2/2/p.Vt);
p.V1 = p.Vapplied-(p.i*p.Lsp/(4*p.C0*p.psp*p.De));

% using method of lines %
tf      = 200;                               % End time
p.nz    = 100;                                  % Number of nodes
x       = linspace(0,p.Lelec,p.nz);                 % Space domain
dx      = diff(x);
p.dx = dx(1);                % Space increment

%% Initial conditions %%
Ceff_0 = p.C0*(p.pma+p.pmi); % Ceff(x,0)= C0
Cch_0 = 0; % Cch(x,0) = 0
Cma_0 = p.C0; % Cma(x,0) = C0
phima_0 = p.V1; %phima(x,0) = V1
Csp = p.C0.*ones(p.nz,1);
Ceff = ones(p.nz,1);
Ceff=Ceff_0.*Ceff;
Cch = zeros(p.nz,1);
Cma = ones(p.nz,1);
Cma=Cma_0.*Cma;
phima = ones(p.nz,1);
phima=phima_0.*phima;
%y0      = [Csp;Ceff;Cch;Cma;phima];
y0      = [Ceff;Cch;Cma;phima];
tspan   = [0 tf];

%% Boundary conditions %%
dCmadx = 0; % dCmadx(Lelec,t) = 0
dphimadx =0; % dphimadx(Lelec,t) =0

%% Mass matrix (for DAE system) %%
M = eye(4*p.nz);
for i=(2*p.nz+1):(4*p.nz)
    M(i,i)=0;
end
options = odeset('Mass',M,'RelTol',1e-3,'AbsTol',1e-10);%'MassSingular','yes');
[t,y] = ode23t(@(t,y)porada1D(t,y,p),tspan,y0,options);

%% Plot %%
Ceff_matrix=y(:,1:p.nz);
Cch_matrix=y(:,p.nz+1:2*p.nz);
Cma_matrix=y(:,2*p.nz+1:3*p.nz);
phima_matrix=y(:,3*p.nz+1:4*p.nz);
plot(t,Ceff_matrix(:,:))
grid
xlabel('time')
ylabel('Ceff')
%title('Dimensionless concentration at ''z = L''')
%figure()
%plot(t, phima_matrix(:,:))


%% function %%
function out = porada1D(~,y,p)
% Variables
Ceff = y(1:p.nz,1);
Cch = y(p.nz+1:(2*p.nz),1);
Cma = y((2*p.nz)+1:(3*p.nz),1);
phima = y((3*p.nz)+1:(4*p.nz),1);
%phima = y((4*p.nz)+1:(5*p.nz),1);

% Space derivatives
dCmadx = zeros(p.nz,1);
d2Cmadx2 = zeros(p.nz,1);
dphimadx = zeros(p.nz,1);
d2phimadx2 = zeros(p.nz,1);
dCmadx(2:end-1) = (Cma(3:end)-Cma(1:end-2))./(2.*p.dx);
d2Cmadx2(2:end-1) = (Cma(3:end)-(2*Cma(2:end-1))+Cma(1:end-2))./(p.dx^2);
dphimadx(2:end-1) = (phima(3:end)-phima(1:end-2))./(2.*p.dx);
d2phimadx2(2:end-1) = (phima(3:end)-(2*phima(2:end-1))+Cma(1:end-2))./(p.dx^2);

% Governing equations
%dCspdt = (2.*p.pma.*p.De/p.psp/p.Lsp).*dCmadx + ((1./p.psp./p.hrt1).*(p.C0-Csp));
dCeffdt = p.pma.*p.De.*d2Cmadx2;
dCchdt = 2.*(p.pma./p.pmi).*p.De.*(dCmadx.*dphimadx+Cma.*d2phimadx2);
phit = phima+(asinh(Cch./(-2.*Cma)))-(p.F.*Cch./(p.Vt.*(p.Cst1+p.alpha.*(Cch.*Cch))))-p.V1;
Cefft = Ceff-p.pma.*Cma-(p.pmi.*Cma.*cosh(Cch./(-2.*Cma)));

% Boundary conditions
Cma(1,1)=p.C0;
dCmadx(end)=0;
dphimadx(end)=0;

%out = [dCspdt;dCeffdt;dCchdt;phit;Cefft];
out = [dCeffdt;dCchdt;phit;Cefft];
end

[Torsten]

% Boundary conditions
Cma(1,1)=p.C0;
dCmadx(end)=0;
dphimadx(end)=0;

The boundary conditions you set at the end of your function "porada1D" are nowhere used in the calculation of your output vector

out = [dCeffdt;dCchdt;phit;Cefft];

This cannot be correct.

[Arpita]

The version I sent you had the boundary conditions after the governing equations. Sorry about that.

Even with boundary conditions before the governing equations, I see the same problem.

With all three boundary conditions, I get the error that the DAE is greater than 1.

When I only use the dCmadx = 0 and dphimadx = 0 as the boundary conditions, I get the error: " Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-19) at time t."

[Torsten]

You should start with a simpler problem from which you know that potentially arising problems with the integrator stem from your programming, not from the difficulty or even unsolvability of the problem itself.

Answer 1

[Sarthak]

Hi,

One way to solve a system of coupled partial differential equations (PDEs) and algebraic equations is to use a numerical method such as finite difference or finite element method. Here is an outline of the steps involved:

1. Discretize the system of PDEs using a numerical method such as finite difference or finite element method. This will transform the PDEs into a system of algebraic equations.
2. Combine the discretized PDEs with the algebraic equations to form a system of nonlinear algebraic equations.
3. Use a numerical solver such as the Newton-Raphson method or a quasi-Newton method to solve the system of nonlinear algebraic equations.
4. Repeat the process for each time step to obtain a time-dependent solution.