# Momentum of slowly spinning (viscous) fluid

If we have a massless cylindrical container (or radius $R$) with a liquid of certain density $\rho$ and viscosity $\mu$ at rest. Then at time zero we impart a constant rotational velocity $\Omega$ on the cylinder and watch as the liquid accelerates from the outer walls inwards (due to the viscosity).

I want to know what is the function of total angular momentum $L$ with time (and consequently the effective mass moment of inertia of the liquid $L=I_{eff} \Omega$). More specifically the radius of gyration $I_{eff} = m \kappa^2$ where $m=\rho \pi R^2 h$.

I looked at concentric cylindrical slices in order to derive the equations of motion but I am stumbling at the shape of the tangential velocity as a function of radius $r$ and time $t$.

Based on of $\nabla$ in cylindrical coordinates (for a Newtonian fluid) I think shear stress is $$\tau_{r\theta} = \mu \left( \frac{\partial v_\theta}{\partial r} - \frac{v_\theta}{r} \right)$$

The cylindrical slice has surface area $A = 2 \pi r h$, or ${\rm d}A = 2 \pi h {\rm d}r$.

The volume is ${\rm d}V = A {\rm d}r = 2 \pi r h {\rm d}r$.

I think the radial force balance acting on the volume due to the shear stress is $$A {\rm d}\tau_{r\theta} + \tau_{r\theta} {\rm d}A = \dot{v_\theta} \rho {\rm d} V$$

In which this leads to a differential equation $$\tau_{r\theta}'= \frac{\partial \tau_{r\theta}}{\partial r} = \rho \dot{v_\theta} - \frac{\tau_{r\theta}}{r}$$

The shear stress slope is (from chain rule) $$\tau_{r\theta}' = \mu \left( v''_{\theta} - \frac{v'_\theta}{r}+ \frac{v_\theta}{r^2} \right)$$

With some algebra I get finally that

$$\dot{v}_{ \theta} = \frac{\mu}{\rho} v''_{\theta}$$

So the acceleration of the cylindrical slice is proportional to the curvature of the velocity profile.

Here is where I am stuck. I am not sure how to proceed to derive $v_\theta (r,t)$.

• What I am really interested in is the shape of the momentum curve over time. – John Alexiou Jan 14 '16 at 2:27
• Well, that should be basically the solution of the last PDE, shouldn't it? – Victor Pira Jan 14 '16 at 7:42
• I can do this for steady-state, but I haven't seen/done the derivations for unsteady. Are you expecting an analytical solution or is something that requires numerical integration acceptable? – tpg2114 Jan 14 '16 at 8:33
• First I'd like to know if there is an analytical solution. If that is not possible, then any numerical results will give me the typical velocity profile and the momentum over time. If it exponential, I'd like to know the coefficient $\beta$ from $\exp(-\beta t)$ terms. – John Alexiou Jan 14 '16 at 13:17
• Hrm... looks like my answer at least got the ball rolling, even if it wasn't a good one! – tpg2114 Jan 20 '16 at 18:29

Assuming the equation you found is correct, it's just the 1-dimension heat equation: $$D^2 \frac{\partial^2 v}{\partial r^2} = \frac{\partial v}{\partial t},$$ where $D^2 \equiv \mu/\rho$. We want to solve it on the domain $r \in [0, R]$, $t \in [0, \infty)$, subject to the boundary conditions $v(r,0) = 0$, $v(R, t) = R \Omega$, and $v(0,t) = 0$. (This last condition is inserted so that the velocity field remains continuous at the origin; remember that $v$ is only the tangential component of the field.)

The steady-state solution $v_\infty(r,t)$ to this equation is pretty obvious: if $\dot{v}_\infty = 0$, then $v''_\infty = 0$ as well, and so $$v_\infty(r,t) = \Omega r.$$ This makes sense: the whole cylinder is rotating rigidly at late times.

Now define $\delta v(r,t) \equiv v(r,t) - v_\infty(r,t)$. By construction, $\delta v$ also satisfies the heat equation, but with different boundary conditions: $$\delta v(R,t) = \delta v (0, t) = 0; \qquad \delta v (r,0) = - \Omega r \equiv \delta v_0(r).$$ This is just a standard heat-diffusion problem with an initial heat distribution $\delta v_0(r)$. To solve this, we use separation of variables; it's not too hard to show that any solution of the form $$f_n(r,t) = \sin \left( \frac{n \pi r}{R} \right) e^{-D^2 \pi^2 n^2 t/R^2}$$ will work. Assuming that $\delta v(r,t) = \sum_n A_n f_n(r,t)$, we have $$\delta v_0 (r) = - \Omega r = \sum_n A_n \sin \left( \frac{n \pi r}{R} \right),$$ i.e., $\delta v_0$ is expressed as a Fourier series. Working through this (I used Mathematica to save time), we get $$A_n = \frac{2 \Omega R (-1)^n}{\pi n}$$

Thus, the full solution for the tangential velocity as a function of time and space is: $$\boxed{ v(r,t) = \Omega r + 2 \Omega R \sum_{n=1}^\infty \frac{(-1)^n}{\pi n} \sin \left( \frac{n \pi r}{R} \right) e^{-D^2 \pi^2 n^2 t/R^2}.}$$ It may be possible to sum this up into a closed-form expression; I'll let you know if I make any progress on this.

• Any series of the form $\Omega r + 2 \Omega R \sum \left( \exp(-\beta n^2 t) K \sin(n \kappa r) \right)$ solves the equation when $\beta = \frac{\mu}{\rho} \kappa^2$. Why did you choose $\kappa = \frac{\pi}{R}$ for the shape function. Intuitively it makes sense to fit the BC, but I wonder if there is any other reason. What about $\kappa = \frac{\pi}{2 R}$ such that $\sin()=1$ when $r=R$? – John Alexiou Jan 20 '16 at 19:01
• Is there an elaborate reason why you use the cartesian 1D heat equation? For example, that the advection terms in Navier Stokes turn the cylindrical laplacian into the cartesian laplacian? Or was it just a mistake in setting up the problem? – ignacio Jan 21 '16 at 17:37
• @ignacio: I just assumed that the equation in the problem statement was correct and proceeded from there. A similar technique could be used if the radial equation was different, and you'd end up with Bessel functions instead of sinusoids. See Thomas's answer for the result; in fact, I expect that his answer is the correct answer to the original problem (while mine is a correct answer to a PDE that's irrelevant to the original problem.) – Michael Seifert Jan 21 '16 at 18:04

First of all, I think the equation of motion is not correct. I believe the equation of motion is $$\dot{v}_\theta = \nu \left( v^{\prime\prime}_\theta + \frac{v^\prime_\theta}{r} - \frac{v_\theta}{r^2} \right) ,$$ where $\nu=\mu/\rho$. This is easiest to derive using the result for the Laplacian of a vector field in cylindrical coordinates.

In order to derive this result from force balance we have to consider a volume element $V=L_z L_\theta L_r$ in cylindrical coordinates. Here, $L_\theta=r\Delta\theta$ and $L_r=\Delta r$. We can ignore $L_z=\Delta z$ as nothing depends on $z$. The force on $L_\theta$ is $F=\tau_{r\theta}\hat{e}_\theta L_\theta$. The net acceleration is due to the difference between the forces on the faces at $r$ and $r+\Delta r$. We get $\dot{v}_\theta=\nu[\tau_{r\theta}'+\tau_{r\theta}/r]$ as explained in the question. The force on $L_r$ is $F=\tau_{r\theta}\hat{e}_r$, and the net acceleration is due to the difference between the forces on $L_r(\theta)$ and $L_r(\theta+\Delta\theta)$. Using $\partial \hat{e}_r/(\partial \theta)=\hat{e}_\theta$ we get $\dot{v}_\theta=\nu\tau_{r\theta}/r$. Combining both we get $$\dot{v}_\theta=\nu\left(\tau_{r\theta}'+2\frac{\tau_{r\theta}}{r}\right),$$ which agrees with the formula above.

The basic time scale is given by vorticity diffusion, so we expect $v_\theta \sim \exp(-c\nu t/R^2)$, but a little more effort is required to get $c$ and the precise $r$ dependence. You can tackle this by Fourier expanions, but because of the cylindrical geometry you should really use a Bessel expansion. Make a separation ansatz $v_\theta(r,t)=g(r)f(t)$. The equation of $g(r)$ is the Bessel $J_1$ differential equation. Then we expand the boundary condition in $J_1(\lambda_n r/R)$, using orthogonality of the Bessel function with respect to the zeros $\lambda_n$. Then (for spinning down) $$v_\theta = 2\Omega R \sum_n \frac{J_1(\lambda_n r/R)}{\lambda_nJ_0(\lambda_n)} \exp\left(-\lambda_n\nu t/R^2\right).$$ The solution for spinning up is just $\exp()\to 1-\exp()$.

Bonus: To convince myself that the solution does indeed satisfy the boundary conditions, I plotted the answer (oscillations near the boundary are due to slow convergence of the Bessel expansion). The figure shows $v_\theta(r)$ for different $t$.

• Shouldn't the solution approach $v_\theta \propto r$ as $t \to \infty$, i.e., the whole cylinder is rotating rigidly? It's not at all obvious to me that your solution does this. – Michael Seifert Jan 20 '16 at 17:30
• Unless I made a mistake the sum should be the Bessel expansion of $r/R$. – Thomas Jan 20 '16 at 18:01
• Makes sense. I did basically the same approach you did, but using the PDE described in the original question (which leads to a straight-up Fourier series rather than a Bessel series.) But I haven't worked with viscous flows enough to be able to say whether the original equation is actually correct. – Michael Seifert Jan 20 '16 at 18:27
• The viscous drag is $\nu\nabla^2 v$. In cylindrical coordinates this must give Bessel functions. – Thomas Jan 20 '16 at 18:31
• I did not arbitrarily drop terms. I plugged $\tau$ and $\tau'$ into the equation of motion $A \tau' {\rm d}r + \tau {\rm d}A = \dot{v} A {\rm d}r$ to get $$\left. r \tau' + \tau = \rho \dot{v} r \right\} \mu r v'' = r \rho \dot{v}$$ – John Alexiou Jan 20 '16 at 18:50

Edit: This use of non-dimensionalization is wrong because it is incompatible with the boundary conditions. More information here.

We're looking for a solution to the fluid velocity $\vec u$ that looks like $\vec u = u(r,t)\hat\theta$. Given the symmetries of the problem, we have to solve Navier Stokes: $$\partial_t u=\nu(\nabla^2 u - \frac u{r^2})$$ The only quantities that appear in this equation, apart from $u$, are $t$, $r$ and $\nu$. Because $u$ has to be a function of non-dimensional variables, we look for non-dimensional combinations of those quantities and find only one: $\frac{r^2}{\nu t}$. We'll call this quantity $z$ and simplify our expression for $\vec u$: $$\vec u = u(z)\hat\theta$$ Now let's look at the derivatives. Because $u(z)$ can only depend on $r$ and $t$ through $z$, we have $$\partial_t u = \partial_t z u' = -\frac z t u'\\ \partial_r u = \partial_r z u' = 2\frac z r u'$$ This leads to the following expression for the laplacian: $$\nabla^2 = \frac 1 r \partial_r(r\partial_r)\\ = \frac{2z}{r^2}\frac d{dz}\left(2z\frac d{dz}\right)\\ = \frac{4z}{r^2}\left(\frac d{dz}+z\frac{d^2}{dz^2}\right)$$ Inserting this into Navier-Stokes, we get $$-\frac z t u' = \frac{4z\nu}{r^2}(u'+zu'')-\nu\frac u{r^2}$$ We multiply everything by $\frac{r^2}{4z^2\nu}$: $$-\frac 1 4 u' = u''+\frac{u'}z-\frac u{4z^2}$$ Rearranging, $$u''+\left(\frac 1 4+\frac 1 z\right)u'-\frac 1{4z^2} u = 0$$ That's all I have so far, but you can see that through non-dimensionalization you can turn the original PDE into a (simpler?) ODE.

I didn't forget about this one! I cannot find fault with your derivation and the resulting PDE. I may have missed something also though, so who knows. But what I can do is come up with a solution for the PDE you do have, which should help if it turns out the original isn't right for one reason or another.

Starting from (where $\nu = \mu/\rho$):

$$\frac{\partial v_\theta}{\partial t} = \nu \frac{\partial^2 v_\theta}{\partial r^2}$$

and taking the Fourier transform from $r\rightarrow k$, we get:

$$\frac{\partial \hat{v}_{\theta}}{\partial t} = \nu (i k)^2 \hat{v}_{\theta}$$

which has an exponential function as a solution. The solution is thus:

$$\hat{v}(k,t) = e^{-\nu k^2 t} \hat{v}(k,0)$$

where $\hat{v}(k,0)$ is the Fourier transform of your initial conditions (the impulsively started wall at $r = R$). For those initial conditions, assuming the wall is started with velocity $U$, the transform of the initial conditions is:

$$\hat{v}_\theta (k,0) = \frac{1}{\sqrt{2\pi}} e^{i k R} U$$

We can now take the inverse Fourier transform of this result and we will get:

$$v_\theta (r,t) = \frac{U}{2\sqrt{\pi \nu t}} e^{-(r-R)^2/4\nu t}$$

Or at least that's what the math says. The part that bothers me now are the units; the RHS seems to have units of $\text{s}^{-1}$ unless there are some hidden units in $\pi$, which given the weirdness of Fourier transforms wouldn't surprise me at all.

• With a cursory look it does not seem this meets the natural boundary conditions of $v_\theta(r=0,t)=0$ and $v_\theta(r=R,t)=\Omega R$. Also the initial conditions are $v_\theta(r<R,0)=0$ I think. – John Alexiou Jan 20 '16 at 14:59
• @ja72 Yeah, I'm still trying to sort that out. I think the initial Fourier transform isn't correct because it assumes an infinite domain. It's possible to do it on a truncated domain, perhaps by using some Heaviside functions, that may enforce the boundary conditions better. – tpg2114 Jan 20 '16 at 15:02
• I agree, it is very tricky to marry the boundary conditions to the solution. – John Alexiou Jan 20 '16 at 15:04
• @ja72 I feel like this approach (Fourier -> solve -> Inverse Fourier) will work to get an analytical solution; there's just some subtleties in the transformation. This isn't an area I am terribly well versed in, but I've been using it recently to analyze the Euler equations. – tpg2114 Jan 20 '16 at 15:07
• If look at with separation of variables I find a solution of the $v(r,t) = \exp(-\nu \kappa^2 t) \left( V_1 \sin(\kappa r) + V_2 \cos(\kappa r)\right)$ form. – John Alexiou Jan 20 '16 at 15:14

Question;

I am not sure how to proceed to derive $vθ(r,t)$

It may be easier to make a clean simple start rather than trying to seek a mistake.

If the fluid is rotating slowly, and it's viscous. You basically have a laminar flow conditions, which makes the whole really simple as the velocity distribution is an exponential curve. Ie parabel. At your case the whole flow will finally rotate as a rigid object, but up to there the velocity change is a cut from this parabel curve.

Here's a helping video about "Rotating Flows".. The needed stuff is at 1-4 min. Basically this picture; $U$ is the velocity. The stuff at 3rd row is the specialties in rotating flows; Centrifugal force and Coriolis force.

I Hope this short answer helps you to proceed.