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

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-');
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-');
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.

• Is the sign in front of "480" in the 2nd equation correct? (See my answer.) – Kirill Jan 7 '15 at 20:12
• Sorry, another thing: please, please, please include a plot in your question. There are so many different ways a solution can blow up! – Kirill Jan 7 '15 at 20:26
• OK, I'll put in a plot showing how it's blowing up. – Derek Jan 7 '15 at 20:50
• @Kirill, I put up a plot of the blown up results. I'm still working on checking over the derivation to see if there is a sign error. – Derek Jan 7 '15 at 21:17
• @DKS sure thing. Here it is: archive.icheme.org/cgi-bin/… – Derek Jun 22 '17 at 8:19

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]]}]
]

• Thanks for taking a look at this. For the analytical derivation Mathematica is certainly much easier to work with than MATLAB. I'll check over my analysis again and get back to you. – Derek Jan 7 '15 at 20:52
• @Derek It's no problem, I actually like this question. Even if my answer is wrong, scientific computing is full of sign errors. – Kirill Jan 7 '15 at 21:10
• It looks like it was a sign error after all. Thanks for the help! – Derek Jan 8 '15 at 1:48

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.
• I can try that. DAE solvers on MATLAB use the symbolic math toolbox, correct? I have that, but mostly only use it for MuPAD... So try something like these examples? – Derek Jan 6 '15 at 18:53
• I don't believe the DAE solvers use the MATLAB symbolic toolbox. – Geoff Oxberry Jan 6 '15 at 19:25
• I think I found it. The help documentation for odeset has a section Mass Matrix and DAE Properties that explains what to do when the mass matrix may be singular. I'll try it out and get back to you. – Derek Jan 6 '15 at 19:47
• I reformulated it as a DAE with a possibly singular mass matrix, but it didn't affect the answer. I put my new code as an update on my original post. – Derek Jan 6 '15 at 20:06