# Can I use Runge-Kutta to solve these equations?

Edit: I'm going to give some more background and derivation to show how I got to these equations. I am basically following the derivation that is found in the appendix of the following paper:

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 is for the laminar fluid flow on the top of a rotating disk, where the liquid is dispensed on the center only.

Parameters: Fluid flow in the $r$-direction is $u(r,z)$, flow in the $\theta$-direction is $v(r,z)$. $\omega$ is the rotation rate of the disk, $h(r)$ is the height of the fluid film. Liquid flow rate is $Q=\int_0^{2\pi}\int_0^{h(r)}u(r,z)\,dz\,r\,d\theta=2\pi\,r\,h\,u_m$, where $u_m$ is the mean velocity in the $r$-direction, or $u_m(r)=\frac{1}{h(r)}\int_0^{h(r)}u(r,z)\,dz$. Similarly $v_m(r)=\frac{1}{h(r)}\int_0^{h(r)}v(r,z)\,dz$.

Assumptions: There is rotational symmetry such that $\frac{\partial}{\partial \theta}=0$. Flow is only horizontal so that flow in the $z$-direction $w=0$. Similarly $\frac{\partial P}{\partial z}=0$.

Assume that flow in the $r$-direction is parabolic and that $u=0$ at $z=0$ and $\frac{\partial u}{\partial z}=0$ at $z=h$. This gives us $u$ in terms of the local mean $u_m$ and fluid height $h$:

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

We do a similar assumption for flow in the $\theta$-direction, so that in this case $v=\omega\,r$ at $z=0$ and $\frac{\partial v}{\partial z}=0$ at $z=h$. This gives us:

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

## Derivation

At this point we take the steady-state Navier-Stokes equations in cylindrical coordinates, assuming that $\frac{\partial P}{\partial r}=0$, keeping the advection terms for $u$ and $v$, and keeping the $\frac{\partial^2}{\partial z^2}$ viscous terms:

$$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)$$

At this point you integrate both equations from $0$ to $h$ in the $z$-direction. Using the above assumptions for $u(r,z)$ and $v(r,z)$, you can perform the integration for the viscous term on the RHS of the N-S equations. For the advection terms on the LHS, Wood and Watts then simply replace $u$ with $u_m$ and $v$ with $v_m$ and throw in an $O(1)$ constant for each, eliminating the need to solve the integral for the advection terms. Finally you use the definition of $Q=2\pi\,r\,h\,u_m$ to eliminate $h$ from the equations with $h=Q/(2\pi\,r\,u_m)$. For Wood and Watts' derivation this gives the following ODE's:

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

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

My goal however, is to eliminate the need for the $K_1$ and $K_2$ fudge factor constants by using the full assumed forms of $u(r,z)$ and $v(r,z)$ and integrating all the terms on the LHS of the NS equations. Finishing the integration and simplification (done using MATLAB symbolic math package) I have the following two equations:

$$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}$$ $$48ru_m\frac{dv_m}{dr}=-21r(v_m-r\omega)\frac{du_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}$$

## Problem

I want to solve them simultaneously for $u(r)$ and $v(r)$, and initial conditions of $u(r_0)$ and $v(r_0)$ are known.

How to use the Runge-Kutta method is not my question. I've used it in the past and know how it works. My question/problem comes from the $\frac{du}{dr}$ term in the 2nd equation. Without that term I could just solve the two equations simultaneously using Runge-Kutta quite easily. But with it I'm having trouble understanding the best way to solve it. I could simply solve the first equation for $\frac{du}{dr}$ and substitute it into the second equation, but that would make the formulation really messy and I'd like to avoid that. Is there a way I could introduce a simple 3rd equation to deal with that? I'm thinking of something similar to how when you solve a 2nd degree ODE like $y''=f(x)$ using the Rung-Kutta method you split it into two equations:

$$y_1'=y_2$$

$$y_2'=f(t)$$

However I'm having trouble coming up with a formulation that would let me do that.

## Update

With some help from the folks over at Computational Science Stackexchange, I have solved the system of equations. The solution can be seen here.

• You didn't finish your algebraic derivation yet. Like you said, you can solve the first equation for du/dr and then substitute in the second. Having said that, the coefficients in your equations look very funny. How did you end up with numbers like 69 and 37? Are you sure you didn't make a mistake somewhere along the way? – CuriousOne Dec 22 '14 at 22:19
• Even if you didn't want to do that substitution (and assuming everything else is correct) -- actually write out the discrete form of your equations as they are and it should be pretty easy to formulate the integration. Particularly if you write the left-hand side as a matrix multiplying a vector of {u,v}. – tpg2114 Dec 22 '14 at 22:38
• can you define what $\omega$ is - is it a constant - I would normally guess it is $d\theta /dt$, but here that seems to be $v$ - would help to answer if that was clear – tom Dec 22 '14 at 22:52
• Crossposted to math.stackexchange.com/q/1078136/11127 Would Computational Science be a better home for this question? – Qmechanic Dec 22 '14 at 22:53
• Is $\omega$ the rate of spin of the disk? – tom Dec 22 '14 at 23:12

The way I see it you have two choices.

Either you introduce a third variable for $du/dr$ - let us call it $y$ and calculate $dy/dr$ using your first equation. The you should be able to do the Runge Kutta with three variables. and you can set $du/dr = y$

Or you substitute in for $du/dr$ from the first equation into the second.

I think the second method, though messy, is probably the best because it looks as if $y$ is not an indepedent variable as it is completely defined by $u$, $v$, $r$ etc.

• If I did that though, wouldn't I essentially be calculating $d^2u/dr^2$? Since my equations are 1st order I would then have to have some initial condition for $du/dr$, which I don't have. – Derek Dec 23 '14 at 16:33
• @Derek yes - that is why I think option 2 is better (substitute in for $du/dr$ in 2nd equation using information from first equation) - the reason I suggest substituting it with $y$ is that it makes it easier to handle the $du/dr$ term on the right hand side of the second equation. + you can calculate the initial value for $du/dr$ with your first equation. – tom Dec 23 '14 at 16:52
• I could do that with a simple Euler method, but I don't think that would work with an automated Runge-Kutta code, since it needs to evaluate the RHS of the equations at several points for each iteration. I would prefer to be able to just use the Runge-Kutta functions in MATLAB to solve this, but if I can't find another way I can always revert to the Euler method. – Derek Dec 23 '14 at 17:08
• @Derek, if it will work with Euler - it can also be programmed with Runge Kutta, but it is just a bit more difficult to programme in my experience. I'm afraid I am not familiar with the functions in MATLAB. – tom Dec 24 '14 at 23:11
• I think this explanation was good enough to get me to the next step - actually getting it to work. Thanks for the help. – Derek Jan 7 '15 at 22:03