# NDSolve does not go past initialisation stage when attempting to solve for Fluid System

BACKGROUND INFO

Following my attempt to solve a fluid system here (I am very much indebted to xzczd for his gracious help), I realised that a few of my equations are physically incorrect since they produces physically inaccurate solutions under the limit that k approaches zero. I have therefore edited the equations to and solved them with the same method as described in my previous post.

Essentially, I included the pressure and inertial terms for the momentum equation, making the assumption that the pressure acting on the fluid is equal to the Laplace pressure acting on the surface as height is very small. For the rim momentum equations, I included a previously ignored term that calculates the momentum of the rim itself. (Details of the rim are included in the previous post).

EQUATIONS

The equations are

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

$$\alpha N[r,t] (\frac{\partial W[r,t]}{\partial t}+W[r,t]\frac{\partial W[r,t]}{\partial r})=-\frac{G1[r,t]^2 N[r,t] }{G2[r,t]} \frac{\partial^3 N[r,t]}{\partial r^3}+\frac{\partial (G1[r,t]+G2[r,t])}{\partial r}-\frac{q W[r,t]}{p};$$

$$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}=\frac{b'[t]}{b[t]}\cdot\frac{d}{dt}(\frac{\pi z[t]^2 b[t]}{4});$$

CODE

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

(*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.000115; c = 5*10^-4; vis =
10^-4; \[Alpha] = 0.001;
g1[x_, t_] := m1 Log[100*B[x, t]] + b1;
g2[x_, t_] := m2 Log[100*B[x, t]] + b2;

unitStepExpand = SimplifyPWToUnitStep@PiecewiseExpand@# &;

(*Defining Equations*)
With[{n = n[r/b[t], t], B = B[r/b[t], t], w = w[r/b[t], t],
g1 = g1[r/b[t], t], g2 = g2[r/b[t], t], \[Epsilon] = 10^-3},
eqnC = Simplify[{(r D[n, t] + D[r*w*n, r] == k*B*r) /. r -> x b[t]}];
eqnM = Simplify[{\[Alpha] n (D[w, t] + w D[w, r]) == -(g1^2/g2) n D[
n, r, r, r] + D[g1 + g2, r] - (q*w/p)}] /. r -> x b[t];
eqnDC =
Simplify[{(D[r n c D[B, r], r] + k B (1 - B) r ==
r unitStepExpand@
If[n < \[Epsilon], \[Epsilon], 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[(\[Pi] z[t]^2 b[t])/4, t]}];
eqnRimM =
Simplify[{\[Alpha] n[1, t] w[1, t] (w[1, t] - b'[t]) + (g0 -
g1[1, t] - g2[1, t]) - (q z[t] b'[t])/p ==
b'[t]/b[t] D[\[Pi] z[t]^2/4 b[t]]}];

(*Setting Initial and Boundary Conditions*)
boundary = 1; \[Sigma] = 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], w = w[x, t]},
ic = {n == nini[x], B == Bini,
w == 0.001 (*computational zero to prevent discontinuities*)} /.
t -> 0;
bc = {D[n, x] == 0, D[B, x] == 0, D[w, x] == 0} /. x -> x0]

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

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

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

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

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


CURRENT PROBLEM

When attempting to solve for the system however, NDSolve leaves a very strange issue. NDSolve returns no errors, nor does it quit from running, but from the Monitor, the timestep of solving is constantly stuck at 0. after 10 minutes of waiting. I interpreted this error that NDSolve does not go past the initialisation stage, but I have no idea how I could possibly go about debugging this code considering that I have no idea of the error at all.

Is there any way if I could resolve this issue, or is the system just too complex for Mathematica to handle? Thanks!

• 1. Your code involves typo, the bracket doesn't match. 2. Your eqnC, eqnM and eqnDC are obviously different from that in $\LaTeX$ code, which one is correct? 3. You make that simple mistake once again. Please, double-check your post. Commented Jun 27, 2022 at 11:32
• (-1) Will retract my downvote after you correct all the simple mistakes. Commented Jun 27, 2022 at 11:38
• You're making code block in a strange way. Usually we just type Ctrl+K, or add  before and after codes. See this for more info: mathematica.stackexchange.com/editing-help#code Commented Jun 27, 2022 at 12:30
• Sorry I just spotted another few mistakes in the code! Thanks for the suggestion and I'll get back to you in a min.
– FLP
Commented Jun 27, 2022 at 12:35
• Yup, now it should be fully edited since its copied directly from my notebook.
– FLP
Commented Jun 27, 2022 at 12:39