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.
Update: Extended Answer
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_] :=
TraditionalForm[
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",
0, ρ*Cp*Velocity.Inactive[Grad][T, X]];
Q = If[Source === "NoSource", 0, Source];
If[ FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
If[ VectorQ[a], a = DiagonalMatrix[a]];
(* Note the - sign in the operator *)
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 \
about Being Convection Dominated *)
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=``"][
ToString@PaddedForm[t, {4, 2}]], 12],
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]
Original Answer
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 \
about Being Convection Dominated *)
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]
