Imagine you have a fluid inside a cylinder that suddenly starts rotating at angular velocity $\omega$:
The fluid is initially stationary but eventually settles to a flow field that does not change further with time. What is the velocity of the fluid in the meantime? This question can be answered using Bessel functions.
Simplifying the Governing Equations
For simplicity, consider the fluid to be incompressible. This is a perfectly reasonable approximation for, say, water. The principle of conservation of mass then leads to a constraint on the velocity $\vec{u}$:
$$\displaystyle \nabla\cdot\vec{u} = 0 $$
Conservation of momentum introduces the density $\rho$, pressure $P$, and kinematic viscosity $\nu$:
$$\displaystyle \frac{\partial\vec{u}}{\partial t} + \vec{u}\cdot\nabla\vec{u} = -\frac{1}{\rho}\nabla P + \nu\Delta\vec{u} $$
Since the region of interest is circular, we will make some further (reasonable) assumptions that will simplify the $\nabla$ differential operators into partial derivatives on only some of the coordinate directions. First, any quantity $q$ will be rotationally symmetric:
$$\displaystyle \frac{\partial q}{\partial \theta} = 0 $$
Second, the component of velocity in the radial (i.e. outward from the center) direction is zero:
$$\displaystyle u_r = 0$$
Then the incompressibility constraint becomes:
$$\displaystyle \nabla\cdot\vec{u} = 0 = \frac{1}{r}\frac{\partial}{\partial r} (ru_r) + \frac{1}{r}\frac{\partial u_\theta}{\partial \theta} \Rightarrow 0 = 0 $$
since $u_r = 0$ everywhere and $\partial/\partial\theta = 0$ everywhere. This means that incompressibility is automatically satisfied so we only really have to worry about the momentum equation now. That equation becomes:
$$\displaystyle \begin{aligned} \frac{\partial\vec{u}}{\partial t} + u_r \frac{\partial\vec{u}}{\partial r} + u_\theta\frac{\partial\vec{u}}{\partial\theta} &= \nu\left( \frac{1}{r}\frac{\partial}{\partial r}\left( r\frac{\partial u_\theta}{\partial r} \right) – \frac{u_\theta}{r^2} \right)\hat{\theta} – \frac{1}{\rho}\frac{\partial P}{\partial r}\hat{r} – \frac{1}{\rho}\frac{\partial P}{\partial \theta}\hat{\theta} \\ \frac{\partial u_\theta}{\partial t}\hat{\theta} &= \nu \frac{\partial}{\partial r}\left( \frac{1}{r}\frac{\partial}{\partial r} (ru) \right)\hat{\theta} – \frac{1}{\rho}\frac{\partial P}{\partial r}\hat{r} \end{aligned}$$
Note that in cylindrical coordinates, the vector Laplacian $\nabla\vec{u}$ is not the scalar Laplacian applied to each component of $\vec{u}$.
Since the only $\hat{r}$ component is the radial pressure gradient on the right-hand side, it is zero so we can care only about the $\hat{\theta}$ component, leaving us with the equation we need to solve:
$$\displaystyle \frac{\partial u_\theta}{\partial t} = \nu\frac{\partial}{\partial r}\left( \frac{1}{r}\frac{\partial}{\partial r} (ru_\theta) \right) $$
Perhaps unsurprisingly, this is a diffusion equation in cylindrical coordinates. A standard problem in introductory fluid mechanics is to solve the steady form of this equation, i.e. with $\partial u_\theta/\partial t = 0$, which is easily seen to have the general solution:
$$\displaystyle (u_\theta)_{steady} = c_1 r + \frac{c_2}{r} $$
The $c_1$ term becomes unbounded as $r\to\infty$ making it unrealistic for solutions outside the cylinder, whereas the $c_2$ term becomes unbounded as $r\to 0$ making it unrealistic for solutions inside the cylinder which we’re interested in. Supposing our cylinder has a radius of $L$, the steady-state solution must therefore be:
$$\displaystyle (u_\theta)_{steady} = \omega r $$
It will simplify our derivation if we subtract off this steady-state solution from $u_\theta$ and solve for the difference $u$:
$$\displaystyle u = u_\theta – (u_\theta)_{steady} = u_\theta – \omega r $$
which satisfies the same PDE we derived just previously:
$$\displaystyle \frac{\partial u}{\partial t} = \nu\frac{\partial}{\partial r}\left( \frac{1}{r}\frac{\partial}{\partial r} (ru) \right) = \nu \left( \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} – \frac{u}{r^2} \right)$$
Solving the PDE
Proceed by separation of variables:
$$\displaystyle u(t,r) = T(t)R(r) $$
Substituting this expression into the PDE for $u$ and dividing by $TR$, we get:
$$\displaystyle T’ R = \nu \left( TR^{”} + \frac{1}{r} TR’ – \frac{1}{r^2} TR \right) $$
$$\displaystyle \frac{T’}{T} = \nu \left( \frac{R^{”}}{R} + \frac{1}{r} \frac{R’}{R} – \frac{1}{r^2} \right) $$
We are justified in equating both sides of the above to a constant because the left-hand side can depend only on $t$ and not $r$ whereas the right-hand side can depend only on $r$ and not $t$, which is only possible if both sides depend on neither $t$ or $r$ meaning they equal some constant. That constant must be negative since otherwise the equation for $T$ would have solutions that grow unbounded as $t$ increases. So we may write:
$$\displaystyle \frac{T’}{T} = \nu \left( \frac{R^{”}}{R} + \frac{1}{r} \frac{R’}{R} – \frac{1}{r^2} \right) = -\sigma_k^2$$
for some real numbers $\sigma_k$ (we will soon see why we can enumerate the possible constants by integers $k$). From this point we can solve for $T(t)$ and $R(r)$ separately. For $T$ we have:
$$\displaystyle T’ = -\sigma_k^2 T \Rightarrow T(t) = Ce^{-\sigma_k^2 t}$$
It is important that these solutions all decrease as $t\to\infty$, because that means the deviations $u$ from the steady-state solution will approach zero which is what we expect. Now for $R$:
$$\displaystyle r^2R^{”} + rR’ + \left( \frac{\sigma_k^2}{\nu}r^2 – 1 \right)R = 0 $$
We can recover the Bessel differential equation
$$\displaystyle x^2 J^{”}_p + x J’_p + (x^2 – p^2)J_p = 0 $$
by making the change of variable $x = \sigma_k r/\sqrt{\nu}$ and letting $p = 1$. So the solution for $R$ is:
$$\displaystyle R(r) = J_1(x) = J_1\left( \frac{\sigma_k r}{\sqrt{\nu}} \right) $$
Recombining the solutions for $T$ and $R$ we get the form of the solution to the original PDE for $u$:
$$\displaystyle u(t,r) = \sum_k C_k e^{-\sigma_k^2 t}J_1\left( \frac{\sigma_k r}{\sqrt{\nu}} \right) $$
It remains to decide the values of $C_k$ and $\sigma_k$. The $\sigma_k$ arise from the boundary conditions, which we have not yet specified. Recall that the unknown $u$ is the difference between the unsteady solution and the steady-state solution. The velocity at the cylinder edge in the unsteady solution must match the velocity of the cylinder itself, as it must also in the steady solution. Therefore the difference $u$ at $r=L$ must be zero:
$$\displaystyle u(t,r=L) = 0 = \sum_k C_k e^{-\sigma_k^2 t}J_1\left( \frac{\sigma_k L}{\sqrt{\nu}} \right) $$
Since the exponential terms are never zero for any $t$, this equation can only be satisfied if $J_1\left( \frac{\sigma_k L}{\sqrt{\nu}} \right) = 0$, meaning that $\frac{\sigma_k L}{\sqrt{\nu}}$ must coincide with a zero $\lambda_{1,k}$ of the Bessel function $J_1$:
$$\displaystyle \frac{\sigma_k L}{\sqrt{\nu}} = \lambda_{1,k} \Rightarrow \sigma_k = \lambda_{1,k}\frac{\sqrt{\nu}}{L} $$
So we have for $u$:
$$\displaystyle u(t,r) = \sum_{k=1}^\infty C_k e^{-\lambda_{1,k}^2\nu t/L^2} J_1\left(\lambda_{1,k}\frac{r}{L}\right) $$
Now the $C_k$ are determined by considering the initial conditions. Since the fluid is initially at rest, $u_\theta = 0 \Rightarrow u = u_\theta – (u_\theta)_{steady} = -\omega r$. So we need $C_k$ such that:
$$\displaystyle u(t=0,r) = -\omega r = \sum_{k=1}^\infty C_k J_1\left(\lambda_{1,k}\frac{r}{L}\right) $$
We are aided by the orthogonality property of Bessel functions:
$$\displaystyle \int_0^L rJ_p\left(\lambda_{p,k}\frac{r}{L}\right) J_p\left(\lambda_{p,m}\frac{r}{L}\right) \,dr = \frac{1}{2} \delta_{k,m} L^2 J_{p+1}(\lambda_{p,k})^2 $$
which allows us to isolate individual $C_k$:
$$\displaystyle \begin{aligned} \int_0^L -\omega r\cdot rJ_1\left( \lambda_{1,m} \frac{r}{L} \right)\,dr &= \int_0^L \sum_{k=1}^\infty C_k J_1\left(\lambda_{1,k}\frac{r}{L}\right)\cdot rJ_1\left( \lambda_{1,m} \frac{r}{L} \right)\,dr \\ &= C_m \frac{1}{2}L^2 J_2(\lambda_{1,m})^2 \end{aligned}$$
The Bessel functions have the property that
$$\displaystyle x^p J_{p-1} = \frac{d}{dx} (x^p J_p) $$
which allows us to write the integral on the left-hand side as:
$$\displaystyle \begin{aligned} \int_0^L r^2 J_1\left( \lambda_{1,m} \frac{r}{L} \right)\,dr &= \int_0^{\lambda_{1,m}} \left( \frac{sL}{\lambda_{1,m}} \right)^2 J_1(s) \frac{L}{\lambda_{1,m}} \,ds\\ &= \left( \frac{L}{\lambda_{1,m}} \right)^3\int_0^{\lambda_{1,m}} s^2 J_1(s) \,ds \\ &= \left( \frac{L}{\lambda_{1,m}} \right)^3 s^2 J_2(s) \Bigg|^{s=\lambda_{1,m}}_{s=0} \\ &= \left( \frac{L}{\lambda{1,m}} \right)^3 \lambda_{1,m}^2 J_2(\lambda_{1,m}) \end{aligned}$$
making the change of variable $s = \lambda_{1,m}r/L$. So for $C_m$ we end up with:
$$\displaystyle C_m = \frac{-2L\omega}{\lambda_{1,m} J_2(\lambda_{1,m})} $$
for a final result of:
$$\displaystyle u_\theta(t,r) = \omega r – 2L\omega\sum_{m=1}^\infty \frac{1}{\lambda_{1,m} J_2(\lambda_{1,m})} e^{-\lambda_{1,m}^2\nu t/L^2} J_1\left(\lambda_{1,m}\frac{r}{L}\right) $$