# NDSolve for system of PDEs -Error: fewer dependent variables!

I am trying to solve the equations below to find the breakthrough curve in an adsorption process:

∂C/∂t+u* ∂C/∂z+ρ/ϵ * ∂q/∂t=  0
∂q/∂t= -Ks * (q-(qm* b*C)/(1+ b*C))

the IC and BCs are:
q(t=0,z)=0
C(t,z=0)=C0
C(t=0,z) = 0


The parameters are given as:

eps = 0.373
ueps = 6.66
C0 = 892
Ks = 0.0872
rho = 390.04
qm = 7.0136
b = 0.05313


I wrote the equation in Mathematica as below:

s=NDSolve[{∂_t C[t,z]==-ueps*∂_z C[t,z]+(rho⁄eps)*Ks*(q[t,z]-(qm*b*(C[t,z])⁄((1+b*C[t,z])))),∂_t q[t,z]==-Ks*(q[t,z]-(qm*b*(C[t,z])⁄((1+b*C[t,z])))),DirichletCondition[{C[t,0]==C0},True],C[0,z]==0,q[0,z]==0},{C,q},{t,0,60},{z,0,22}]


After evaluating the notebook, I receive this error: NDSolve: There are fewer dependent variables, {q[t,z]}, than equations, so the system is overdetermined.

I have been trying to fix this problem by rearranging the equations, but no success yet.

Thanks

• Maybe my post on Chemical adsorption in fixed beds on the Wolfram Community website will give you some ideas. Commented Jul 27, 2020 at 21:00
• Keep in mind, C is a reserved letter. Commented Jul 27, 2020 at 21:17

# Update: Mass Transport Tutorial

In my answer, I adapted the Heat Transfer Tutorial to handle Mass Transport. There is a new and fairly extensive Mass Transport Tutorial that you can use directly. In my answer 227821, I extended the Mass Transport functions in the tutorial to include axisymmetric and optional porosity.

To use the FEM method it is good to cast your equations into coefficient form as shown FEM Tutorial.

$$\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

The benefits of coefficient form include:

• Standardization so others can quickly interpret your equations
• Avoidance of sign errors
• Nice one-to-one mapping with other FEM codes such as COMSOL

For advective-diffusive problems with and without source terms, the Heat transfer tutorial is a good place to start to learn how to set up these types of problems. It shows how to properly setup the PDE operators in coefficient form. The TimeHeatTransferModel function in the Heat Transfer Tutorial can help guide you on how to setup the problem.

## Create Fluid and Solid PDE Operators Using Heat Transfer Tutorial

The following shows how one can setup your coupled system of interphase transport of the fluid and solid using the TimeHeatTransferModel :

(* From Vitaliy Kaurov for nice display of operators *)
pdConv[f_] :=
f /. Derivative[inds__][g_][vars__] :>
Apply[Defer[D[g[vars], ##]] &,
Transpose[{{vars}, {inds}}] /. {{var_, 0} :>
Sequence[], {var_, 1} :> {var}}]]
(*PDEModels/tutorial/HeatTransfer/HeatTransferVerificationTests#\
463435833*)
ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, ρ_, Cp_, Velocity_, Source_] :=
Module[{V, Q, a = k},
V = If[Velocity === "NoFlow",
Q = If[Source === "NoSource", 0, Source];
If[ FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
If[ VectorQ[a], a = DiagonalMatrix[a]];
a = PiecewiseExpand[Piecewise[{{-a, True}}]];
Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]
TimeHeatTransferModel[T_, TimeVar_, X_List, k_, ρ_, Cp_,
Velocity_, Source_] := ρ*Cp*D[T, {TimeVar, 1}] +
HeatTransferModel[T, X, k, ρ, Cp, Velocity, Source]
(* Create Parametric PDE Operators for Fluid and Solid *)
(* Include Small Diffusive Term to Fluid Otherwise MMA Might Complain \
parmfop =
TimeHeatTransferModel[c[t, x], t, {x}, {d}, 1,
1, {u}, -(ρ/ϵ) Q];
parmsop =
TimeHeatTransferModel[q[t, x], t, {x}, {0}, 1, 1, "NoFlow", Q];
parmfop // pdConv
parmsop // pdConv


Slight rearrangement shows the operators follow coefficient form:

$$\begin{matrix} m\frac{{{\partial ^2}}}{{\partial {t^2}}}u & +\ \ d\frac{\partial }{{\partial t}}u & + \ \ \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) u & + \ \ \beta \cdot\nabla u & + \ \ au & -\ \ f & =0 & \ \ (1) \\ & & & & & & \\ & +\ \ \frac{{\partial c(t,x)}}{{\partial t}} & + \ \ {\nabla _{\left\{ x \right\}}}\cdot\left( {\left( {\begin{array}{*{20}{c}} { - d} \end{array}} \right).{\nabla _{\left\{ x \right\}}}c\left( {t,x} \right)} \right) & + \ \ \left\{ u \right\}.{\nabla _{\left\{ x \right\}}}c\left( {t,x} \right) & & +\ \ \frac{{\rho Q}}{\varepsilon } & =0 & \ \ (2) \\ & & & & & & \\ & +\ \ \frac{{\partial q(t,x)}}{{\partial t}} & & & & -\ \ Q & =0 & \ \ (3) \end{matrix}$$

## New Workflow

The following workflow uses the PDE operators defined above and produces an annotated plots that can be scrolled frame-by-frame with manipulate.

(* Define Parameters *)
eps = 0.373;
ueps = 6.66;
c0 = 892;
Ks = 0.0872;
rho = 390.04;
qm = 7.0136;
b = 0.05313;
tend = 60;
xend = 22;
Qsource = -Ks*(q[t, x] - (qm*b*c[t, x])/(1 + b*c[t, x]));
timeinc = 0.25; (* Output times *)
(* Critical Breakthrough Fraction *)
btfrac = 0.05;
(* Use CombinePlots For Secondary Axis *)
cp = ResourceFunction["CombinePlots", "Function"];
(* Setup PDE System *)
pdef = (parmfop == 0) /. {ϵ -> eps, u -> ueps, ρ -> rho,
d -> 1, Q -> Qsource};
pdes = (parmsop == 0) /. {Q -> Qsource};
dc = DirichletCondition[c[t, x] == c0 (1 - Exp[-40 t]), x == 0];
icf = c[0, x] == 0;
ics = q[0, x] == 0;
(* Solve PDE *)
{cifn, qifn} =
NDSolveValue[{pdef, pdes, dc, icf, ics}, {c, q}, {t, 0, tend}, {x,
0, xend},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.5}}}];
(* Plot Solutions *)
(* Define Plot Function Using Combine Plots *)
bfn = Function[{t},
If[cifn[t, xend] >= btfrac c0, xend,
x /. FindRoot[cifn[t + 0.01, x] - btfrac c0, {x, 0, xend},
Method -> "Brent"]]];
bttemp = With[{xbt = #},
If[xbt == xend, "Sat",
StringTemplate["x="][ToString@PaddedForm[xbt, {4, 2}]]]] &;
pltfn = Function[{t}, cp[
With[{xbt = bfn[t]},
Plot[Callout[cifn[t, x], bttemp[xbt], xbt, RoundingRadius -> 5,
Appearance -> "SlantedLabel"], {x, 0, xend},
PlotRange -> {{-2, xend + 1}, {-10, 1.1 c0}}, Frame -> True,
FrameLabel -> {"x", "c"},
PlotLabel ->
Style[StringTemplate[
"Fluid and Solid Concentrations @ time="][
Epilog -> {Green, Thick, Dashed,
Line[{{xbt, -10}, {xbt, 1.1 c0}}]}]],
Plot[qifn[t, x], {x, 0, xend}, PlotPoints -> 100,
PlotRange -> {0, 7},
Frame -> True, FrameStyle -> Red, PlotStyle -> Red,
FrameLabel -> {"x", "q"}],
"AxesSides" -> "TwoY"
], Listable];
(* Create Image List and Manipulate *)
imgs = Rasterize[#, ImageSize -> 340] & /@
pltfn[Range[0, tend, timeinc]];
Manipulate[imgs[[i]], {i, 1, tend/timeinc, 1},
ControlPlacement -> Top]


I would cast the problem using the following workflow:

(* Define Parameters *)
eps = 0.373;
ueps = 6.66;
c0 = 892;
Ks = 0.0872;
rho = 390.04;
qm = 7.0136;
b = 0.05313;
tend = 60;
xend = 22;
Qsource = -Ks*(q[t, x] - (qm*b*c[t, x])/(1 + b*c[t, x]));
(* Use CombinePlots For Secondary Axis *)
cp = ResourceFunction["CombinePlots"];
(* Define Parametric PDE operators *)
(* Include Small Diffusive Term to Fluid Otherwise MMA Might Complain \
parmfop =
D[c[t, x], t] +  D[-d D[c[t, x], x], x] + u  D[c[t, x], x] + (
Q ρ)/ϵ;
parmsop = D[q[t, x], t] - Q;
(* Setup PDE System *)
pdef = (parmfop == 0) /. {ϵ -> eps, u -> ueps, ρ -> rho,
d -> 1, Q -> Qsource};
pdes = (parmsop == 0) /. {Q -> Qsource};
dc = DirichletCondition[c[t, x] == c0 (1 - Exp[-40 t]), x == 0];
icf = c[0, x] == 0;
ics = q[0, x] == 0;
(* Solve PDE *)
{cifn, qifn} =
NDSolveValue[{pdef, pdes, dc, icf, ics}, {c, q}, {t, 0, tend}, {x,
0, xend},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.5}}}];
(* Plot Solutions *)
(* Define Plot Function Using Combine Plots *)
pltfn = Function[{t}, cp[
Plot[cifn[t, x], {x, 0, xend}, PlotRange -> {0, 1.1 c0},
Frame -> True, FrameLabel -> {"x", "c"}],
Plot[qifn[t, x], {x, 0, xend}, PlotPoints -> 100,
PlotRange -> {0, 7},
Frame -> True, FrameStyle -> Red, PlotStyle -> Red,
FrameLabel -> {"x", "q"}],
"AxesSides" -> "TwoY"
], Listable];
(* Create Image List and Animate *)
imgs = Rasterize /@ pltfn[Range[0, tend, 0.25]];
ListAnimate[imgs, ControlPlacement -> Top]