Error in Attempting Moving Boundary Fluid System

Question

Recently I was attempting to solve a moving boundary fluid system on mathematica, which I have managed to convert into a coupled PDE-ODE system based on this helpful reference over here.

The equations I am solving are

$$G1[r,t]=m_1 \ln(100*B[r,t])+b1;$$ $$G2[r,t]=m_2 \ln(100*B[r,t])+b2;$$

$$W[r,t]=\frac{p}{q}\frac{\partial (G1[r,t]+G2[r,t])}{\partial r};$$ $$r\frac{\partial N[r,t]}{\partial t}+\frac{\partial N[r,t] W[r,t]}{\partial r}=kB[r,t];$$

$$\frac{\partial}{\partial r}(r N[r,t] \times c \frac{\partial B[r,t]}{\partial r})+k B[r,t] (1-B[r,t])=r N[r,t] (\frac{\partial B[r,t]}{\partial t}+W[r,t] \frac{\partial B[r,t]}{\partial r});$$

Where the maximum value of r, $b[t]$, is solved by the following equation

$$N[b[t],t]\cdot b[t]\cdot(W[b[t],t]-b'[t])=\frac{d}{dt}(\frac{\pi z[t]^2 b[t]}{4});$$

$$\alpha\cdot W[b[t],t] \cdot (W[b[t],t]-b'[t])+(g_0-G1[b[t],t]-G2[b[t],t])-\frac{q\space z[t] \space b'[t]}{p};$$

So in total there are 4 base variables $N[r,t]$, $B[r,t]$, $b[t]$ and $z[t]$.

The approach I took for this problem is to create normalised variable $x=r/b[t]$ to remove the free boundary and use the pdetoode function to convert all PDEs into a massive ode system that NDSolve is able to handle. The code is here:

(*Constants and known functions*)
b1 = 0.0698; m1 = -0.0119; b2 = 0.0526; m2 = -0.0128; g0 = 0.0298;
q = 3*10^-5; p = 3; k = -0.115; c = 5*10^-4; vis = 10^-4; α = 0.001;
g1[x_, t_] := m1 Log[100*B[x, t]] + b1;
g2[x_, t_] := m2 Log[100*B[x, t]] + b2;
w[x_, t_] := Simplify[(D[g1[r/b[t], t] + g2[r/b[t], t], r]*p/q) /. r -> x b[t]]

(*Defining Equations*)
With[{n = n[r/b[t], t], B = B[r/b[t], t], w = w[r/b[t], t]},
 eqnC = Simplify[{(r D[n, t] + D[r*w*n, r] == vis D[ n, r, r] + k*B*r) /. r -> x b[t]}];
 eqnDC = Simplify[{(D[ r n c D[B, r], r] + k B (1 - B) r == r n (D[B, t] + w D[B, r])) /. r -> x b[t]}];]

eqnRimC = Simplify[{n[1, t] b[t] (w[1, t] - b'[t]) == D[ (π z[t]^2 b[t])/4,t]}];
eqnRimM = Simplify[{ α n[1, t] w[1, t] (w[1, t] - b'[t]) + (g0 - g1[1, t] - g2[1, t]) - (q z[t] b'[t])/p == 0}];

(*Setting Initial and Boundary Conditions*)
boundary = 1; σ = 0.2; n0 = 0.4; x0 = 0.001;
nini[x_] := n0*(-x^2) + n0*(1.05);
Bini = 0.42;

With[{n = n[x, t], B = B[x, t]},
 ic = {n == nini[x], B == Bini} /. t -> 0;
 bc = {D[n, x] == 0, D[B, x] == 0} /. x -> x0]

(*Conversion into ODE system*)
points = 6;
domain = {x0, 1};
grid = Array[# &, points, domain];
difforder = 4;
timesolve = 0.03;

tfunc = pdetoode[{n[x, t], B[x, t]}, t, grid, difforder];
removeredundant = #[[2 ;; -1]] &;

odeqnC = Flatten[eqnC // tfunc] // removeredundant;
odeqnDC = Flatten[eqnDC // tfunc] // removeredundant;
odeqnRimC = eqnRimC // tfunc;
odeqnRimM = eqnRimM // tfunc;

int = removeredundant /@ tfunc@ic;
int1 = Append[int, {b[t] == 10} /. t -> 0];
odic = Append[int1, {z[t] == 0} /. t -> 0];
odbc = bc // tfunc;

(*Solve*)
time = 0;
Monitor[{soln, solB, solb, solz} = NDSolveValue[{odeqnC, odeqnDC, odeqnRimC, odeqnRimM, odic, odbc}, {n /@ grid, B /@ grid, b, z}, {t, 0, timesolve}, SolveDelayed -> True, EvaluationMonitor :> (time = t)], time];

However, when I attempt to solve the equations, NDSolve keeps returning me the error that

NDSolveValue::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

My questions are: why would the system be classified by NDSolveValue a differential-algebric system, since every variable has at least one differential term? How could I resolve this issue and obtain my solutions?

Thanks!

EDIT:

It turns out after increasing the grid to a more reasonable value of 500, the initial error disappears but the code stops solving at 5*10^-6 according to the Monitor function.

BACKGROUND INFO*

I was interested in the work of this contributor over here but realised that the assumption made were not realistic. I therefore added a rim equation that is described in detail by Villermaux and Bossa (2010).

In brief I added a cylindrical rim to the edge of the mother droplet, where $z[t]$ is the diameter of the circular cross section of the rim. The first equation describes continuity; the second equation describes conservation of momentum where the effects of surface tension and viscous drag are considered.

Answer 1

… why would the system be classified by NDSolveValue a differential-algebric system, since every variable has at least one differential term?

This isn't quite right, because there's no derivative of t in odbc. Do remember derivatives of x have all disappear after discretization.

So the resulting system is a DAE system. As we know, DAE system is generally harder to solve compared with ODE system, and it's clear that, the DAE solver of NDSolve is having trouble with your system (to be precise, in v12.3.1 it crashes), so let's transform the odbc to something involving derivative of t.

Since the transformation is frequently required when dealing with troublesome DAE system, I've just included a diffbc function here. (To learn more about why we need diffbc, you may want to read this post. ) With diffbc, definition of odbc is modified as follows:

(* Definition of diffbc isn't included in this post,
   please find it in the link above. *)
With[{sf = 1}, odbc = diffbc[t, sf]@bc // tfunc]

To recover the original b.c. from the new odbc, proper i.c. is needed, so we modify the definition of int to

int =tfunc@ic;

Avoiding DAE system isn't enough. Trial-and-error suggests the i.c. z[0] == 0 has lead to a removable singularity at t==0, so we modify it to

With[{zero = 10^-3}, odic = Append[int1, {z[t] == zero} /. t -> 0];]

Still, this isn't the end. With these 3 modifications and points = 40; difforder=2 (it's not necessary, but you can take away SolveDelayed -> True in this case), we

NDSolveValue::ndsz: At t == 0.015003215536343043`, step size is effectively zero; singularity or stiff system suspected.

…encounter ndsz warning around t == 0.015. After observing the solution and the original PDE carefully, I believe the stiffness shows up because $N$ is too close to 0 when $t>0.015$, which leads to another removable singularity. Thus I modify eqnDC to

unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;

With[{n = n[r/b[t], t], B = B[r/b[t], t], w = w[r/b[t], t], ϵ = 10^-3},
 eqnDC = Simplify[{(D[r n c D[B, r], r] + k B (1 - B) r == 
        r unitStepExpand@If[n < ϵ, ϵ, n] (D[B, t] + w D[B, r])) /. 
      r -> x b[t]}];]

Alternatively, eqnDC can be defined as

forcepositive[a_List] := With[{ϵ = 10^-3}, If[# < ϵ, ϵ, #] & /@ a];

With[{n = n[r/b[t], t], B = B[r/b[t], t], w = w[r/b[t], t]}, 
 eqnDC = Simplify[{(D[r n c D[B, r], r] + k B (1 - B) r == 
        r forcepositive[n] (D[B, t] + w D[B, r])) /. r -> x b[t]}];]

Now even timesolve = 0.5; isn't a problem! The following is obtained with points = 100; difforder=2:

tend = solb["Domain"][[1, -1]];

solBtst = rebuild[solB, grid, 2];
solntst = rebuild[soln, grid, 2];

ListLinePlot[#, PlotRange -> All] & /@ {solb, solz} // GraphicsRow

Manipulate[
 Plot[{solBtst[r/solb[t], t], solntst[r/solb[t], t]}, {r, x0 solb[t], solb[t]}, 
  PlotRange -> {{0, solb[tend]}, {0, 0.5}}], {t, 0, tend}]

Remark

1. As to the usage of unitStepExpand: PiecewiseExpand can, as documented, transform a If[…] to mathematically equivalent Piecewise[…]; Simplify`PWToUnitStep can transform Piecewise[…] to mathematically equivalent combination of UnitStep[…]. (This undocumented function is used quite a bit in this site, just have a search if you want to read more examples. ) In other words, unitStepExpand is a function that transforms If[…] to mathematically equivalent combination of UnitStep[…]. I define unitStepExpand because pdetoode can only handle Listable function, but If is not, while UnitStep is Listable. The underlying idea of forcepositive is similar. Another related post can be found here.

2. As to the meaning of last argument of rebuild: as suggested by the source code of rebuild, last argument of rebuild specifies position of time i.e. it specifies where the independent variable representing time (in your case, t) is placed. By default it's 1 i.e. 1st argument of the generated InterpolatingFunction[…] is t, but you've chosen the order …[x, t] i.e. in your code the function with independent variables is e.g. $B(x,t)$ rather than $B(t,x)$, so we set the last argument of rebuild to 2 to keep the order of independent variables consistent.

3. With points = 40 a ListInterpolation::inddp warning pops up, and rebuild[…] fails. It seems that the DiscontinuityProcessing of NDSolve is triggered in this case. The underlying reason isn't immediately clear to me, but adding Method -> {"DiscontinuityProcessing" -> False} to NDSolveValue, or simply choosing a denser grid e.g. points = 100 resolves the problem.