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][1] (thanks to xzczd)
Clear[fdd, pdetoode, tooderule, pdetoae, diffbc, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@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,
MapThread[
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[Developer`ToPackedArray@#["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 = Simplify`PWToUnitStep@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] :=
MapThread[
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[NDSolve`FiniteDifferenceDerivative[#, 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}];
[![Figure 1][2]][2]
**Update 1. FEM code.**
Using FEM we can construct function analogue to `ju` as follows
Needs["NDSolve`FEM`"]
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}]
[![Figure 2][3]][3]
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.
[1]: https://mathematica.stackexchange.com/questions/127980/dynamic-euler-bernoulli-beam-equation/127997#127997
[2]: https://i.sstatic.net/j9Ve5mFd.gif
[3]: https://i.sstatic.net/9efWnlKN.png