# Error in Attempting Moving Boundary Fluid System

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.

New contributor
FLP is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• The issue with setting a boundary condition like 'D[B,x]==0/.x-> boundary' is that it implicitly assumes that there is no velocity inflow (due to the definition of w) into the rim, which is quite unlikely. So I thought that if NDsolve is able to figure out the solution without such an assumption it would be a more physically accurate solution.
– FLP
Jun 22 at 12:59
• Well I think a second order boundary condition holds promise at x = boundary, currently Im setting D[B,x,x], D[n,x,x] ==0 at x-> boundary and it is able to solve past x = 0.1.
– FLP
Jun 22 at 13:41
• Do you still encounter the NDSolveValue::ivres warning after correcting the aforementioned mistakes? If so, please add your version info; if not, please update the question a bit. (In v12.3.1 NDSolve crashes. ) 2 days ago

… 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,
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 = SimplifyPWToUnitStep@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[…]; SimplifyPWToUnitStep 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.

• Thank you so much for your solution! I'll see if I can find any boundary conditions to resolve the stiffness issue :D
– FLP
2 days ago
• @FLP Aha, I made it. See my update. 2 days ago
• Wow, theres really a lot to learn... @xzczd could you explain what the function unitStepExpand does? I have never saw it before and it seems that a cursory search on Google reveals very little info about it.
– FLP
2 days ago
• Thanks for the reply! In addition, why is there a 2 at the end of the rebuild function? I tried the rebuild function normally (without the 2 at the end of the input) and obtained rather weird results.
– FLP
2 days ago
• When I tried to solve the equation beyond time=0.5, the rebuild function strangely returns me the error that the point 0.5429501840077648' in dimension 2 is duplicated`.
– FLP
2 days ago