# Peeter Joot's (OLD) Blog.

• ## Archives

 Adam C Scott on avoiding gdb signal noise… Ken on Scotiabank iTrade RESP …… Alan Ball on Oops. Fixing a drill hole in P… Peeter Joot's B… on Stokes theorem in Geometric… Exploring Stokes The… on Stokes theorem in Geometric…

• 287,565

## Steady state velocity profile of stirred cup of non-bottomless coffee.

Posted by peeterjoot on April 29, 2012

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Motivation.

Having tackled the problem of the spin down of a bottomless cup of coffee, lets try setting up the harder problem with a non-infinite cup. It turns out that this is a lot harder. Even the steady state solution requires a superposition of Bessel functions (whereas we only had that in the bottomless problem when we tried to find the time evolution after ceasing the stirring).

I can find the form of the solution, but don’t know how to actually apply the boundary value constraints. I also find the no-slip constraints themselves become problematic, because they lead to an inconsistency at the point of contact of a moving and a static interface. Before continuing, I need to find out how to deal with that inconsistency.

# Navier-Stokes for the problem.

Working in cylindrical coordinates is the only sensible option. Let’s look first at the steady state problem, with the cup being stirred at at constant rate with the unrealistic rotating cylinder stir stick used previously. We’ll assume an azimuthal velocity profile

\begin{aligned}\mathbf{u} = u(r, z, t) \hat{\boldsymbol{\phi}}.\end{aligned} \hspace{\stretch{1}}(2.1)

Let’s first verify that this satisfies our incompressible condition

\begin{aligned}\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u}&=\left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi + \hat{\mathbf{z}} \partial_z \right) \cdot ( u \hat{\boldsymbol{\phi}} ) \\ &=\not{{\hat{\mathbf{r}} \cdot \hat{\boldsymbol{\phi}}}} \partial_r u + \frac{\hat{\boldsymbol{\phi}}}{r} \cdot (u \partial_\phi \hat{\boldsymbol{\phi}} + \hat{\boldsymbol{\phi}} \not{{\partial_\phi u}})+ \not{{\hat{\mathbf{z}} \cdot \hat{\boldsymbol{\phi}}}} \partial_z u.\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.2)

The only term that survives is the $\partial_\phi \hat{\boldsymbol{\phi}} = -\hat{\mathbf{r}}$ but that’s perpendicular to $\hat{\boldsymbol{\phi}}$ so we have $0 = \boldsymbol{\nabla} \cdot \mathbf{u}$ as desired. This leaves us with

\begin{aligned}\frac{u}{r} \partial_\phi ( u \hat{\boldsymbol{\phi}} )= -\frac{1}{{\rho}} \left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi + \hat{\mathbf{z}} \partial_z \right) p+ \nu \left( \frac{1}{{r}} \partial_r ( r \partial_r ) + \frac{1}{r^2} \partial_{\phi \phi} + \partial_{zz} \right) u \hat{\boldsymbol{\phi}} + \mathbf{g}\end{aligned} \hspace{\stretch{1}}(2.3)

Splitting into $\hat{\boldsymbol{\phi}}, \hat{\mathbf{r}}, \hat{\mathbf{z}}$ coordinates respectively we have

\begin{aligned}0 = -\frac{1}{{r \rho}} \partial_\phi p + \nu \left( \frac{1}{{r}} \partial_r ( r \partial_r u ) - \frac{u}{r^2} + \partial_{zz} u \right) \end{aligned} \hspace{\stretch{1}}(2.4a)

\begin{aligned}- \frac{u^2}{r} = -\frac{1}{{\rho}} \partial_r p \end{aligned} \hspace{\stretch{1}}(2.4b)

\begin{aligned}0 = -\frac{1}{{\rho}} \partial_z p - g.\end{aligned} \hspace{\stretch{1}}(2.4c)

Demanding a symmetrical solution kills off the pressure term in 2.4a, and we can proceed with separation of variables to find the allowable solutions for the velocity before imposing our boundary value conditions. With $u = R(r) Z(z)$ we have

\begin{aligned}\frac{1}{{r R}} ( R' + r R'' ) - \frac{1}{r^2} = -\frac{Z''}{Z} = -\frac{1}{{a^2}} = \frac{1}{{b^2}}\end{aligned} \hspace{\stretch{1}}(2.5)

Should we pick a negative or a positive constant here (ie: trig or hyperbolic functions for $Z$)? Allowing for both temporarily, we find for $R$

\begin{aligned}0 = r R' + r^2 R'' + R \left( -1 + \frac{r^2}{a^2} \right) \end{aligned} \hspace{\stretch{1}}(2.6a)

\begin{aligned}0 = r R' + r^2 R'' + R \left( -1 - \frac{r^2}{b^2} \right).\end{aligned} \hspace{\stretch{1}}(2.6b)

Using Mathematica (phy454/coffeeCupWithBottom.cdf) to look up the solutions, we find that this was a Bessel equation with solutions

\begin{aligned}R(r) \in \text{span} \{ J_1(r/a), Y_1(r/a) \} \end{aligned} \hspace{\stretch{1}}(2.7a)

\begin{aligned}R(r) \in \text{span} \{ J_1(i r/b), Y_1(i r/b) \}.\end{aligned} \hspace{\stretch{1}}(2.7b)

However, the second set are not real valued for $r > 0$. This means we want the hyperbolic solutions for $Z$. Because we also have a boundary value constrain of $u = 0$ on the bottom of the cup, we can pick only hyperbolic sine solutions. As before, we’ll “stir” the coffee with a cylinder at radius $s$ (with cup radius $R$), our solution has the form

\begin{aligned}u(r, z, 0) = \left\{\begin{array}{l l}s \Omega\sum C_a J_1(r/a) \sinh(z/a)& \quad \mbox{latex r < s} \\ s \Omega\sum \left( D_a J_1(r/a) + E_a Y_1(r/a) \right) \sinh(z/a) & \quad \mbox{$r \in [s, R]$} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.8)

Observe that, regardless of $a$ we have $u(0, z, 0) = 0$ because $J_1(0) = 0$. We also can’t scale $a$ as $\lambda_i/s$ (where $\lambda_i$ are the zeros of $J_1$) when the stirring isn’t at the edge since we need $u(s, z, 0) = \Omega s$. That suggests that we probably want to write our system as

\begin{aligned}u(r, z, 0) = \left\{\begin{array}{l l}s \Omega\sum C_i J_1(\lambda_i r/R) \sinh(\lambda_i z/R)& \quad \mbox{latex r < s} \\ s \Omega\sum \left( D_i J_1(\lambda_i r/R) + E_i Y_1(\lambda_i r/R) \right) \sinh(\lambda_i z/R) & \quad \mbox{$r \in [s, R]$} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.9)

The boundary value constraints of $u(r = s) = \Omega s$ and $u(r = R) = 0$ require the solution of

\begin{aligned}1 = \sum C_i J_1(\lambda_i s/R) \sinh(\lambda_i z/R) \end{aligned} \hspace{\stretch{1}}(2.10a)

\begin{aligned}1 = \sum \left( D_i J_1(\lambda_i s/R) + E_i Y_1(\lambda_i s/R) \right) \sinh(\lambda_i z/R)\end{aligned} \hspace{\stretch{1}}(2.10b)

\begin{aligned}0 = \sum E_i Y_1(\lambda_i) \sinh(\lambda_i z/R).\end{aligned} \hspace{\stretch{1}}(2.10c)

We know how to solve a system like 2.10a if we did not have the $\sinh$ term in the mix. Can we look for a numerical (least squares?) solution to this more general problem. Will it be possible to find something that’s a good fit regardless of the value of $z$?

We have other troubles too. We can’t simultaneously satisfy the boundary value condition of $u(s, z) = \Omega s$ and also $u(s, 0) = 0$. Should we attempt something like a least squares solution and start going near $z = 0$ the small $\sinh$ values near there will start forcing the $C_i$‘s arbitrarily high.

We have to somehow modify the no-slip condition so that when we have a moving interface in contact with a static interface we somehow deal with the fact that the no-slip constraints cannot simultaneously require the velocity to match both the moving and static interfaces at that point of contact.