Solution blows up when using Runge-Kutta to solve simultaneous ODEs for liquid film equations

Question

This is a continuation of a problem I asked over at physics exchange and math exchange. Basically I have two ODEs that I am solving in order to calculate the radial and tangential velocity of liquid dispensed on the center of a disk rotating at rate of $\omega$. $u_m$ is the average flow velocity in the $r$-direction, and the $v_m$ is the average flow velocity in the $\theta$-direction. They are defined as

$$u_m(r)=\frac{1}{h(r)}\int_0^{h(r)}u(r,z)\,dz$$

$$v_m(r)=\frac{1}{h(r)}\int_0^{h(r)}v(r,z)\,dz$$

Where $h(r)$ is the height of the liquid film on the spinning disk. Using the continuity equation, the height of the film is

$$h(r)=\frac{Q}{2\pi\,r\,u_m}$$

Where $Q$ is the flow rate of liquid onto the spinning disk. Using appropriate boundary conditions of no-slip at $z=0$ and no shear at $z=h(r)$, the velocities $u$ and $v$ are:

$$u(r,z)=3u_m\left[ \frac{z}{h}-\frac{1}{2}\left(\frac{z}{h} \right)^2 \right]$$

$$v(r,z)=r\,\omega+\frac{3}{2}(r\,\omega- v_m)\left[\left(\frac{z}{h} \right)^2 -2 \frac{z}{h}\right]$$

Then taking the Navier Stokes equations in cylindrical coordinates and assuming that $\partial P/\partial r=0$, the remaining relevant terms are then

$$u\frac{\partial u}{\partial r}-\frac{v^2}{r}=\nu \left(\frac{\partial^2 u}{\partial z^2} \right)$$

$$u\frac{\partial v}{\partial r}+\frac{vu}{r}=\nu \left(\frac{\partial^2 v}{\partial z^2} \right)$$

Then I substitute $u(r,z)$ and $v(r,z)$ into the N-S equations, integrate them from $0$ to $h(r)$ in the $z$-direction, and then substitute $h(r)$ in terms of $Q$ using the equation above. This then gives me two coupled ODEs where $u_m$ and $v_m$ are the dependent variables and $r$ is the independent variable. The ODEs are

$$69ru_m\frac{du_m}{dr}=8r^2\omega^2-16r\omega v-21{u_m}^2+48{v_m}^2-\frac{480\pi^2\nu r^3{u_m}^3}{Q^2}$$

$$21r(v_m-r\omega)\frac{du_m}{dr}+48ru\frac{dv_m}{dr}=-69u_m v_m+37r\omega u_m+\frac{480\pi^2\nu r^3(v_m-r\omega){u_m}^2}{Q^2}$$

To solve these equations I rearranged them in matrix form:

$$\textbf{A}=\begin{bmatrix}69\,r\,u_m & 0\\ 21(r\,v_m-r^2\omega) & 48\,r\,u_m \end{bmatrix}$$

$$\textbf{u}=\begin{bmatrix}u_m \\ v_m \end{bmatrix}$$

$$\textbf{f}=\begin{bmatrix} 8\,r^2\omega^2-16\,r\,\omega\,v_m-21{u_m}^2+48{v_m}^2-\frac{480\pi^2\nu}{Q^2}r^3{u_m}^3 \\ -69\,u_m\,v_m+37\,r\,\omega\,u_m-\frac{480\pi^2\nu}{Q^2}r^3{u_m}^2(v_m-r\,\omega) \end{bmatrix}$$

Then the equations are of the form

$$\textbf{A} \frac{d \textbf{u}}{dr}=\textbf{f}$$

I can rearrange this as

$$\frac{d \textbf{u}}{dr}=\textbf{A}^{-1}\textbf{f}$$

And solve it using a Runge-Kutta method. I have tried solving it using the following code in MATLAB:

function output = improvedPigford(flow,spin)
if nargin<2
    spin = 75;          % [RPM]
end
if nargin<1
    flow = 935;         % [mL/min]
end

w = spin*pi*2/60;       % rotation of wafer [rad/s]
d = 0.002;              %diameter of injection nozzle [m]

r0 = 0.5*d;             % begin solving here [m]
q = flow/(100^3*60);    % volumetric flow rate of liquid [m^3/s]
kv = 8.926780E-07;      % kinematic viscosity of water @ 25C [m^2/s]
u0 = q/(pi*r0^2);       % initial r-velocity (assume same as nozzle velocity) [m/s]
v0 = 0;                 % initial theta-velocity [m/s]

%% solve improved Pigford model
rRange = [r0 0.1];
initCond = [u0 v0];

options = odeset('Stats','on','RelTol',1e-6,'MStateDependence','strong');
dudr = @(r,u) improvedModel(r,u,q,kv,w);
sol = ode15s(dudr,rRange,initCond,options);
r = sol.x;
u = sol.y(1,:);
v = sol.y(2,:);
h = q./(2*pi*r.*u);

%% plot section
figure(1);clf;
plot(r*1000,u,'r-',r*1000,v,'b-');
xlabel('Radius (mm)');
ylabel('Fluid Velocity (m/s)');

%% output section
output=[r,u,v,h];

end

% --------------------------------------------------------------------------
function dudr = improvedModel(r,u,q,kv,w)
%{
Function is of the form A*u=f, so u=inv(A)*f
%}

A = [69*r*u(1), 0; 21*r*(u(2)-r*w), 48*r*u(1)];
f = [8*r^2*w^2 - 16*r*w*u(2) - 21*u(1)^2 + 48*u(2)^2 - (480*pi^2*kv*r^3*u(1)^3)/q^2;...
    -69*u(1)*u(2) + 37*r*w*u(1) + (480*pi^2*kv*r^3*u(1)^2*(u(2)-r*w))/q^2];
dudr = A\f;
end

However when I run this code the solution always blows up at about %r=0.0441%, no matter which RK function I use: ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb. The place where the solution blows up does change when I change the parameters like flow rate, spin speed, viscosity, etc., but it does always blow up.

Am I perhaps doing something wrong in my implementation of the Runge-Kutta method in MATLAB, or is it possible that the equations are simply ill-posed? Maybe there is simply a singularity at the location where it blows up or something? I've checked the integration of the N-S equations several times using symbolic mathematics software, and I'm fairly confident that I did it correctly. Otherwise though I'm kind of out of ideas of what to look at next. Any suggestions or ideas would be appreciated.

Update

Following the suggestion of @Geoff Oxberry I changed the MATLAB formulation to include the possibility that the mass matrix $\textbf{A}$ may be singular. Doing that the new code is this:

function output = improvedPigford(flow,spin)
if nargin<2
    spin = 200;         % [RPM]
end
if nargin<1
    flow = 1000;        % [mL/min]
end

w = spin*pi*2/60;       % rotation of wafer [rad/s]
d = 0.002;              %diameter of injection nozzle [m]

r0 = 0.5*d;             % begin solving here [m]
q = flow/(100^3*60);    % volumetric flow rate of liquid [m^3/s]
kv = 8.926780E-07;      % kinematic viscosity of water @ 25C [m^2/s]
u0 = q/(pi*r0^2);       % initial r-velocity (assume same as nozzle velocity) [m/s]
v0 = 0;                 % initial theta-velocity [m/s]

%% solve improved Pigford model
rRange = [r0 0.1];
initCond = [u0 v0];

dudr = @(r,u) improvedModel(r,u,q,kv,w);
M = @(r,u) massMatrix(r,u,w);
options = odeset('Stats','on','RelTol',1e-6,'Mass',M);
sol = ode15s(dudr,rRange,initCond,options);
r = sol.x;
u = sol.y(1,:);
v = sol.y(2,:);
h = q./(2*pi*r.*u);

%% plot section
figure(1);clf;
plot(r*1000,u,'r-',r*1000,v,'b-');
xlabel('Radius (mm)');
ylabel('Fluid Velocity (m/s)');

%% output section
output=[r,u,v,h];

end

function f = improvedModel(r,u,q,kv,w)
%{
Function is of the form M*u'=f
%}
f = [8*r^2*w^2 - 16*r*w*u(2) - 21*u(1)^2 + 48*u(2)^2 - (480*pi^2*kv*r^3*u(1)^3)/q^2;...
    -69*u(1)*u(2) + 37*r*w*u(1) + (480*pi^2*kv*r^3*u(1)^2*(u(2)-r*w))/q^2];
end

function M = massMatrix(r,u,w)
%{
Function calculates the mass matrix M for M*u'=f
%}
M = [69*r*u(1), 0; 21*r*(u(2)-r*w), 48*r*u(1)];
end

However, the result is the same: the solution still blows up at the same location.

Update 2

As @Kirill asked, I am showing a plot of what the solution looks like when it blows up.

$v_m$ goes to $-\infty$, while $u_m$ goes to $+\infty$. What makes me suspicious that something may be wrong is the fact that due to boundary conditions on $v(r,z)$ [$v=r\,\omega$ @ $z=0$ and $\frac{\partial v}{\partial z}=0$ @ $z=h(r)$], $v_m$ should almost certainly be positive. @Kirill found a sign discrepancy between my derivation and his derivation done using Mathematica, so I will be checking my derivation again, then getting back with my results.

Update 3

I went over the calculations again, and I think I did have a sign error as @Kirill postulated. After working through derivation one more time, I ended up with the following equations:

$$69r\,u_m\frac{du_m}{dr}=8r^2\omega^2-16r\omega v_m-21{u_m}^2+48{v_m}^2-\frac{480\pi^2\nu r^3{u_m}^3}{Q^2}$$

$$21r(v_m-r\omega)\frac{du_m}{dr}+48r\,u_m\frac{dv_m}{dr}=-69u_m v_m+37r\omega u_m-\frac{480\pi^2\nu r^3(v_m-r\omega){u_m}^2}{Q^2}$$

This is equivalent to the equations that @Kirill obtained doing the derivation in Mathematica. When solving this system of equations (using the above code but with the corrected term) I get the following plots for $u_m$ and $v_m$:

This appears to be identical to the plot that @Kirill showed.

For a further comparison, this derivation is an improvement of the Pigford model, which can be found here:

R.M. Wood and B.E. Watts, "The Flow, Heat, and Mass Transfer Characteristics of Liquid Films on Rotating Disks", Trans. Instn Chem. Engrs., Vol 51, 1973.

The derivation of the Pigford model is identical to mine, except that the final integration $\int_0^h(r) dz$ is simply skipped over, with the $u(r,z)$ and $v(r,z)$ simply being replaced by $u_m(r)$ and $v_m(r)$, with an $O(1)$ correction factor added to each equation. From this the Pigford model equations are:

$$r\,u_m \frac{du_m}{dr}={v_m}^2-\frac{12 K_1 \nu\, \pi^2\,r^3\,{u_m}^3}{Q^2}$$

$$r\,u_m \frac{dv_m}{dr}=-u_m v_m-\frac{12 K_2 \nu\, \pi^2\,r^3\,(v_m-r\,\omega){u_m}^2}{Q^2}$$

If I solve those two equations simultaneously, I get the following solution:

Which is qualitatively close to the equations that I derived. For a final comparison, I'll look at the height $h(r)$ calculated by both of the models:

The shapes of the two curves are a little different, but mostly only differ by a $O(1)$ constant, as one might expect because of the addition of the fitting constants in the Pigford model.

Answer 1

The first thing I can think of with your derivation is the $\mathbf{A}$ term. If at any point that term fails to be invertible, then your derivation is no longer valid. You should see the the value of $\mathbf{A}$ is at your singularity, and if it is still invertible.

Another way to diagnose this problem would be to use a DAE solver on $\textbf{A} \frac{d \textbf{u}}{dr}=\textbf{f}$ directly, instead of inverting $\mathbf{A}$; you could use ode15i for this purpose. If the DAE solver works, then there's a good chance singularity of $\mathbf{A}$ is to blame; if not, then you'll have to look for other alternative solutions.

EDIT: Here are some other possible causes of the behavior you're seeing that you may want to investigate:

• Is there a numerical stability issue with the integration methods you're using? Based on what you've said, it seems a little unlikely; you've tried methods that are generally geared toward stiff problems (ode15s: NDF/BDF methods, ode23s: Rosenbrock methods, ode23tb: TR-BDF2) and nonstiff problems (ode45 & ode23: explicit Runge-Kutta). Neither seems to help. Sometimes, there are telltale instabilities that can be seen when trying to solve a stiff problem with explicit methods (e.g., sometimes growing oscillations precede the solution blowing up).
• Is there a scaling issue with units? One way to diagnose this sort of error would be to look at the conditioning of the linear systems being solved within the implicit ODE (or DAE) solvers (so, ode15s, ode23s, ode23tb, etc.). Poorly conditioned systems also tend to be stiff (even if nondimensionalizing them would yield a nonstiff system), which would explain why nonstiff integrators also blow up. Nondimensionalization usually improves the conditioning in these cases, so you may also try to solve a nondimensionalized version of your problem instead. I'd check scaling next.
• Is your right-hand side regular enough over the interval of integration? By that, I mean, is it $k$ times continuously differentiable, where $k$ is greater than or equal to the order of time discretizations you're using? I also doubt that the problems you're seeing are due to derivatives failing to exist, because your right-hand side and mass matrix consist of polynomials.
• Is there something about the system analytically that causes it to blow up? I sort of doubt it, but exploring solved test problems in the literature and asymptotics can be helpful here. If there's a test problem in the literature that maps closely to what you're modeling here, you should use it as a benchmark for validation, because it can be helpful in debugging. You could find out that there's a "blow up time" beyond which something about your approximation breaks down.

Answer 2

Bearing in mind that I have just a vague idea of what results you are expecting, I want to point out that I do not get the same ODE system as you do. I'm using Mathematica.

Using your equations for $h(r)$, $u(r,z)$, $v(r,z)$ and the two NS equations, the following code generates the ODE system:

Module[{uu, vv, h},
 h = q/(2 Pi r u[r]);
 uu = 3 u[r] (z/h - (z/h)^2/2);
 vv = r w + 3/2 (r w - v[r]) ((z/h)^2 - 2 z/h);
 Integrate[
  {uu D[uu, r] - vv^2/r - kv D[uu, z, z],
   uu D[vv, r] + vv uu/r - kv D[vv, z, z]}
  , {z, 0, h}]
 ]
Thread[{0, 0} == (Numerator[Together[#/q]] & /@ %)]

$$ 0=480 \pi ^2 \nu r^3 u(r)^3-8 q^2 r^2 w^2+69 q^2 r u(r) u'(r)+21 q^2 u(r)^2+16 q^2 r w v(r)-48 q^2 v(r)^2\\ 0=-480 \pi ^2 \nu r^4 w u(r)^2+480 \pi ^2 \nu r^3 u(r)^2 v(r)-21 q^2 r^2 w u'(r)+21 q^2 r v(r) u'(r)+48 q^2 r u(r) v'(r)+69 q^2 u(r) v(r)-37 q^2 r w u(r) $$

The only difference seems to be the sign on the "480"-term in the second equation, I get it as "$-$".

Now, if I solve the system of ODEs that you have in the question using Mathematica's NDSolve, I also get a blow-up. But if I solve the system I generated above from the same equations, I get the following plot:

This looks okay to me, there is no blow-up, although I can't be sure this is right. So having said all this, are you absolutely sure your ODE system is actually correct?

Here is the ode solver invocation (the commented out sections (*...*) are two versions of your ODE system):

Module[{
  spin = 75,
  flow = 935,
  w, d, r0, rRange, q, kv, u0, v0, r, u, v
  , slv
  },
 w = spin Pi*2/60;
 d = 0.002;
 r0 = 0.5 d;
 q = flow/(100^3*60);
 kv = 8.926780*^-7;
 u0 = q/(Pi r0^2);
 v0 = 0;
 rRange = {r0, 0.1};
 (*slv=NDSolve[
 Join[{u[r0]==u0,v[r0]==v0},
 Thread[{u'[r],v'[r]}== Inverse[{{69r u[r],0},{21r (v[r]-r w),48r \
u[r]}}]
 .{8*r^2*w^2-16*r*w*v[r]-21*u[r]^2+48*v[r]^2-(480*Pi^2*kv*r^3*u[r]^3)/\
q^2,-69*u[r]*v[r]+37*r*w*u[r]+(480*Pi^2*kv*r^3*u[r]^2*(v[r]-r*w))/q^2}\
]]
 ,{u,v},{r,r0,rRange[[2]]}];*)
 (*slv=NDSolve[{
 69r u[r]u'[r]==8r^2w^2-16r w v[r]-21u[r]^2+48 v[r]^2-480Pi^2 kv \
r^3u[r]^3/q^2,
 21r (v[r]-r w)u'[r]+48 r u[r] v'[r]==-69 u[r]v[r]+37r w u[r]+480Pi^2 \
kv r^3(v[r]-r w)u[r]^2/q^2,
 u[r0]==u0,
 v[r0]==v0
 },{u,v},{r,r0,rRange[[2]]}];*)
 slv = NDSolve[{
    0 == -8 q^2 r^2 w^2 + 21 q^2 u[r]^2 + 480 kv \[Pi]^2 r^3 u[r]^3 + 
      16 q^2 r w v[r] - 48 q^2 v[r]^2 + 69 q^2 r u[r] 
\!\(\*SuperscriptBox["u", "\[Prime]",
MultilineFunction->None]\)[r], 
    0 == -37 q^2 r w u[r] - 480 kv \[Pi]^2 r^4 w u[r]^2 + 
      69 q^2 u[r] v[r] + 480 kv \[Pi]^2 r^3 u[r]^2 v[r] - 21 q^2 r^2 w 
\!\(\*SuperscriptBox["u", "\[Prime]",
MultilineFunction->None]\)[r] + 21 q^2 r v[r] 
\!\(\*SuperscriptBox["u", "\[Prime]",
MultilineFunction->None]\)[r] + 48 q^2 r u[r] 
\!\(\*SuperscriptBox["v", "\[Prime]",
MultilineFunction->None]\)[r],
    u[r0] == u0,
    v[r0] == v0
    }, {u, v}, {r, r0, rRange[[2]]}];
 If[Head[slv] =!= NDSolve, slv = slv[[1]], Abort[]];
 Print@Plot[Evaluate[{u[r], v[r]} /. slv], {r, r0, rRange[[2]]}]
 ]