Peeter Joot's Blog.

Math, physics, perl, and programming obscurity.

«
»

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.

About these ads

This entry was posted on April 29, 2012 at 9:18 pm and is filed under Math and Physics Learning.. Tagged: , , . You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

Leave a Reply

Fill in your details below or click an icon to log in:

Gravatar
WordPress.com Logo

You are commenting using your WordPress.com account. ( Log Out / Change )

Twitter picture

You are commenting using your Twitter account. ( Log Out / Change )

Facebook photo

You are commenting using your Facebook account. ( Log Out / Change )

Cancel

Connecting to %s

«
»
 
Follow

Get every new post delivered to your Inbox.

%d bloggers like this: