Building on the system of equations in this post, I attempted to solve an additional convection-diffusion equation describing the concentration of solute in the lens, which affects its spreading.

enter image description here

However, I noticed that the LHS is multiplied by h2-h1. As h2-h1 could become very small or even negative, it results in the concentration field destabilising and producing sharp peaks (ie a sudden gain of solute).

enter image description here

I have attempted to resolve the issue using unitStepExpand@If[(h2 - h1) < hinf, hinf, (h2 - h1)] as in this post (Error in Attempting Moving Boundary Fluid System). However, NDSolveValue fails to initialise after a long time and the kernel crashes.

The code can be found below, and the definition of pdetoode can be found here. I am currently using a graded mesh that is finer in the region with the lens. I have also tried a uniform grid, as can be created below:

domain = {-20, 20}; 
points = 120; 
grid = Array[# &, points, domain]; 

This initialises but uses extremely small timesteps of 10^-7 s.

Thus, how could I ensure that h2-h1 remains positive and non-zero? Thank you!

With[{h1 = h1[r, t], h2 = h2[r, t], c = c[r, t]},
 p1[c_] := -(\[Gamma]ow[c]/r)*D[r*D[h1, r], r] - \[Gamma]w[c]/r*
    D[r*D[h2, r], r];
 p2[c_] := -\[Gamma]w[c]/r*D[r*D[h2, r], r] - \[CapitalPi]dep2[
    h2 - h1, c];]

unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;
With[{h1 = h1[r, t], h2 = h2[r, t], c = c[r, t], P1 = P1[r, t], 
  P2 = P2[r, t], q1 = q1[r, t], q2 = q2[r, t]},
 Acoeff = 
  Simplify[-h1*D[P1, r] - (h2 - h1)*D[P2, r] + D[\[Gamma]w[c], r] + 
    D[\[Gamma]ow[c], r]];
 Ccoeff = D[\[Gamma]w[c], r] - h2*D[P2, r];
 Dcoeff = 
  Simplify[Mrat*(-(1/2)*h1^2*D[P1, r] - h1*(h2 - h1)*D[P2, r] + 
       h1*(D[\[Gamma]w[c], r] + D[\[Gamma]ow[c], r])) - 
    1/2*h1^2*D[P2, r] - h1*D[\[Gamma]w[c], r] + h1*h2*D[P2, r]];
 eqp1 = P1 == p1[c];
 eqp2 = P2 == p2[c];
 eqq1 = q1 == Simplify[1/\[Mu]o*(1/6*h1^3*D[P1, r] + 1/2*h1^2*Acoeff)];
 eqq2 = q2 == 
    1/\[Mu]w*(1/6*D[P2, r]*(h2^3 - h1^3) + Ccoeff/2*(h2^2 - h1^2) + 
       Dcoeff*(h2 - h1))];
 eqh1 = D[h1, t] == -1/r*D[r*q1, r];
 eqh2 = D[h2, t] == -1/r*D[r*q1, r] - 1/r*D[r*q2, r] - J[c];
 eqc =(*(h2-h1)*)
  unitStepExpand@If[(h2 - h1) < hinf, hinf, (h2 - h1)]*
    D[c, t] == (diff/r*D[r*(h2 - h1)*D[c, r], r] - q2/r*D[r*c, r] - 
 ic = {h1 == Hbath, h2 == hinit, c == cinit} /. t -> 0;

hinf = 0.00001; 
ccrit = 0.32152;(*from simulation data used here, where S=0*)

\[CapitalPi]dep2[h_, c_] := (2*Sfunc[c])/
  hinf*(hinf/h)^3*(1 - Sfunc[0]/Sfunc[c]*(hinf/h)^3);

(*interfacial tensions*)
\[Gamma]ow[c_] := 0.0312 - 0.0177*c + 7.97*10^-3*c^2;
\[Gamma]o = 0.0309;
\[Gamma]w[c_] := 4.94*10^-3 - 8*10^-4*c - 1.12*10^-3*c^2;
Sfunc[c_] := \[Gamma]o - \[Gamma]ow[c] - \[Gamma]w[c];
Plot[Sfunc[c], {c, 0, 1}]

\[Mu]o = 4.88*10^-5;(*https://downloads.hindawi.com/archive/2014/\
\[Mu]w = 2.75*10^-6;(*https://www.rheosense.com/applications/\
Mrat = \[Mu]w/\[Mu]o;

diff = 0.4*10^-3;(*diffusion coefficient*)
J[c_] := 7.3*10^-4*c;(*evap rate*)

(*initial conditions*)
hinit = hinf + Hbath + 0.25 E^(-r^2);
cinit = 0.45*0.5 (1 + Tanh[1.5*(1.2 - Abs[r])]);

(*Generate graded grid*)

Hbath = 3;
ldrop = 10; ndrop = 90;(*number of grid points in the region*)
loil = 10; noil = 40;
dropgrid = Subdivide[ldrop, ndrop];(*non-graded mesh for drop*)
size = dropgrid[[2]]
oilgridmid = Table[size + (1/8*i)^2, {i, 0, loil, loil/noil}];
oilgrid = 
 Table[ldrop + size*i, {i, 0, Length[oilgridmid] - 1}] + oilgridmid
combined = Union[dropgrid, oilgrid];
reflected = DeleteCases[Union[combined, -combined], 0];
toplot = Transpose[{reflected, Table[1, {i, 1, Length[reflected]}]}];
ListPlot[toplot(**), PlotStyle -> PointSize[Small]]

(*Discretisation and Solving*)

grid = reflected;
difforder = 2;
domain = {Min[grid], Max[grid]}

tfunc = pdetoode[{h1, h2, c, P1, P2, q1, q2}[r, t], t, grid, 
   difforder, True];
odemid = Map[tfunc, {eqp1, eqp2, eqq1, eqq2}, {2}];
ode = Block[{P1, P2, q1, q2}, Set @@@ odemid;
   tfunc@{eqh1, eqh2, eqc}];
odeic = tfunc@ic;
var = Outer[#[#2] &, {h1, h2, c}, grid];

timer = 0.11;
time = 0;

Monitor[sol = 
   NDSolveValue[{ode, odeic}, var, {t, 0, timer}, 
    Method -> {"EquationSimplification" -> "Residual"}, 
    EvaluationMonitor :> (time = t)], time];

{h1sol, h2sol, csol} = rebuild[#, grid, 2] & /@ sol;
 Plot[{h1sol[r, t], h2sol[r, t], csol[r, t]}, {r, ##}, 
    PlotRange -> {{-5, 5}, {0, 3.3}}, 
    PerformanceGoal -> "Quality"] & @@ domain, {t, 0, timer}]


Using a much finer graded grid, with 500 points from {0,10} and 50 points from {10,30}, we observe that the droplet concentration forms a clear step function before eventually exploding. This verifies that the instability is likely caused by the sharp discontinuity in the concentration field.

Since the numerical instability is caused by discontinuities in the solution field itself, similar to the Burger's equation, are there any potential strategies implementable for we to continue using the finite difference scheme? Alternatively, would discretising our problem in a separate domain i.e. finite element scheme allow us to maintain a stable solution despite the discontinuity in the solution field?

enter image description here

  • $\begingroup$ I doubt if you're in the right direction. 1. According to my (limited) experience, non-uniform grid for FDM can easily lead to instability, that's the reason I myself don't often use it. 2. If you use the eqc without unitStepExpand and choose a coarse grid e.g. domain = {-20, 20}; points = 24; grid = Array[# &, points, domain];, the problem is solved, but problem that can only be solved with coarse grid is sometimes a sign of ill-posed problem. (Example: mathematica.stackexchange.com/a/69918/1871 ) Are you sure the problem itself is correct? Can you add a reference? $\endgroup$
    – xzczd
    Commented Dec 20, 2022 at 13:43
  • $\begingroup$ While there isn't an exactly identical model in literature, the research conducted by this paper sci-hub.wf/10.1021/la9019469 is quite close. The main difference is that while their spreading is surfactant driven and on solid surfaces, while our spreading is concentration driven and on another fluid. They have managed to solve the diffusion convection equation in the same form, though granted the convection of their fluid is much slower. $\endgroup$
    – FLP
    Commented Dec 20, 2022 at 14:03
  • $\begingroup$ Another paper which is quite similar to ours is this one users.auth.gr/gkarapetsas/Publications/…, where they spread a surfactant on a liquid surface. Similarly the original authors managed to solve the diffusion convection equation in the same form, although they had to use finite element as their discretisation scheme instead. $\endgroup$
    – FLP
    Commented Dec 20, 2022 at 14:04
  • $\begingroup$ Btw if our problem is really ill-posed, is there any other ways we could potentially solve it, for instance by discretisation via finite element? $\endgroup$
    – FLP
    Commented Dec 21, 2022 at 9:54
  • $\begingroup$ FEM is not magic. For problem defined in regular domain, I don't think FEM will show any advantage (as far as the methods have been properly implemented). For any FEM implementation, when element with regular shape is chosen, one can always find an equivlant FDM implementation, I believe. $\endgroup$
    – xzczd
    Commented Dec 21, 2022 at 10:07


Your Answer

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