# Method of lines - Dirichlet and mixed BC

I have a dissolution problem to solve with two equations (everything is in dimensionless form - concentration, time and distance - EDIT: that came from the second Fick's law, where the distance was normalized by L, time by L^2/D_H and S by S_h^0) :

that represents diffusion equation for proton transport across the fluid boundary layer in dimensionless form;

and

that represents the transport equation for the porous layer.

There are three layers (rock, porous media and liquid media), three boundary conditions, and one initial condition:

My main challenge is to solve all the equations simultaneously to obtain the rate:

.

EDIT: where J_H is the first Fick law, Bi is the Biot number

In summary, that is the representation of my problem

and this is what I'm trying to obtain: Dissolution rate as a function of layer thickness for different values of ε; and Dissolution rate as a function of layer thickness for different Biot numbers (Bim).

and

EDIT: Im tryng to reproduce those results: https://www.sciencedirect.com/science/article/abs/pii/S0009254122001620?via%3Dihub

Here is my code, but unfortunately it is not working.

EDIT: I change my code but is still not working

(* Problem parameters *)
epsilon = 0.1;
Bi = 2.0;
lValues = Range[0, 10, 1];  (* Different values of l for the R plot *)

(* Discretization of space and time *)
Nx = 100;  (* Number of points in x *)
Nt = 200;  (* Number of points in time *)
dx = 1.0/Nx;
dt = 1.0/Nt;

(* Matrix initialization *)
SH = ConstantArray[0, {Nx, Nt + 1}];

(* Initial condition *)
SH[[All, 1]] = ConstantArray[0, Nx];

(* Main loop - Finite difference method *)
For[n = 1, n <= Nt, n++,
For[i = 2, i <= Nx - 1, i++,
SH[[i, n + 1]] = SH[[i, n]] +
dt*((SH[[i + 1, n]] - 2 SH[[i, n]] + SH[[i - 1, n]])/(dx^2) -
epsilon*(SH[[i + 1, n]] - SH[[i - 1, n]])/(2*dx)
)
];

(* Central boundary condition *)
SH[[1, n + 1]] = 1;

(* Boundary condition at the end *)
SH[[-1, n + 1]] = Bi*SH[[-1, n]]*dx + SH[[-2, n + 1]];
];

(* Plotting the SH graph over x *)
xValues = Range[0, 1, dx];
ListLinePlot[SH[[All, -1]], DataRange -> {0, 1}, PlotLabel -> "Distribution of SH along x",
AxesLabel -> {"x", "SH"}, ImageSize -> Medium]

(* Calculating R for different l values *)
RValues = Table[Bi*SH[[Round[l*Nx], -1]], {l, lValues}];

(* Plotting the R graph over l *)
ListPlot[Transpose[{lValues, RValues}], PlotStyle -> PointSize[Medium],
PlotLabel -> "R values for different l values",
AxesLabel -> {"l", "R"}, ImageSize -> Medium]

• 1. Please double check those formulas, for example, some $S_H$ is mistakenly written as $S_h$ in your picture. 2. Why do you write a difference formula in i.sstatic.net/c0QOC.png and claim it's a boundary condition? What does it mean? 3. What's $J_H$? 4. Why not directly use NDSolve? Dec 3, 2023 at 2:16
• 1) Please also present your equations as Mathematica code. 2) Have you tried NDSolve rather than writing your own solver? Dec 3, 2023 at 2:40
• I can't use NDSolve.. I really need to write my own code Dec 3, 2023 at 15:52
• In what sense does your code not work? Errpr messages? Clearly wrong answers? Please be specific. Dec 4, 2023 at 4:12

We can solve this problem using FDM and NDSolve as well to compare solutions. First code with FDM implementation

R[L_, Bi_, \[CurlyEpsilon]_] :=
Module[{SpatialStep = 0.1, TemporalStep = 0.001, FinalTime = 1.0},
nTemporalSteps = Round[FinalTime/TemporalStep];
grid1 = Range[0, 1, SpatialStep]; nx1 = Length[grid1];
grid2 = Range[1, L + 1, SpatialStep]; nx2 = Length[grid2];

(*Initialize the grid*)
SH0[x_] := Exp[-(x)^2]; (*Initial condition for SH*)
Sh0[x_] :=
Exp[-1]*(1 -
1/Pi Sin[2 \[Pi] (x - 1)/L]); (*Initial condition for Sh*)

(*Lists to store the results*)
solutionSH = Table[0, {nx1}, {nTemporalSteps}];
solutionSh = Table[0, {nx2}, {nTemporalSteps}];

(*Initial conditions*)
solutionSH[[All, 1]] = SH0[#] & /@ grid1;
solutionSh[[All, 1]] = Sh0[#] & /@ grid2;
M21 = NDSolveFiniteDifferenceDerivative[Derivative[2], grid1,
DifferenceOrder -> 2]@"DifferentiationMatrix";
M22 = NDSolveFiniteDifferenceDerivative[Derivative[2], grid2,
DifferenceOrder -> 2]@"DifferentiationMatrix";

(*Solve the equations*)
For[n = 1, n < nTemporalSteps, n++,
solutionSH[[2 ;; -2, n + 1]] =
solutionSH[[2 ;; -2, n]] +
TemporalStep (M21 . solutionSH[[All, n]])[[2 ;; -2]];
(*Dirichlet boundary condition at x=0*)
solutionSH[[1, n + 1]] = 1;

solutionSh[[All, n + 1]] =
solutionSh[[All, n]] +
TemporalStep \[CurlyEpsilon] M22 . solutionSh[[All, n]];
(*Boundary conditions at x=1*)
solutionSH[[-1, n + 1]] =
solutionSH[[-2,
n + 1]] + \[CurlyEpsilon] (solutionSh[[2, n + 1]] -
solutionSh[[1, n + 1]]);
solutionSh[[1, n + 1]] = solutionSH[[-1, n + 1]];
(*Mixed boundary condition at x=l+1*)
solutionSh[[-1,
n + 1]] = (solutionSh[[-2, n + 1]])/(1 +
Bi /\[CurlyEpsilon] SpatialStep );];
Bi solutionSh[[-1, -1]]];


Visualization $$R_{AB}(\epsilon)$$

RAB = Table[
Table[{L, R[L, 1, eps]}, {L, 1, 10, .1}], {eps, {.5, .1, .01}}];
p1 = ListLinePlot[RAB,
PlotLegends ->
Table[Row[{"\[CurlyEpsilon] =", eps}], {eps, {0.5, .1, .01}}],
PlotRange -> All, PlotLabel -> Row[{"Bi =", 1}], Frame -> True,
FrameLabel -> {"Layer thickness, L",
"Dissolution rate, \!$$\*SubscriptBox[\(R$$, $$AB$$]\)"},
PlotStyle -> {{Blue, Dashed}, {Red, Dashed}, {Green, Dashed}}]


Note, that for initial condition used above the time integration FinalTime = 1. is too small. It is why we see waves in the plot p1. Second code

d[x_, e_] := Piecewise[{{1, 0 <= x < 1}, {e, True}}];
ini[x_, L_] :=
Piecewise[{{Exp[-x],
0 <= x < 1}, {Exp[-1]*(1 - 1/Pi Sin[2 \[Pi] (x - 1)/L]), True}}];
R1[L_, Bi_, eps_] :=
Module[{},
sol = NDSolveValue[{D[u[t, x], t] - D[d[x, eps] D[u[t, x], x], x] ==
NeumannValue[-Bi u[t, x], x == L + 1],
DirichletCondition[u[t, x] == 1, x == 0], u[0, x] == ini[x, L]},
u, {t, 0, 10}, {x, 0, L + 1}]; sol];


Visualization $$R_{AB}(\epsilon)$$

p2 = ListLinePlot[
Table[Table[{L, R1[L, 1, eps][1, L + 1]}, {L, 1,
10, .1}], {eps, {0.5, .1, .01}}],
PlotLegends ->
Table[Row[{"\[CurlyEpsilon] =", eps}], {eps, {0.5, .1, .01}}],
PlotRange -> All]


Now we can compare two solutions

Show[p1,p2]


We can see a slight difference due to errors of $$h^2$$ (h=SpatialStep) in the FDM solution. Finally we can reproduce Figure 2,a from the paper as follows

list10 =
Table[Table[{L, R1[L, 1, eps][10, L + 1]}, {L, .01,
10, .1}], {eps, {0.5, .1, .01}}];

ListLinePlot[list10,
PlotLegends ->
Table[Row[{"\[CurlyEpsilon] =", eps}], {eps, {0.5, .1, .01}}],
PlotRange -> All, Frame -> True,
FrameLabel -> {"Layer thickness, L",
"Dissolution rate, \!$$\*SubscriptBox[\(R$$, $$AB$$]\)"}]