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

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.

• When I run your code I do not get that output. ![Output of what happens when I run your code](i.stack.imgur.com/P8hRN.png) I think you have defined variables earlier that has caused this error. Either that or the matrix is singular. May 28, 2022 at 4:22
• 1. Which version are you in? 2. Adding background information of the system to the question might help. (Is it from certain literature? Or you deduce it yourself? If so, how you decuce it? ) 3. I suggest using With to rewrite the code (to something like With[{u=u[t,x]},D[u,t]+D[u,x]==0]) to make it easier to read. 3. The highest differential order of c in r direction is 2, and that of h is 1, usually this suggests you need 2 b.c. for c and 1 b.c. for h, are you sure your b.c.s are correct? May 28, 2022 at 4:39
• Hi, thanks for the comment! Sorry for being unclear in the post, I am very new to StackExchange. I am using Mathematica 12. 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. May 28, 2022 at 5:01
• 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 May 28, 2022 at 5:11
• James, could you share a literature reference to the model? May 28, 2022 at 13:00

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.

• Wow, thanks so much!!! Also, when I run the code I obtain the error "warning: boundary and initial conditions are inconsistent", but it does not seem to affect the solving of the equations. Can you help explain this? Thanks again : ) May 28, 2022 at 11:37
• @JamesHe Check the update. May 28, 2022 at 11:52
• Okay, thank you! May 28, 2022 at 12:53
• (+1) I always appreciate your careful distinction pointing out in what cases to use TPG vs. FEM. May 28, 2022 at 13:01