Droplet of water-alcohol mixture spreading due to evaporation-caused surface tension gradient

Question

I am trying to solve a coupled system of PDEs for 2 functions h[t,r] and c[t,r], with initial conditions of h[0,r] as a normal distribution over r and c[0,r] set to a fixed value of 0.42 (c[r,t] is supposed to always be bounded between 0 and 1), and fixed boundary values for h[t,r] at two points r = 0.1 and r = -0.1 (solve domain boundaries). However, when I put the equations and boundary conditions into NDSolve, I get the warning

Dif = 10^(-5);
cinitial = 0.42;
μ = 0.10;
H = 0.003;
k = 0.001;
boundary = 0.1;
u = H*0.025*D[c[t, r], r]/(μ*c[t, r]);
NDSolve[{k r (1 - c[t, r]) c[t, r] + h[t, r] Dif D[c[t, r], r] + 
    r D[h[t, r], r] Dif D[c[t, r], r] + 
        r h[t, r] Dif D[c[t, r], {r, 2}] == 
      r h[t, r] (D[c[t, r], t] + u D[c[t, r], r]), 
  r D[h[t, r], t] + (r u D[h[t, r], r] + D[u, r]*h[t, r]) == k c[t, r] r, 
    h[0, r] == E^(-(r^2/(2 0.02^2)))/(0.02 10 Sqrt[2 π]), 
    c[0, r] == cinitial, 
    h[t, boundary] == E^(-(boundary^2/(2 0.02^2)))/(
        0.02 10 Sqrt[2 π]), 
    h[t, -boundary] == E^(-(boundary^2/(2 0.02^2)))/(
        0.02 10 Sqrt[2 π])}, {h[t, r], c[t, r]}, {t, 0, 
    10}, {r} ∈ ImplicitRegion[r^2 <= boundary^2, {r}]]

NDSolve::femcnsd: The PDE coefficient 0. -0.001 r c[r]+0.001 r c[r]^2-0.00001 h[r] c'[r]+(0.00075 r h[r] c'[r]^2)/c[r]-0.00001 r c'[r] h'[r] does not evaluate to a numeric scalar at the coordinate {-0.1}; it evaluated to Indeterminate instead.

I am really not sure if the cause is my equations or NDSolve. I have checked the equations and boundary conditions and was unable to find any errors, and I do not have enough understanding of NDSolve and Mathematica. Is anyone able to identify the problem? Thanks!

I am using Mathematica 12.

---

Update: background info

The equations describe a droplet of water-alcohol mixture spreading on a surface due to a surface tension gradient which is the result of evaporation (concentration at the edge of droplet decreases faster, surface tension ). The first equation is the diffusion-convection equation for alcohol transport in the droplet, while the second equation is the continuity equation accounting for evaporative loss. u is the radial velocity of each fluid element in the drop.

The model is 1D axisymmetric. The initial condition for h is defined to be non-zero for all values of r to make the solving easier, or an additional moving boundary will have to be added into the equation for the radius of spreading of the droplet. The boundary conditions for h are equivalent to setting the bc that the height of the droplet at a radial coordinate approaching infinity to be 0. I think it is also possible to include the condition D[c[r,t],r] at r = 0.1 and -0.1 (again at the radial boundaries) to be 0, if required.

Answer 1

At least 3 issues here.

1. You're handling the cylindrical coordinates in an unsual way i.e. setting domain of definition of r to $r\in(-\text{boundary},\text{boundary})$. At least in v12.3.1 this leads to LinearSolve::parpiv, which isn't too surprising, because there's a removable singularity at $r=0$. This isn't hard to fix, we just need to change the domain of definition to $r\in(0,\text{boundary})$. (In the code below I have to set the left boundary to a eps that's close to 0. Actually if one chooses FiniteElement method for spatial discretization, the singularity at $r=0$ will be automatially handled. I don't choose it because the old good TensorProduct method is still generally a better choice compared with FiniteElement method for time-dependent PDE in regular domain. )

2. After transforming to the required standard form of NDSolve, h term appears at denominator, so h that's too close to 0 will cause trouble, we need to modify the i.c. to a milder one. I'll archieve this by adding a small constant (epsic below) to the i.c.

3. The troublesome feature of hyperbolic conservation law seems to show up in the solution. This might be circumvented by transforming the equation system to something close to conservation form, but it's a bit too hard for me without knowing more about the background info of the system, so I'll simply add artificial viscosity (vis below).

The following is the solution.

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}
Dif = 1/10^5;
cinitial = 0.42;
μ = 0.1;
H = 0.003;
k = 0.001;
boundary = 0.1;
u = (H 0.025 D[c[t, r], r])/(μ c[t, r]);
eps = boundary/100;
vis = 10^-7;
epsic = 10^-2;
hic = 1/(E^(r^2/(2 0.02^2)) (0.02 10 Sqrt[2 π])) + epsic;
With[{h = h[t, r], c = c[t, r]}, 
  eq = {k r (1 - c) c + Dif D[D[c, r] h r, r] == r h (D[c, t] + u D[c, r]), 
        r D[h, t] + (r u D[h, r] + D[u, r] h) == vis D[h, r, r] + k c r};
  ic = {h == hic, c == cinitial} /. t -> 0;
  bc = {{D[c, r] == 0, D[h, r] == 0} /. r -> eps, 
        {D[c, r] == 0, h == hic} /. r -> boundary}];
points = 400; difforder = 2;
{tsth, tstc} = 
  NDSolveValue[{eq, ic, bc}, {h, c}, {t, 0, 10}, {r, eps, boundary}, 
   Method -> mol[points, difforder]];
(* Timing:1.4 second *)

NDSolve spits out ibcinc warning, but it's not too big a problem in this case. If you don't like it, change the b.c. D[h, r] == 0 to D[h, r] == D[hic, r] and set a denser grid with e.g. points = 1000. If you want to learn more about ibcinc warning, read the following posts:

Boundary condition with spatial derivative is ignored by NDSolve

2D Heat equation: inconsistent boundary and initial conditions

Finally let's visualize:

Plot[#, {r, eps, boundary}, PlotRange -> All, PlotStyle -> Dashed] & /@ 
  Transpose@Table[{tsth[#, r], tstc[#, r]} &[t], {t, 0, 10, 10/5}] // GraphicsRow

You may adjust epsic, vis, etc. further to see how the solution changes.

---

If you prefer FEM:

molfem[measure_: Automatic] := {"MethodOfLines", 
   "SpatialDiscretization" -> {"FiniteElement", 
     "MeshOptions" -> MaxCellMeasure -> measure}};
{tsthfem, tstcfem} = 
  NDSolveValue[{eq, ic, bc[[-1, -1]]}, {h, c}, {t, 0, 10}, {r, 0, boundary}, 
    Method -> molfem[boundary/200]]
(* Timing: 481.155 second *)

The discrepency between 2 approaches is ignorable.

I've only explicitly set one b.c. bc[[-1, -1]] because the default setting of FiniteElement method is zero Neumann value.