I’m trying to model the evaporation of a binary droplet in a square well by solving a couple of PDEs for its height profile and composition. I can manage it for a single liquid pretty well, but when I introduce the PDE for composition NDSolve quickly finds a singularity in the composition variable (importantly the singularity moves slightly if I change the number of points). I can’t seem to find any way to stabilise the system, can anyone help me? Is there something fundamental about NDSolve that I’m missing?

(I asked a similar question a while ago but I’ve restructured my code to make it easier to share here)

The equations are:

$ $ \frac{\partial h}{\partial t} = – C \frac{\partial}{\partial x}\Big[\frac{1}{3} \frac{\partial^3h}{\partial x^3}h^3 + \frac{B1}{2} \frac{\partial X}{\partial x}\Big] – (1-AX) $ $ $ $ \frac{\partial X}{\partial t} = -\frac{C}{3} h^2 \frac{\partial^3h}{\partial x^3}\frac{\partial X}{\partial x} – \frac{hB1}{2} \Big(\frac{\partial X}{\partial x}\Big)^2 + \frac{A}{h}X(1-X) $ $

with height $ h$ and composition $ X$ .

And the boundary conditions are

$ h(t)=1$ and $ X(t)=X(t=0)$ at $ x=1$ , $ \frac{\partial h}{\partial x}=\frac{\partial X}{\partial x}=0$ at $ x=0$ , and $ \frac{\partial^3h}{\partial x^3}=0$ at $ x=0$ and $ x=1$ .

I discretise all this using finite differences.

This is what I get for a single liquid droplet (when I suppress equation 2) and I expect something qualitatively similar for a binary droplet.

Here is a minimum replicable example:

` N1 = 100; dx = 1/N1; (*Discretise*) C1 = 1; (*parameter C*) a = 1.5; b = 1 - a; hInitial = Table[a + b i^2 dx^2, {i, 0, N1}]; (*Initial height profile*) E1 = 1;(*evaporation on/off switch*) Cc = C1/(24 * dx^3); A1 = 0.01; B1 = 1; (*Control constants in the equations above*) CcX = (C1 B1)/(2 dx^2); CcX2 = C1/(6 dx^2); X0 = 0.4; XInitial = Table[X0, {i, 0, N1}]; (*Initial composition ratio*) Sc = 1; (*Composition on/off switch*) TGap = 0.085;(*for spacing curves in the plot*) dv[v_List] := With[{h = Take[v, Length[v]/2 - 1], hN = v[[(Length[v]/2)]], X1 = v[[Length[v]/2 + 1 ;; Length[v] - 1]], XN = v[[-1]]}, With[{ dh1 = ListCorrelate[{0, 0, 1 , 1, 0}, #] &, (*derivatives are discretised as finite differences. ListCorrelate achieves this.*) dh2 = ListCorrelate[{0, 1, 1, 0, 0}, #] &, dh3 = ListCorrelate[{0, -1, 3, -3, 1}, #] &, dh4 = ListCorrelate[{-1, 3, -3, 1, 0}, #] &, hi = ListCorrelate[{0, 0, 1, 0, 0}, #] &, dh5 = ListCorrelate[{-2, 1, 0, -1, 2}, #] &, dx1 = ListCorrelate[{0, -1, 1}, #] &, dx2 = ListCorrelate[{-1, 1, 0}, #] &, Xi = ListCorrelate[{0, 1, 0}, #] & }, Flatten@{ -Cc*(dh1[#1]^3*dh3[#1] - dh2[#1]^3*dh4[#1]) - CcX*(dh1[#1]^2*dx1[#2] - dh2[#1]^2*dx2[#2]) - E1*(1 - A1 Xi[#2]), (*ODE for heights*) 0,(*height is pinned at edge, (1, 1)*) Sc*(-CcX2* hi[#1]^2*dh5[#1]*dx2[#2] - CcX/C1* hi[#1]*dx2[#2]^2 + A1 *(hi[#1])^-1*Xi[#2]*(1 - Xi[#2])), (*ODE for composition*) 0 (*composition ratio is fixed at edge, XN = 0.4*) } &[ Flatten@Join[{h[[3]]}, {h[[2]]}, {h}, {1}, {3 - 3 h[[-1]] + h[[-2]]}], Flatten@Join[{X1[[2]]}, {X1}, {XN}]]] (*this contains the boundary conditions*) ]; v0 = Join[hInitial, XInitial]; (*Initial conditions*) system2d = NDSolve[{v'[t] == dv[v[t]], v[0] == v0, WhenEvent[ Min@Table[ Take[values[[tt]], Length[values[[tt]]]/2], {tt, 1, Length[values]}] < 0, "StopIntegration"](*WhenEvent to stop integration at touchdown*)}, v, {t, 0, 2}][[1, 1, 2]]; Needs@"DifferentialEquations`InterpolatingFunctionAnatomy`"; values = InterpolatingFunctionValuesOnGrid@system2d; times = Flatten@InterpolatingFunctionGrid@system2d; T1 = 0; timesteps = {}; For[i = 1, i < Length[times], i++, If[times[[i]] > T1, {T1 = T1 + TGap, timesteps = Append[timesteps, {i, times[[i]]}]}]]; finaltime = {{-1}}; firstandfinaltime = {{1}, {-1}}; heights = Table[Take[values[[tt]], Length[values[[tt]]]/2], {tt, 1, Length[values]}]; comps = Table[ Take[values[[tt]], -(Length[values[[tt]]]/2)], {tt, 1, Length[values]}]; Show[Table[ ListPlot[Table[{i*dx, heights[[ttt[[1]], i + 1]]}, {i, 0, N1}], PlotRange -> {{0, 1}, {0, 1.05*a}}, Joined -> True, AxesLabel -> {x, h}], {ttt, finaltime}]] Show[Table[ ListPlot[Table[{i*dx, comps[[ttt[[1]], i + 1]]}, {i, 0, N1}], PlotRange -> {{0, 1}, {0, 1}}, Joined -> True, AxesLabel ->{x, X}], {ttt, finaltime}]] `

Thank you for any help you can offer me!