Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Inclined flow without constant height assumption. Setting up the problem, but not actually solving it.

Posted by peeterjoot on March 18, 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.)]


In an informal discussion after class, it was claimed that the steady state flow down a plane would have constant height, unless you bring surface tension effects into the mix. Part of that statement just doesn’t make sense to me. Consider the forces acting on the fluid in the figure (\ref{fig:inclinedFlowWithoutConstantHeightAssumption:inclinedFlowWithoutConstantHeightAssumptionFig1})

\caption{Gravitational force components acting on fluid flowing down a plane.}


In the inclined reference frame we have a component of the force acting downwards (in the negative y-axis direction), and have a component directed down the x-axis. Wouldn’t this act to both push the fluid down the plane and push part of the fluid downwards? I’d expect this to introduce some vorticity as depicted.

While we are just about to start covered surface tension, perhaps this is just allowing the surface to vary, and then solving the Navier-Stokes equations that result. Let’s try setting up the Navier-Stokes equation for steady state viscous flow down a plane without any assumption that the height is constant and see how far we can get.

The equations of motion.

We’ll use the same coordinates as before, with the directions given as in figure (\ref{fig:inclinedFlowWithoutConstantHeightAssumption:inclinedFlowWithoutConstantHeightAssumptionFig2}). However, this time, we let the height h(x) of the fluid at any distance x down the plane from the initial point vary.

\caption{Diagram of coordinates for inclined flow problem.}


For viscous incompressible flow down the plane our equations of motion are


\begin{aligned}\rho \frac{\partial {u}}{\partial {t}} + \rho (u \partial_x + v \partial_y) u = -\partial_x p + \mu \left(\partial_{xx} + \partial_{yy}\right) u + \rho g \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.1a)

\begin{aligned}\rho \frac{\partial {v}}{\partial {t}} + \rho (u \partial_x + v \partial_y) v = -\partial_y p + \mu \left(\partial_{xx} + \partial_{yy}\right) v - \rho g \cos\alpha \end{aligned} \hspace{\stretch{1}}(2.1b)

\begin{aligned}0 = -\partial_z p \end{aligned} \hspace{\stretch{1}}(2.1c)

\begin{aligned}0 = \partial_x u + \partial_y v.\end{aligned} \hspace{\stretch{1}}(2.1d)


Now, can we kill the time dependent term? Even allowing for u to vary with x and introducing a non-horizontal flow component, I think that we can. If the flow at x = 0 is constant, not varying at all with time, I think it makes sense that we will have no time dependence in the flow for x \ne 0. So, I believe that our starting point is as above, but with the time derivatives killed off. That is


\begin{aligned}\rho (u \partial_x + v \partial_y) u = -\partial_x p + \mu \left(\partial_{xx} + \partial_{yy}\right) u + \rho g \sin\alpha \\ \end{aligned} \hspace{\stretch{1}}(2.2a)

\begin{aligned}\rho (u \partial_x + v \partial_y) v = -\partial_y p + \mu \left(\partial_{xx} + \partial_{yy}\right) v - \rho g \cos\alpha \\ \end{aligned} \hspace{\stretch{1}}(2.2b)

\begin{aligned}0 = -\partial_z p \\ \end{aligned} \hspace{\stretch{1}}(2.2c)

\begin{aligned}0 = \partial_x u + \partial_y v.\end{aligned} \hspace{\stretch{1}}(2.2d)


These don’t look particularily easy to solve, and we haven’t even set up the boundary value constraints yet. Let’s do that next.

The boundary value constraints.

One of out constraints is the no-slip condition for the velocity components at the base of the slope

\begin{aligned}\boxed{u(x, 0) = v(x, 0) = 0.}\end{aligned} \hspace{\stretch{1}}(3.3)

We should also have a zero tangential component to the traction vector at the interface. We need to consider some geometry, and refer to figure (\ref{fig:inclinedFlowWithoutConstantHeightAssumption:inclinedFlowWithoutConstantHeightAssumptionFig3})

\caption{Differential vector element.}


A position vector on the surface has the value

\begin{aligned}\mathbf{r} = x \hat{\mathbf{x}} + h \hat{\mathbf{y}}\end{aligned} \hspace{\stretch{1}}(3.4)

so that a differential element on the surface, tangential to the surface is proportional to

\begin{aligned}d\mathbf{r} = \left( \hat{\mathbf{x}} + \frac{d{{h}}}{dx} \hat{\mathbf{y}} \right) dx\end{aligned} \hspace{\stretch{1}}(3.5)

so our unit tangent vector in the direction depicted in the figure is

\begin{aligned}\hat{\boldsymbol{\tau}} = \frac{1}{{\sqrt{1 + (h')^2}}} \left( 1, h' \right).\end{aligned} \hspace{\stretch{1}}(3.6)

The outwards facing normal has a value, up to a factor of plus or minus one, of

\begin{aligned}\hat{\mathbf{n}} = \frac{1}{{\sqrt{1 + (h')^2}}} \left( h', -1 \right).\end{aligned} \hspace{\stretch{1}}(3.7)

We can fix the orientation by considering the unit bivector

\begin{aligned}\hat{\boldsymbol{\tau}} \wedge \hat{\mathbf{n}} &= \frac{1}{{1 + (h')^2}} \left( 1, h' \right) \wedge\left( h', -1 \right) \\ &=\begin{vmatrix}1 & h' \\ h' & -1\end{vmatrix}\mathbf{e}_1 \mathbf{e}_2 \\ &=( -1 - (h')^2 )\mathbf{e}_1 \mathbf{e}_2.\end{aligned}

So we really want the other orientation

\begin{aligned}\hat{\mathbf{n}} = \frac{1}{{\sqrt{1 + (h')^2}}} \left( -h', 1 \right).\end{aligned} \hspace{\stretch{1}}(3.8)

Our traction vector relative to the normal \hat{\mathbf{n}} is

\begin{aligned}\mathbf{t} &= \mathbf{e}_i T_{ij} n_j \\ &= \mathbf{e}_i \left( -p \delta_{ij}+ \mu e_{ij}\right) n_j \\ &=-p \hat{\mathbf{n}} + \mu \mathbf{e}_i e_{ij} n_j\end{aligned}

So the component in the tangential direction is

\begin{aligned}\mathbf{t} \cdot \hat{\boldsymbol{\tau}} &=-p \not{{\hat{\mathbf{n}} \cdot \hat{\boldsymbol{\tau}}}} + \mu \mathbf{e}_i e_{ij} n_j tau_i \\ &=\frac{\mu}{1 + (h')^2}\begin{bmatrix}1 & h' \end{bmatrix}\begin{bmatrix}e_{11} & e_{12} \\ e_{21} & e_{22}\end{bmatrix}\begin{bmatrix}-h'  \\ 1\end{bmatrix} \\ &=\frac{\mu}{1 + (h')^2}\begin{bmatrix}1 & h' \end{bmatrix}\begin{bmatrix}-h' e_{11} + e_{12} \\ -h' e_{21} + e_{22}\end{bmatrix} \\ &=\frac{\mu}{1 + (h')^2}\left(-h' e_{11} + e_{12} + h'( -h' e_{21} + e_{22} )\right) \\ &=\frac{\mu}{1 + (h')^2}\left(h' (e_{22} - e_{11} )+ e_{12} (1 - (h')^2 )\right) \end{aligned}

Our strain tensor components, for a general 2D flow, are

\begin{aligned}e_{11} &= \frac{\partial {u}}{\partial {x}} \\ e_{22} &= \frac{\partial {v}}{\partial {y}} \\ e_{12} &= \frac{1}{{2}} \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right).\end{aligned} \hspace{\stretch{1}}(3.9)

So, a zero tangential traction vector component at the interface requires

\begin{aligned}\boxed{0 = h' \left( {\left.{{\left( \frac{\partial {v}}{\partial {y}} - \frac{\partial {u}}{\partial {x}} \right) + \frac{1}{{2}} \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right)}}\right\vert}_{{y = h}}\right)(1 - (h')^2 ).}\end{aligned} \hspace{\stretch{1}}(3.12)

What is the normal component of the traction vector at the interface? We can calculate

\begin{aligned}\mathbf{t} \cdot \hat{\mathbf{n}} &=-p \hat{\mathbf{n}} \cdot \hat{\mathbf{n}} + \mu \mathbf{e}_i e_{ij} n_j n_i \\ &=-p+\frac{\mu}{1 + (h')^2}\begin{bmatrix}-h' & 1\end{bmatrix}\begin{bmatrix}-h' e_{11} + e_{12} \\ -h' e_{21} + e_{22}\end{bmatrix} \\ &=-p+\frac{\mu}{1 + (h')^2}\left(-h'(-h' e_{11} + e_{12}) -h' e_{21} + e_{22}\right)  \\ &=-p+\frac{\mu}{1 + (h')^2}\left(-2 h' e_{12} + (h')^2 e_{11} + e_{22}\right)\end{aligned}

So this component of the traction vector is

\begin{aligned}\mathbf{t} \cdot \hat{\mathbf{n}} =-p+\frac{\mu}{1 + (h')^2}\left(- h' \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right)+ (h')^2 \frac{\partial {u}}{\partial {x}}  + \frac{\partial {v}}{\partial {y}} \right)\end{aligned} \hspace{\stretch{1}}(3.13)

For the purely recilinear flow, with h' = 0 and v = 0, as a a consequence of NS and our assumptions, all but the pressure portion of this component of the traction vector was zero. The force balance equation for the interface was therefore just a matching of the pressure with the external (ie: air) pressure.

In this more general case we have the same thing, but the non-pressure portions of the traction vector are all zero in the medium. Outside of the fluid (in the air say), we have assumed no motion, so this force balance condition becomes

\begin{aligned}{\left.{{\mathbf{t} \cdot \hat{\mathbf{n}}}}\right\vert}_{{\text{fluid}}}={\left.{{\mathbf{t} \cdot \hat{\mathbf{n}}}}\right\vert}_{{\text{air}}}.\end{aligned} \hspace{\stretch{1}}(3.14)

Again assuming no motion of the air, with air pressure p_A, this is

\begin{aligned}\boxed{-p(x, h)+{\left.{{\frac{\mu}{1 + (h')^2}\left(- h' \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right)+ (h')^2 \frac{\partial {u}}{\partial {x}}  + \frac{\partial {v}}{\partial {y}} \right)}}\right\vert}_{{y = h}}= -p_A}\end{aligned} \hspace{\stretch{1}}(3.15)

Observe that for the horizontal flow problem, where h' = 0 and v = 0, this would have been nothing more than a requirement that p(h) = p_A, but now that we introduce downwards motion and allow the height to vary, the pressure matching condition becomes a much more complex beastie.

Laplacian of Pressure and Vorticity.

Supposing that we are neglecting the non-linear term of the Navier-Stokes equation. For incompressible steady state flow, without any external forces, we would then have


\begin{aligned}0 = -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u} \end{aligned} \hspace{\stretch{1}}(4.16a)

\begin{aligned}0 = \boldsymbol{\nabla} \cdot \mathbf{u}\end{aligned} \hspace{\stretch{1}}(4.16b)


How do we actually solve this beastie?

Separation of variables?

Considering this in 2D, assuming no z-dependence, with \mathbf{u} = \mathbf{u}(x, y) = (u, v) we have

\begin{aligned}0 &= -\partial_x p + \mu (\partial_{xx} + \partial_{yy} )u \\ 0 &= -\partial_y p + \mu (\partial_{xx} + \partial_{yy} )v \\ 0 &= \partial_x u + \partial_y v.\end{aligned} \hspace{\stretch{1}}(4.17)

Attempting separation of variables seems like something reasonable to try. With

\begin{aligned}u &= X(x) Y(y) \\ v &= R(x) S(y) \\ p &= P(x) Q(y)\end{aligned} \hspace{\stretch{1}}(4.20)

we get

\begin{aligned}0 &= -P' Q + \mu (X'' Y + X Y'') \\ 0 &= -P Q' + \mu (R'' S + R S'') \\ 0 &= X' Y + R S'\end{aligned} \hspace{\stretch{1}}(4.23)

I couldn’t find a way to substitute any of these into the other that would allow me to separate them, but perhaps I wasn’t clever enough.

In terms of vorticity?

The idea of substuting the zero divergence equation \boldsymbol{\nabla} \cdot \mathbf{u} will clearly lead to something a bit simpler. Treating the Laplacian as a geometric (clifford) product of two gradients we have

\begin{aligned}0 &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u} \\ &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla} ( \boldsymbol{\nabla} \mathbf{u} ) \\ &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla} ( \not{{\boldsymbol{\nabla} \cdot \mathbf{u}}} + \boldsymbol{\nabla} \wedge \mathbf{u} ) \\ &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla} ( \boldsymbol{\nabla} \wedge \mathbf{u} ) \\ &= \boldsymbol{\nabla} \left( -p + \mu \boldsymbol{\nabla} \wedge \mathbf{u} \right)\end{aligned}

Writing out the vorticity (bivector) in component form, and writing i = \mathbf{e}_1 \wedge \mathbf{e}_2 = \mathbf{e}_1 \mathbf{e}_2 for the 2D pseudoscalar, we have

\begin{aligned}\boldsymbol{\nabla} \wedge \mathbf{u} &= (\mathbf{e}_1 \partial_x + \mathbf{e}_2 \partial_y) \wedge (\mathbf{e}_1 u + \mathbf{e}_2 v) \\ &= \mathbf{e}_1 \mathbf{e}_2 (\partial_x v - \partial_y u ) \\ &= i (\partial_x v - \partial_y u )\end{aligned}

It seems natural to write

\begin{aligned}\Theta = \partial_x v - \partial_y u,\end{aligned} \hspace{\stretch{1}}(4.26)

so that Navier-Stokes takes the form

\begin{aligned}0 = \boldsymbol{\nabla} (-p + i \mu \Theta ).\end{aligned} \hspace{\stretch{1}}(4.27)

Operating on this from the left with another gradient we find that this combination of pressure and vorticity must satisfy the following multivector Laplacian equation

\begin{aligned}0 = \boldsymbol{\nabla}^2 (-p + i \mu \Theta ).\end{aligned} \hspace{\stretch{1}}(4.28)

However, note that \boldsymbol{\nabla}^2 is a scalar operator, and this zero identity has both scalar and pure bivector components. Both must separately equal zero


\begin{aligned}0 = \boldsymbol{\nabla}^2 p \end{aligned} \hspace{\stretch{1}}(4.29a)

\begin{aligned}0 = \boldsymbol{\nabla}^2 \Theta.\end{aligned} \hspace{\stretch{1}}(4.29b)


Note that we can obtain 4.29 much more directly, if we know that’s what we want to do. Just operate on 4.16a with the gradient from the left right off the bat. We find

\begin{aligned}0 &= -\boldsymbol{\nabla}^2 p + \mu \boldsymbol{\nabla}^3 \mathbf{u} \\ &= -\boldsymbol{\nabla}^2 p + \mu \boldsymbol{\nabla}^2 (\boldsymbol{\nabla} \mathbf{u}) \\ &= -\boldsymbol{\nabla}^2 p + \mu \boldsymbol{\nabla}^2 (\boldsymbol{\nabla} \wedge \mathbf{u}) \\ \end{aligned}

Again, we have a multivector equation scalar and bivector parts, that must separately equal zero. With the magnitude of the vorticity \Theta given by 4.26, we once again obtain 4.29. This can be done in plain old vector algebra as well by operating on the left not by the gradient, but separately with a divergence and curl operator.

Pressure and vorticity equations with the non-linear term retained.

If we add back in our body force, and assume that it is constant (i.e. gravity), then this this will get killed with the application of the gradient. We’ll still end up with one Laplacian for pressure, and one for vorticity. That’s not the case for the inertial (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} term of Navier-Stokes. Let’s take the divergence and curl of this and see how we have to modify the Laplacian equations above.

Starting with the divergence, with summation implied over repeated indexes, we have

\begin{aligned}\boldsymbol{\nabla} \cdot ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u})&=\partial_k ( \mathbf{u} \cdot \boldsymbol{\nabla} u_k ) \\ &=\partial_k ( u_m \partial_m u_k ) \\ &=(\partial_k u_m) (\partial_m u_k ) + u_m \partial_m \partial_k u_k \\ &=\sum_k (\partial_k u_k)^2+\sum_{k \ne m} (\partial_k u_m) (\partial_m u_k ) + (\mathbf{u} \cdot \boldsymbol{\nabla}) ( \boldsymbol{\nabla} \cdot \mathbf{u} )\end{aligned}

So we have

\begin{aligned}\boldsymbol{\nabla} \cdot ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u})=\sum_k (\partial_k u_k)^2+2 \sum_{k < m} (\partial_k u_m) (\partial_m u_k ) + (\mathbf{u} \cdot \boldsymbol{\nabla}) ( \boldsymbol{\nabla} \cdot \mathbf{u} )\end{aligned} \hspace{\stretch{1}}(4.30)

We are working with the \boldsymbol{\nabla} \cdot \mathbf{u} = 0 incompressibility assumption so we kill off the last term. Our velocity ends up introducing a non-homogeneous forcing term to the Laplacian pressure equation and we’ve got something trickier to solve

\begin{aligned}\boxed{\rho \sum_k (\partial_k u_k)^2+2 \rho \sum_{k < m} (\partial_k u_m) (\partial_m u_k ) =-\boldsymbol{\nabla}^2 p.}\end{aligned} \hspace{\stretch{1}}(4.31)

Now let’s see how our vorticity Laplacian needs to be modified. Taking the curl of the implusive term we have

\begin{aligned}\boldsymbol{\nabla} \wedge ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}) &=\mathbf{e}_k \partial_k \wedge ( u_m \partial_m u_n \mathbf{e}_n ) \\ &=(\mathbf{e}_k \wedge \mathbf{e}_n) \partial_k ( u_m \partial_m u_n ) \\ &=(\mathbf{e}_k \wedge \mathbf{e}_n) \left( (\partial_k u_m) (\partial_m u_n ) +u_m \partial_m \partial_k u_n \right) \\ &=(\mathbf{e}_k \wedge \mathbf{e}_n) \left( (\partial_k u_m) (\partial_m u_n ) + (\mathbf{u} \cdot \boldsymbol{\nabla}) \partial_k u_n\right) \end{aligned}

So we have

\begin{aligned}\boldsymbol{\nabla} \wedge ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}) =(\boldsymbol{\nabla} u_m) \wedge (\partial_m \mathbf{u})+ (\mathbf{u} \cdot \boldsymbol{\nabla}) (\boldsymbol{\nabla} \wedge \mathbf{u})\end{aligned} \hspace{\stretch{1}}(4.32)

Putting things back together, our vorticity equation is

\begin{aligned}\rho (\boldsymbol{\nabla} u_m) \wedge (\partial_m \mathbf{u})+ \rho (\mathbf{u} \cdot \boldsymbol{\nabla}) (\boldsymbol{\nabla} \wedge \mathbf{u})= \mu \boldsymbol{\nabla}^2 (\boldsymbol{\nabla} \wedge \mathbf{u})\end{aligned} \hspace{\stretch{1}}(4.33)

Or, with

\begin{aligned}\boldsymbol{\omega} = \boldsymbol{\nabla} \wedge \mathbf{u},\end{aligned} \hspace{\stretch{1}}(4.34)

this is

\begin{aligned}\boxed{(\boldsymbol{\nabla} u_m) \wedge (\partial_m \mathbf{u}) + (\mathbf{u} \cdot \boldsymbol{\nabla}) \boldsymbol{\omega} = \nu \boldsymbol{\nabla}^2 \boldsymbol{\omega}.}\end{aligned} \hspace{\stretch{1}}(4.35)

It is this and 4.31 that we really have to solve. Before moving on, let’s write out all the non-boundary condition equations in coordinate form for the 2D case that we are interested in here. We have


\begin{aligned}\rho \left( \left(\frac{\partial {u}}{\partial {x}}\right)^2 +\left(\frac{\partial {u}}{\partial {y}}\right)^2 + 2 \frac{\partial {u}}{\partial {y}} \frac{\partial {v}}{\partial {x}} \right) = -\boldsymbol{\nabla}^2 p\end{aligned} \hspace{\stretch{1}}(4.36a)

\begin{aligned}2 \left( \frac{\partial {u}}{\partial {x}} \frac{\partial {v}}{\partial {y}}+\frac{\partial {v}}{\partial {x}} \frac{\partial {v}}{\partial {y}}\right)+ \left( u \frac{\partial {}}{\partial {x}} + v \frac{\partial {}}{\partial {y}} \right) \Theta=\nu \boldsymbol{\nabla}^2 \Theta\end{aligned} \hspace{\stretch{1}}(4.36b)

\begin{aligned}\Theta = \frac{\partial {v}}{\partial {x}} - \frac{\partial {u}}{\partial {y}}\end{aligned} \hspace{\stretch{1}}(4.36c)


Our solution has to satisfy these equations, as well as still satisfying the original Navier-Stokes system 2.2 that includes our gravitational term, and also has to satisfy both of our boundary value constraints 3.12, 3.15, and have u(x, 0) = v(x, 0) = 0. Wow, what a mess! And this is all just the steady state problem. Imagine adding time into the mix too!

Now what?

The strategy that I’d thought to attempt to tackle a problem with is something along the following lines

\item First ignore the non-linear terms. Find solutions for the homogenous vorticity and pressure Laplacian equations that satisfy our boundary value conditions, and use that to find a first solution for h(x).
\item Use this to solve for u and v from the vorticity.

However, after having gotten this far, just setting up the problem for solution, I think that it’s perhaps best not to try to solve it yet, and study some more first. I have a feeling that there are likely a number of techniques that have been developed that will likely be useful to know before I try to plow my way through things. It’s interesting to see just how tricky the equations of motion become when one doesn’t make unrealistic assumptions. I have a feeling that to actually attempt this specific problem, I may very well need a computer and numerical techniques.


Leave a Reply

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

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

Google+ photo

You are commenting using your Google+ 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 )


Connecting to %s

%d bloggers like this: