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

Question

Arpita 2023-1-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1900755-how-to-numerically-solve-a-system-of-coupled-partial-differential-and-algebraic-equations

回答： Sarthak 2023-3-9

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?

14 个评论
显示 12更早的评论隐藏 12更早的评论

Arpita 2023-2-2

编辑：Torsten 2023-2-2

在 MATLAB Online 中打开

Porada1D_2023_01_31.m

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

Arpita 2023-2-3

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 2023-2-3

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 2023-3-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1900755-how-to-numerically-solve-a-system-of-coupled-partial-differential-and-algebraic-equations#answer_1189025

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:

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.
Combine the discretized PDEs with the algebraic equations to form a system of nonlinear algebraic equations.
Use a numerical solver such as the Newton-Raphson method or a quasi-Newton method to solve the system of nonlinear algebraic equations.
Repeat the process for each time step to obtain a time-dependent solution.