Peeter Joot's Blog.

Math, physics, perl, and programming obscurity.

«
»

Steady state inclined flow down a plane of two layers of incompressible viscous fluids.

Posted by peeterjoot on March 4, 2012

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

Motivation.

Here’s a generalization of one of the problems from section 2 of [1], itself a slight variation on what we did in class.

In the previous calculation we did the calculation for two incompressible fluids of the same densities flowing down an inclined plane. Now, let’s generalize this slightly, allowing for different densities.

I’m curious how much the air in the neighborhood of some flowing water gets dragged by that flow. It never occurred to me that this would occur, and I’d like to plug in some numbers and see what the results are. This should be something that can be modeled with two layers like this, one of fluid, one of air of a specific thickness (allowing pressure to vary due to the velocity gradient), and one final layer of air at a fixed pressure (atmospheric). I wouldn’t expect that problem to be much harder than this one, although it may end up being worthwhile to let a computer algebra system do some of the grunt work to solve all the resulting equations.

In the end, when we get to putting in some numbers for this problem, we can probably also get an idea how deep the region where the air gets dragged by the fluid can get.

Steady state inclined flow down a plane of two layers of incompressible viscous unequal density fluids.

Setup of the equations of motion for this system.

Our problem is illustrated in figure (\ref{fig:twoLayerInclinedFlowDifferentDensities:twoLayerInclinedFlowDifferentDensitiesFig1}) with a plane set at angle \alpha, fluid depths of h^{(1)} and h^{(2)} respectively, and viscosities \mu^{(1)} and \mu^{(2)}. We’ll write H = h^{(1)} + h^{(2)}. We have a pair of Navier-Stokes equations to solve

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{twoLayerInclinedFlowFig1}
\caption{Two fluids layers in inclined flow.}
\end{figure}

\begin{aligned}\rho^{(i)} \frac{d{{\mathbf{u}^{(i)}}}}{dt} = \rho^{(i)} \frac{\partial {\mathbf{u}^{(i)}}}{\partial {t}} + \rho^{(i)} (\mathbf{u}^{(i)} \cdot \boldsymbol{\nabla}) \mathbf{u}^{(i)} = - \boldsymbol{\nabla} p^{(i)} + \mu^{(i)} \boldsymbol{\nabla}^2 \mathbf{u}^{(i)} + \mu^{(i)} \boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{u}^{(i)}) + \rho^{(i)} \mathbf{g}\end{aligned} \hspace{\stretch{1}}(2.1)

Our steady state and incompressibility constraints break this into a few independent equations

\begin{aligned}\rho^{(i)} \frac{\partial {\mathbf{u}^{(i)}}}{\partial {t}} &= 0 \\ \boldsymbol{\nabla} \cdot \mathbf{u}^{(i)} &= 0 \\ \rho^{(i)} (\mathbf{u}^{(i)} \cdot \boldsymbol{\nabla}) \mathbf{u}^{(i)} &= - \boldsymbol{\nabla} p^{(i)} + \mu^{(i)} \boldsymbol{\nabla}^2 \mathbf{u}^{(i)} + \rho^{(i)} \mathbf{g}\end{aligned} \hspace{\stretch{1}}(2.2)

Let’s require no components of the flows in the y, or z directions initially. As the equivalent of Newton’s law for fluid flows, conservation of linear momentum requires that for our steady state problem we have u_y = u_z = 0 for the flows in both fluid layers.

FIXME: should go back and think about momentum conservation in the context of fluids. I’m now so used to thinking about this as a symmetry issue from a Lagrangian context. With no Lagrangian here, what would be a robust way treat this in fluids? What is the Lagrangian for a fluid mechanics system? I’d imagine that it would be possible to set up a field Lagrangian with velocity fields.

Our problem is now reduced to a problem in four quantities (two velocities and two pressures). With \mathbf{u}^{(i)} = (u^{(i)}, 0, 0) we can restate Navier-Stokes in coordinate form as

\begin{subequations}

\begin{aligned}\partial_x u^{(i)}(x,y,z) = 0\end{aligned} \hspace{\stretch{1}}(2.5a)

\begin{aligned}\rho^{(i)} (u^{(i)} \partial_x u^{(i)} = - \partial_x p^{(i)} + \mu^{(i)} (\partial_{xx} + \partial_{yy} + \partial_{zz}) u^{(i)} + \rho^{(i)} g \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.5b)

\begin{aligned}0 = - \partial_y p^{(i)} - \rho^{(i)} g \cos\alpha \end{aligned} \hspace{\stretch{1}}(2.5c)

\begin{aligned}0 = - \partial_z p^{(i)} \end{aligned} \hspace{\stretch{1}}(2.5d)

\end{subequations}

In order to solve this, we have eight simultaneous non-linear PDEs, four unknown functions, plus boundary conditions!

What are the boundary conditions? One is the “no-slip” condition, the experimental observation that velocities match at the interfaces. So we should have zero velocity for the fluid lying against the plane, and velocity matching between the fluids. The air above the fluid will also be flowing along at the rate of the uppermost portion of the top layer, but we’ll neglect that effect (i.e. considering two layers of equal density and not three, with one having a separate density). We also have matching of the traction vectors at the interfaces.

Writing this, it occurred to me that I didn’t fully understand what motivated the traction vector matching boundary value condition. Talking to our Prof about this, the matching of the traction vectors at any point can be thought of as an observational issue, but this is also a force balance issue. There is an induced velocity in the direction of the traction vector at any given point. For example, when we have unidirectional flow, we must have no normal component of the traction vector, and only a tangential component, because we have only the tangential flow. It is probably reasonable to think about this roughly as the equivalent of matching both acceleration and velocity at the boundary, but because densities and viscosities vary, we have to match the traction vectors and not the acceleration itself.

Before continuing to solve our Navier-Stokes equations let’s express the condition that the tangential component of the traction vectors match algebraically.

Dropping indexes temporarily, for the normal to the surface \hat{\mathbf{n}} = (n_1, n_2, n_3) = (0, 1, 0) we want to compute

\begin{aligned}T_1 &= \sigma_{1 k} n_k \\ &= \sigma_{1 k} \delta_{2 k} \\ &= \sigma_{1 2} \\ &= -p \not{{\delta_{1 2}}} + 2 \mu e^{1 2} \\ &= \mu \left( \frac{\partial {u_1}}{\partial {y}} + \not{{\frac{\partial {u_2}}{\partial {x}}}} \right) \end{aligned}

So the tangential component of the traction vector is

\begin{aligned}\mathbf{T}^{(i)} = \mu^{(i)} \frac{\partial {u^{(i)}}}{\partial {y}} \hat{\mathbf{x}}.\end{aligned} \hspace{\stretch{1}}(2.6)

As noted above, this is in fact, the only component of the traction vector, since we don’t have any non-horizontal flow.

Our boundary value conditions, what we need in addition to the Navier-Stokes equations of 2.5, to solve our problem, are the matching at any interface of the following conditions

\begin{subequations}

\begin{aligned}u^{(i)} = u^{(j)} \end{aligned} \hspace{\stretch{1}}(2.7a)

\begin{aligned}p^{(i)} = p^{(j)} \end{aligned} \hspace{\stretch{1}}(2.7b)

\begin{aligned}\mu^{(i)} \frac{\partial {u^{(i)}}}{\partial {y}} = \mu^{(i)} \frac{\partial {u^{(j)}}}{\partial {y}}.\end{aligned} \hspace{\stretch{1}}(2.7c)

\end{subequations}

There are actually three interfaces to consider, that of the lower layer liquid with the inclined plane, the interface between the two fluid layers, and the interface between the upper layer fluid and the air above it.

Solving our equations of motion.

Starting with the simplest, the z-coordinate equation, of Navier-Stokes 2.5d, we can conclude that each of the pressures is not a function of z, so that we have

\begin{aligned}p^{(i)} = p^{(i)}(x, y).\end{aligned} \hspace{\stretch{1}}(2.8)

Using this, we can integrate our y-coordinate Navier-Stokes equation 2.5c, to find

\begin{aligned}p^{(i)} = - \rho^{(i)} g y \cos\alpha + f^{(i)}(x)\end{aligned} \hspace{\stretch{1}}(2.9)

At this point we can introduce the first boundary value constraint, that the pressures must match at the interfaces. In particular, on the upper surface, where we have atmospheric pressure p_A our pressure is

\begin{aligned}p^{(2)}(H) = - \rho^{(2)} g H \cos\alpha + f^{(2)}(x) = p_A,\end{aligned} \hspace{\stretch{1}}(2.10)

so f^{(2)} is constant with value

\begin{aligned}f^{(2)}(x) = p_A + \rho^{(2)} g H \cos\alpha,\end{aligned} \hspace{\stretch{1}}(2.11)

which fully determines the density of the upper surface

\begin{aligned}p^{(2)}(y) = \rho^{(2)} g \cos\alpha (H - y) + p_A.\end{aligned} \hspace{\stretch{1}}(2.12)

Matching the pressure between the two layers of fluids we have

\begin{aligned}p^{(1)}(h_1) &= - \rho^{(1)} g h_1 \cos\alpha + f^{(1)}(x) \\ &= p^{(2)}(h_1) \\ &= \rho^{(2)} g \cos\alpha (H - h_1) + p_A \\ &= \rho^{(2)} g h_2 \cos\alpha + p_A,\end{aligned}

so that our undetermined function f^{(1)}(x)

\begin{aligned}f^{(1)}(x) = \left(\rho^{(1)} h_1 + \rho^{(2)} h_2 \right) g \cos\alpha + p_A.\end{aligned} \hspace{\stretch{1}}(2.13)

With the densities not equal, we no longer find that the pressure is dependent only on the total height y, independent of the velocities and viscosities

\begin{aligned}p^{(1)}(y) = g \cos\alpha \left(\rho^{(1)} (h_1 - y) + \rho^{(2)} h_2\right)+ p_A.\end{aligned} \hspace{\stretch{1}}(2.14)

However, this is still a fairly satisfying result. The pressure on the bottom layer is the total pressure due to the layer above it (the contribution due to the total height h_2 of that layer of the fluid). To that we add the pressure at our specific height, a linear function of the difference from the interface above it. Specified piecewise our pressure is now fully determined

\begin{aligned}\boxed{p(y) =\left\{\begin{array}{l l}g \cos\alpha \left(\rho^{(1)} (h_1 - y) + \rho^{(2)} h_2\right)+ p_A & \quad \mbox{latex y H$}\end{array}\right.}\end{aligned} \hspace{\stretch{1}}(2.15)$

Observe that we have the usual \rho g h form in all the terms of the pressure above, just scaled by the cosine of the angle since only a portion of the gravitational force is pushing normally on the fluids.

Having solved for the pressure, we are now set to return to the remaining Navier-Stokes equations 2.5a, and 2.5b for this system. From 2.5a we see that the non-linear term on the LHS of 2.5b is killed and also see that our velocities can only be functions of y and z

\begin{aligned}u^{(i)} = u^{(i)}(y, z).\end{aligned} \hspace{\stretch{1}}(2.16)

While more general solutions can likely be found, we will limit ourselves to looking only for solutions that are functions of y. From our solution to the pressure part of the problem p^{(i)} = p^{(i)}(y), we also see that the pressure term \partial_x p^{(i)} of 2.5b is killed. We are left with just

\begin{aligned}0 &= \mu^{(i)} (\partial_{xx} + \partial_{yy} + \partial_{zz}) u^{(i)} + \rho^{(i)} g \sin\alpha \\ &= \mu^{(i)} \partial_{yy} u^{(i)} + \rho^{(i)} g \sin\alpha \\ &= \mu^{(i)} \frac{d^2}{dy^2} u^{(i)} + \rho^{(i)} g \sin\alpha \end{aligned}

This is directly integrable, and we find for the velocities and traction vectors respectively

\begin{subequations}

\begin{aligned}u^{(i)} = -\frac{\rho^{(i)} g \sin\alpha }{2 \mu^{(i)}} y^2 + A^{(i)} y + B^{(i)}.\end{aligned} \hspace{\stretch{1}}(2.17a)

\begin{aligned}T_x^{(i)} = \mu^{(i)} \frac{\partial {u^{(i)}}}{\partial {y}} = - \rho^{(i)} g y \sin\alpha + \mu^{(i)} A^{(i)} \end{aligned} \hspace{\stretch{1}}(2.17b)

\end{subequations}

The boundary conditions left to exploit are

\begin{aligned}u^{(1)}(0) &= 0 \\ u^{(1)}(h_1) &= u^{(2)}(h_1) \\ T_x^{(1)}(h_1) &= T_x^{(2)}(h_1) \\ T_x^{(2)}(H) &= 0,\end{aligned} \hspace{\stretch{1}}(2.18)

The first is the no-slip condition with the plane. The last is an approximation that assumes the liquid isn’t producing a measurable force on the air above it. The other two are for the interfaces between the two fluids.

From u^{(1)}(0) = 0 we see immediately that we have B^{(1)} = 0. From the traction vector equality in the atmosphere, we have

\begin{aligned}0 = - \rho^{(2)} g H \sin\alpha + \mu^{(2)} A^{(2)},\end{aligned} \hspace{\stretch{1}}(2.22)

or

\begin{aligned}A^{(2)} = \frac{\rho^{(2)} g H \sin\alpha }{ \mu^{(2)} }.\end{aligned} \hspace{\stretch{1}}(2.23)

These reduce the problem to solving for two last integration constants, where our velocities are

\begin{aligned}u^{(1)} &= -\frac{\rho^{(1)} g \sin\alpha }{2 \mu^{(1)}} y^2 + A^{(1)} y \\ u^{(2)} &= \frac{\rho^{(2)} g \sin\alpha }{2 \mu^{(2)}} \left( 2 H y -y^2 \right) + B^{(2)}.\end{aligned} \hspace{\stretch{1}}(2.24)

and our traction vectors are

\begin{aligned}T_x^{(1)} &= - \rho^{(1)} g y \sin\alpha + \mu^{(1)} A^{(1)} \\ T_x^{(2)} &= \rho^{(2)} g \sin\alpha \left( H - y \right).\end{aligned} \hspace{\stretch{1}}(2.26)

Matching both at the interface (y = h_1) gives us

\begin{aligned}-\frac{\rho^{(1)} g \sin\alpha }{2 \mu^{(1)}} h_1^2 + A^{(1)} h_1 &= \frac{\rho^{(2)} g \sin\alpha }{2 \mu^{(2)}} h_1 \left( 2 h_2 + h_1 \right) + B^{(2)} \\ - \rho^{(1)} g h_1 \sin\alpha + \mu^{(1)} A^{(1)} &= \rho^{(2)} g h_2 \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.28)

We find

\begin{aligned}A^{(1)} = \frac{1}{{ \mu^{(1)} }}(\rho^{(1)} h_1 + \rho^{(2)} h_2 )g \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.30)

Let’s substitute this back for our first fluid velocity

\begin{aligned}u^{(1)} &= -\frac{\rho^{(1)} g \sin\alpha }{2 \mu^{(1)}} y^2 + \frac{y}{ \mu^{(1)} }(\rho^{(1)} h_1 + \rho^{(2)} h_2 )g \sin\alpha \\ &=g \sin\alpha \left( -\frac{\rho^{(1)} }{2 \mu^{(1)}} y^2 + \frac{y}{ \mu^{(1)} }(\rho^{(1)} h_1 + \rho^{(2)} h_2 )\right) \\ &=\frac{g y \sin\alpha}{2 \mu^{(1)}} \left( \rho^{(1) } (2 h_1 - y) +2 \rho^{(2)} h_2 \right) \\ \end{aligned}

As a check we see this is consistent with the previous calculation when \rho^{(1)} = \rho^{(2)}. For our final integration constant we now find

\begin{aligned}B^{(2)} &=\frac{g h_1 \sin\alpha}{2 \mu^{(1)}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) -\frac{\rho^{(2)} g \sin\alpha }{2 \mu^{(2)}} h_1 \left( 2 h_2 + h_1 \right) \\ &=\frac{g h_1 \sin\alpha}{2}\left(\frac{1}{{\mu^{(1)}}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) -\frac{\rho^{(2)} }{\mu^{(2)}} \left( 2 h_2 + h_1 \right) \right) \\ &=\frac{g h_1 \sin\alpha}{2 \mu^{(2)}}\left(\frac{\mu^{(2)}}{\mu^{(1)}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) - \rho^{(2)} \left( 2 h_2 + h_1 \right) \right) \\ \end{aligned}

So, finally, we have for the velocities

\begin{subequations}

\begin{aligned}u^{(1)}(y) = \frac{g y \sin\alpha}{2 \mu^{(1)}} \left( \rho^{(1) } (2 h_1 - y) +2 \rho^{(2)} h_2 \right) \end{aligned} \hspace{\stretch{1}}(2.31a)

\begin{aligned}u^{(2)}(y) = \frac{g \sin\alpha }{2 \mu^{(2)}} \left(\rho^{(2)} \left( 2 H y -y^2 \right) + h_1 \left(\frac{\mu^{(2)}}{\mu^{(1)}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) - \rho^{(2)} \left( 2 h_2 + h_1 \right) \right) \right)\end{aligned} \hspace{\stretch{1}}(2.31b)

\end{subequations}

The final result looks reasonable. If the viscosities and densities are equal then we have the same velocity profile in both layers. That makes sense given the equal densities, since there would really be nothing that would then distinguish the two layers.

Numerical application.

To try this out numerically, I’ve created a Mathematica worksheet. The results are fairly suprising, showing either an error in this calculation, an error programming the worksheet, or the folly of even considering a steady state flow of this form for anything that is not extremely viscous. Some validation is required to see what’s up.

References

[1] D.J. Acheson. Elementary fluid dynamics. Oxford University Press, USA, 1990.

About these ads

This entry was posted on March 4, 2012 at 12:42 am 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.

4 Responses to “Steady state inclined flow down a plane of two layers of incompressible viscous fluids.”

  1. Toby Haynes said

    March 7, 2012 at 9:56 am

    My first thought when I saw the title was “Kelvin Helmholtz instability”. Instantly compounded when I saw the Numerical application section. I suspect that any attempt to flow two different densities of fluid down a slope quickly degenerates, even in the non-turbulent phases.

    http://hmf.enseeiht.fr/travaux/CD0001/travaux/optmfn/hi/01pa/hyb72/kh/kh_theo.htm

    Reply

    • peeterjoot said

      March 13, 2012 at 2:37 pm

      Thanks for the reference Toby, we’ve not yet gotten to any stability issues. What’s your background? I was suprised to see a comment on fluid dynamics out of an IBMer.

      Reply

      • Toby Haynes said

        March 13, 2012 at 2:50 pm

        Ah, my distant past… Natural Sciences, specializing in Experimental Physics. Oh, and some mucking about with magnetic fields for clusters of galaxies while I was playing with Observational Astronomy. So I went from watching Large Objects to coding Large Objects (database humour if there are watchers in the galleries).

        Anyway, I wonder if there is a steady state solution to the problem you pose if the two layers were to move at different speeds.

      • peeterjoot said

        March 14, 2012 at 7:52 am

        I think that the solution I found is valid for layers that are thin enough that the height does not vary significantly over the interval of interest. Once the height starts to vary (which I think it should, due to the gravitational force on the fluid), then you’ll start getting a y component to the flow, and the solution is invalid. I’m trying this calculation to see if it works a bit better, but it’s harder, and my Professor has suggested a nice approach to attempt a pertubation solution.

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: