# How do I couple an 1D PDE to a 2D Laplace equation?

Recently I've come across the following system governing the spreading of an evaporating droplet. The height of the droplet $$h$$ is defined by the following equation

$$\frac{\partial h}{\partial t}=-\frac{1}{r}\frac{\partial }{\partial r}\left(- \frac{r h^3}{3\mu}\frac{\partial p}{\partial r}\right)-j(r,t)$$

Which is a clear, solvable 1D Partial-Differential Equation solvable by the Method of Lines.

The evaporation rate, $$j(r,t)$$, however is given as a solution to the following 2D Laplace Equation:

$$\nabla^2 c(r,z,t)=0$$ $$j(r,t)=\alpha \frac{\partial c(r,0,t)}{\partial z}$$

With boundary conditions $$c(r,0,t)=h(r,t) \ for \ r < R(t)$$ $$\frac{\partial c(r,0,t)}{\partial z}=0 \ for \ r> R(t)$$ $$c(r,z,t)=0 \ for \ r\rightarrow \infty, z \rightarrow \infty$$

While I'm able to solve the spreading of the droplet, as well as the evaporation due to an abitrary concentration field independently, I have no idea how to couple the two solutions together. The two way coupling is particularly annoying as

• The solutions of the height field is used as the boundary condition of the Laplace equation;
• The solution of the Laplace equation gives the evaporation rate which affects the height field.

I'll appreciate any advice to go forward on this problem. Below are my current code are both the spreading of the Droplet and the Evaporation of an arbitrary concentration field. Thanks!

Note that the pdetoode package by xzczd is used. For sake of clarity of code that part is not included.

(*Definitions*)

(*Constants*)
\[Mu] = 2.75*10^-6 (*mN mm^-2 s*);(*Viscosity of water*)
\[Gamma] = 0.072(*mN mm^-1*);(*Air-water surface tension*)
\[Beta] = 1(*mm s^-1*);(*Slip constant*)

(*Grid specification*)
difforder = 2;
lb = 1/100; rb = 1; points = 100;
grid = Array[# &, points, {lb, rb}];

unitStepExpand = SimplifyPWToUnitStep@PiecewiseExpand@# &;

With[{h = h[r/R[t] , t], q = q[r/R[t], t], p = p[r/R[t], t], R = R[t]},

(*Equations neglecting evaporation*)
EqnP = Simplify[
p == \[Gamma] (-D[h, r, r] -
unitStepExpand@If[x < 1/20, D[h, r, r], D[h, r]/r])] /.
r -> x R;
EqnQ = Simplify[q == (r h^3)/(3 \[Mu]) (-D[p, r])] /. r -> x R;
EqnC = Simplify[ D[h, t] == -(1/r) D[ q, r]] /. r -> x R;
EqnR = D[R, t] == \[Beta] (-D[h, r]^3 - 1) /. r -> R;]

(*Boundary conditions*)
dr = R[t] lb;

EqnBC1 = {D[r h[lb, t] dr, t] + r q[lb, t] == 0} /. r -> lb R[t];
EqnBC2 = {Derivative[1, 0][h][lb, t] == 0};
EqnBC3 = {h[rb, t] == 0};

(*Initial conditions*)
EqnIC = {h[x, 0] == (1 - x^2)};

(*Equation Discretisation*)
removeredundant = #[[3 ;; -2]] &;

tfunc = pdetoode[{p, q, h}[x, t], t, grid, difforder];
odemid = Map[tfunc, {EqnP, EqnQ}, {2}];
odeC = Block[{p, q}, Set @@@ odemid; tfunc@EqnC] // removeredundant;

odeBC1 = Block[{p, q}, Set @@@ odemid; tfunc@EqnBC1];
odeBC2 = With[{sf = 100}, diffbc[t, sf]@EqnBC2 // tfunc];
odeBC3 = With[{sf = 100}, diffbc[t, sf]@EqnBC3 // tfunc];

odeR = Block[{p, q}, Set @@@ odemid; tfunc@EqnR];

odeIC = Append[EqnIC // tfunc, R[0] == 10];

(*Solving the equations*)
var = {h /@ grid, R};

Monitor[sollst =
NDSolveValue[{odeC, odeR, odeBC1, odeBC2, odeBC3, odeIC},
var, {t, 0, 6}, EvaluationMonitor :> (time = t),
Method -> {"EquationSimplification" -> "Solve"},
WorkingPrecision -> 20, MaxSteps -> Infinity], time];

{hsol} = rebuild[#, grid, 2] & /@ {sollst[[1]]}
Rsol = sollst[[2]]

(*Visualising Solutions*)
Manipulate[
Plot[hsol[Abs[r]/Rsol[t], t], {r, 1/10, Rsol[t]},
PlotRange -> {{0, 10}, {0, 3}}], {t, 0, 1}]


Evaporation Code

Needs["NDSolveFEM"]

(*Mesh Generation*)
L = 10;
H = 5;

Line[{{0}, {L}}], <|"Alignment" -> "Left", "ElementCount" -> 40|>]
Line[{{0}, {H}}], <|"Alignment" -> "Left", "ElementCount" -> 20|>]

mesh = ElementMeshRegionProduct[l1, h1]

(*Arbitrary height function*)
h[x_] := 1 - x^2;

(*Solving Equations*)
EqnC = 1/x D[x Derivative[1, 0][c][x, y], x] +
Derivative[0, 2][c][x, y] ==
NeumannValue[0, y == 0 && Abs[x] > 1] + NeumannValue[0, y == H] +
NeumannValue[0, Abs[x] == L];
EqnB = DirichletCondition[c[x, y] == h[x], Abs[x] <= 1 && y == 0];
EqnB2 = DirichletCondition[c[x, y] == 0, y == H];
EqnB3 = DirichletCondition[c[x, y] == 0, Abs[x] == L];

csol = NDSolveValue[{EqnC, EqnB, EqnB2, EqnB3},
c, {x, y} \[Element] mesh]
j[r_] := \[Alpha] (csol[r, 0] - csol[r, 0.001])/0.001

(*Visualising Results*)
Plot[j[x] /. \[Alpha] -> 1, {x, 0, 10}, PlotRange -> All]

• Do you mean conjugate problem on time dependent mesh like in my answer on physics.stackexchange.com/questions/509623/… ? Commented Apr 6 at 14:21
• Hi! In my case the mesh motion isn't that big of an issue - it's more that the solution to the 1D PDE need to be applied to the boundary of the 2D Laplace equation at every timestep.
– FLP
Commented Apr 6 at 14:36
• I mean that problem has moving interface. Commented Apr 6 at 14:43
• Oh! The moving interface is handled using a height function h which evolves via the equation above. The edge of the droplet is modelled using R(t) which evolves via an ODE. I therefore solved this "moving interface" problem by converting it into a 1D PDE.
– FLP
Commented Apr 6 at 14:46
• Thanks, see my answer. :) Commented Apr 9 at 19:25

This problem can be solve with FEM and FDM as well. First we show FDM code since pdetoode is implementation of FDM. We copy this function as well from here (thanks to xzczd)

Clear[fdd, pdetoode, tooderule, pdetoae, diffbc, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolveFiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_ : False] :=
pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer,
periodic : True | False | {(True | False) ..} : False] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
Piecewise[{{1, PossibleZeroQ[# - #2[[1]]]}, {-1,
PossibleZeroQ[# - #2[[-1]]]}}, All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};

fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o,
PeriodicInterpolation -> periodic]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]

pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=
Module[{t},
Function[
pde, #[pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]]

diffbc[rst__][a : _List | _Equal] := diffbc[rst] /@ a
diffbc[dvar : {t_, order_} | (t_) .., sf_ : 0][a_] /; sf =!= t :=
sf a + D[a, dvar]

rebuild[funcarray_, grid_?VectorQ, timeposition_ : 1] :=
rebuild[funcarray, {grid}, timeposition]

rebuild[funcarray_, grid_, timeposition_?Negative] :=
rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]

rebuild[funcarray_, grid_, timeposition_ : 1] /;
Dimensions@funcarray === Length /@ grid :=
With[{depth = Length@grid},
ListInterpolation[
Transpose[
Map[DeveloperToPackedArray@#["ValuesOnGrid"] &, #, {depth}],
Append[Delete[Range[depth + 1], timeposition], timeposition]],
Insert[grid, Flatten[#][[1]]["Coordinates"][[1]],
timeposition]] &@funcarray]

(*Definitions*)(*Constants*)\[Mu] =
2.75*10^-6 (*mN mm^-2 s*);(*Viscosity of water*)
\[Gamma] = 0.072(*mN mm^-1*);(*Air-water surface tension*)
\[Beta] = 1(*mm s^-1*);(*Slip constant*)

(*Grid specification*)
difforder = 4;
lb = 1/100; rb = 1; points = 20;
grid = Array[# &, points, {lb, rb}];

unitStepExpand = SimplifyPWToUnitStep@PiecewiseExpand@# &;

With[{h = h[r/R[t], t], q = q[r/R[t], t], p = p[r/R[t], t],
R = R[t]},(*Equations neglecting evaporation*)
EqnP = Simplify[
p == \[Gamma] (-D[h, r, r] -
unitStepExpand@If[x < 1/20, D[h, r, r], D[h, r]/r])] /.
r -> x R;
EqnQ = Simplify[q == (r h^3)/(3 \[Mu]) (-D[p, r])] /. r -> x R;
EqnC = Simplify[D[h, t] == -(1/r) D[q, r]] /. r -> x R;
EqnR = D[R, t] == \[Beta] (-D[h, r]^3 - 1) /. r -> R;]

(*Boundary conditions*)
dr = R[t] lb;

EqnBC1 = {D[r h[lb, t] dr, t] + r q[lb, t] == 0} /. r -> lb R[t];
EqnBC2 = {Derivative[1, 0][h][lb, t] == 0};
EqnBC3 = {h[rb, t] == 0};

(*Initial conditions*)
EqnIC = {h[x, 0] == (1 - x^2)};

(*Equation Discretisation*)
removeredundant = #[[3 ;; -2]] &;

tfunc = pdetoode[{p, q, h}[x, t], t, grid, difforder];
odemid = Map[tfunc, {EqnP, EqnQ}, {2}];
odeC = Block[{p, q}, Set @@@ odemid; tfunc@EqnC] // removeredundant;

odeBC1 = Block[{p, q}, Set @@@ odemid; tfunc@EqnBC1];
odeBC2 = With[{sf = 100}, diffbc[t, sf]@EqnBC2 // tfunc];
odeBC3 = With[{sf = 100}, diffbc[t, sf]@EqnBC3 // tfunc];

odeR = Block[{p, q}, Set @@@ odemid; tfunc@EqnR];

odeIC = Append[EqnIC // tfunc, R[0] == 10];


FDM solution of the Laplace equation in general form

delx = Differences[grid][[1]];
XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom, pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]]}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{domain, pts,
difforder} = {{{lb, 5 points delx + lb }, {0, H}}, {5 points + 1,
20}, 4};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid,
difforder]; boundaries = Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1]; bot1 =
Take[bot, points]; bot2 = Complement[bot, bot1];
u = Table[uu[i], {i, nx ny}]; X =
sgrid[[All, 1]]; eqs = (dx2 + dy2) . u + (dx . u)/X;
eqs[[left]] = (dx . u)[[left]]; eqs[[right]] = u[[right]];
eqs[[top]] = u[[top]];
varh = h[#][t] & /@ grid; eqs[[bot1]] = u[[bot1]] - varh;
eqs[[bot2]] = (dy . u)[[bot2]];

{vec, mat} = CoefficientArrays[eqs, u];

inv = Inverse[mat // N];

U = -inv . vec;

ju = (dy . U)[[bot1]]; ju1 = (ju // removeredundant)/R[t];

odeC1 = Thread[odeC[[All, 1]] == ju1];


Please note, that ju is the evaporation rate expressed as a function of h[t] on the grid. We use odeC1 to compute solution as follows

Monitor[sollst =
NDSolveValue[{odeC1, odeR, odeBC1, odeBC2, odeBC3,
odeIC}, {h /@ grid, R}, {t, 0, 6},
EvaluationMonitor :> (time = t),
Method -> {"EquationSimplification" -> "Solve"}], time];

{hsol} = rebuild[#, grid, 2] & /@ {sollst[[1]]};
Rsol = sollst[[2]];


Animation

frames =
Table[Plot[hsol[Abs[r]/Rsol[t], t], {r, 1/10, Rsol[t]},
PlotRange -> {{0, 10}, {0, 1}}], {t, 0, 6, .05}];


Update 1. FEM code. Using FEM we can construct function analogue to ju as follows

Needs["NDSolveFEM"]
L = 5;
H = 5;
mesh = ToElementMesh[Rectangle[{10^-2, 0}, {L, H}],
MaxCellMeasure -> .01];
jF[h_, g_] :=
Module[{c, x, y, csol, hh}, hh = Interpolation[Transpose[{g, h}]];
csol =
NDSolveValue[{1/x  D[x  Derivative[1, 0][c][x, y], x] +
Derivative[0, 2][c][x, y] == 0,
DirichletCondition[c[x, y] == hh[x], x <= 1 && y == 0],
DirichletCondition[c[x, y] == 0, y == H],
DirichletCondition[c[x, y] == 0, x == L]},
c, {x, y} \[Element] mesh]; -Derivative[0, 1][csol][#, 0] & /@ g]


Now we can compare FEM and FDM code using numerical solution shown above

var = Join[h /@ grid, {R}]; sol =
NDSolve[{odeC1, odeR, odeBC1, odeBC2, odeBC3, odeIC}, var, {t, 0, 6},
Method -> {"EquationSimplification" -> "Solve"}];

Table[ListLinePlot[{-ju /. sol[[1]], jF[varh /. sol[[1]], grid]},
PlotStyle -> {{Blue}, {Red, Dashed}}, PlotLabel -> Row[{"t = ", t}],
PlotLegends -> {"FDM", "FEM"}], {t, 0, 6, .5}]


As we can see from picture above there are discrepancies FDM and FEM function due to difference in option DifferenceOrder which is 4 for FDM differentiation matrices and 2 for FEM matrices.

• Thank you so much!!! :D
– FLP
Commented Apr 9 at 23:54
• You are welcome! Commented Apr 10 at 0:09
• Hi! Regarding the Laplace Equation FDM implementation - is there a reason why the Dirichlet boundary condition is implemented on the right boundary while the Neumann boundary condition is implemented on the left (instead of implementing both Dirichlet and Neumann on both)? thanks!!!
– FLP
Commented Apr 11 at 14:16
• @FLP Left boundary is axis of symmetry by definition, therefore we should put zero Neumann at x==0 ( x==lb` in numerical model), while right boundary is "Infinity" in numerical model, therefore we should put zero Dirichlet. Please note, that evaporation rate in your model is given by $j=- \alpha (\nabla c.\vec {n})$, where $\vec {n}$ is unit normal, and in my code above we use $\frac{\partial c}{\partial y}$ only as for disk. But actually it is not a disk. Commented Apr 11 at 14:59