2
$\begingroup$

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) :

enter image description here

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

and enter image description here

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: enter image description here

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

enter image description here.

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

In summary, that is the representation of my problem

enter image description here

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).

enter image description here and enter image description here

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]
$\endgroup$
4
  • $\begingroup$ 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.stack.imgur.com/c0QOC.png and claim it's a boundary condition? What does it mean? 3. What's $J_H$? 4. Why not directly use NDSolve? $\endgroup$
    – xzczd
    Dec 3, 2023 at 2:16
  • 1
    $\begingroup$ 1) Please also present your equations as Mathematica code. 2) Have you tried NDSolve rather than writing your own solver? $\endgroup$
    – MarcoB
    Dec 3, 2023 at 2:40
  • $\begingroup$ I can't use NDSolve.. I really need to write my own code $\endgroup$ Dec 3, 2023 at 15:52
  • 1
    $\begingroup$ In what sense does your code not work? Errpr messages? Clearly wrong answers? Please be specific. $\endgroup$
    – bbgodfrey
    Dec 4, 2023 at 4:12

1 Answer 1

3
$\begingroup$

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 = NDSolve`FiniteDifferenceDerivative[Derivative[2], grid1, 
      DifferenceOrder -> 2]@"DifferentiationMatrix"; 
   M22 = NDSolve`FiniteDifferenceDerivative[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}}]

Figure 1 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]  

Figure 2 Now we can compare two solutions

Show[p1,p2]

Figure 2 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\)]\)"}]

Figure 4

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.