Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

PHY454H1S Continuum mechanics. Problem Set 2. Rectilinear flow with pressure gradient. Two layer flow driven by moving surface.

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


This problem set is as yet ungraded.

Problem Q1.


Imagine a steady rectilinear blood flow of the form \mathbf{u} = u(y) \hat{\mathbf{y}} through a two dimensional artery. It is driven by a constant pressure gradient G = -dp/dx maintained by an external `heart’. The top and bottom walls of the artery are 2h distance apart and the fluid satisfies no-slip boundary conditions at the walls. Assuming that the fluid is Newtonian,

\item Show that the Navier-Stokes equation reduces to

\begin{aligned}\frac{d^2 u}{dy^2} = -\frac{G}{\mu}\end{aligned} \hspace{\stretch{1}}(2.1)

where \mu is the viscosity of the blood.

\item Show that the velocity profile of the fluid inside the artery is a parabolic profile.
\item What is the maximum speed of the fluid? Draw the velocity profile to show where the maximum speed occurs inside the artery.
\item If due to smoking etc., the viscosity of the blood gets doubled, then what should be the new pressure gradient to be maintained by the `heart’ to keep the liquid flux through the artery at the same level as the non-smoking one?

Solution. Part 1. Navier-Stokes equation for the system.

The Navier-Stokes equation, for an incompressible unidirectional fluid \mathbf{u} = (u, 0, 0), assuming that there is no z dependence, takes the form


\begin{aligned}\rho \frac{\partial {u}}{\partial {t}} + u \frac{\partial {u}}{\partial {x}} = - \frac{\partial {p}}{\partial {x}} + \mu \left( \frac{\partial^2 {{}}}{\partial {{x}}^2} + \frac{\partial^2 {{}}}{\partial {{y}}^2} \right) u\end{aligned} \hspace{\stretch{1}}(2.2a)

\begin{aligned}0 = - \frac{\partial {p}}{\partial {y}} \end{aligned} \hspace{\stretch{1}}(2.2b)

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

\begin{aligned}0 = \frac{\partial {u}}{\partial {x}}.\end{aligned} \hspace{\stretch{1}}(2.2d)


With a steady state assumption we kill the {\partial {u}}/{\partial {t}} term, and 2.2d kills of the x-component of the Laplacian and our non-linear inertial term on the LHS, leaving just


\begin{aligned}0 = - \frac{\partial {p}}{\partial {x}} + \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2} \end{aligned} \hspace{\stretch{1}}(2.3a)

\begin{aligned}0 = - \frac{\partial {p}}{\partial {y}} \end{aligned} \hspace{\stretch{1}}(2.3b)

\begin{aligned}0 = - \frac{\partial {p}}{\partial {z}}.\end{aligned} \hspace{\stretch{1}}(2.3c)


With \frac{\partial {p}}{\partial {z}} = \frac{\partial {p}}{\partial {y}} = 0, we have \frac{\partial {p}}{\partial {x}} = dp/dx = -G, so 2.3a is reduced to

\begin{aligned}0 = G + \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.4)

Finally, since we have an assumption of no z-dependence ({\partial {u}}/{\partial {z}} = 0) and from the incompressibility assumption 2.2d ({\partial {u}}/{\partial {x}} = 0), we have

\begin{aligned}\frac{\partial^2 {{u}}}{\partial {{y}}^2} = \frac{d^2 u}{dy^2} = -\frac{G}{\mu},\end{aligned} \hspace{\stretch{1}}(2.5)

as desired.

Solution. Part 2. Velocity profile.

Integrating 2.5 twice, we have

\begin{aligned}u = -\frac{G}{2 \mu} y^2 + A y + B.\end{aligned} \hspace{\stretch{1}}(2.6)

Application of the no-slip boundary value condition u(\pm h) = 0, we have

\begin{aligned}0 &= -\frac{G}{2 \mu} h^2 + A h + B \\ 0 &= -\frac{G}{2 \mu} h^2 - A h + B \end{aligned} \hspace{\stretch{1}}(2.7)

Adding and subtracting these, we find


\begin{aligned}A = 0 \end{aligned} \hspace{\stretch{1}}(2.9a)

\begin{aligned}B = \frac{G}{2 \mu} h^2,\end{aligned} \hspace{\stretch{1}}(2.9b)


so the velocity is given by the parabolic function

\begin{aligned}u(y) = \frac{G}{2 \mu} \left( h^2 - y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.10)

Solution. Part 3. Maximum speed of the flow.

It is clear that the maximum speed of the fluid is found at y = 0

\begin{aligned}u(0) = \frac{G h^2}{2 \mu}\end{aligned} \hspace{\stretch{1}}(2.11)

The velocity profile for this flow is drawn in figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig1}).

\caption{Velocity profile for 1D constant pressure gradient steady state flow.}

Solution. Part 4. Effects of viscosity doubling.

With our velocity being dependent on the G/\mu ratio, it is clear, even without calculating the flux, that we will need twice the pressure gradient if the viscosity is doubled to maintain the same flux through the artery and veins. To demonstrate this more thoroughly we can calculate this mass flux. For an element of mass leaving a portion of the conduit, bounded by the plane normal to \hat{\mathbf{x}} we have

\begin{aligned}\frac{dm}{dt} &= \rho \frac{dV}{dt} \\ &= \rho dz dy \mathbf{u} \cdot \hat{\mathbf{x}}\end{aligned}

Integrating this over a width \Delta z, our flux through the plane is

\begin{aligned}\text{Flux} &= \int_0^{\Delta z} dz\int_{-h}^h dy \frac{G}{2 \mu} \left( h^2 - y^2 \right) \\ &= \Delta z \frac{G}{2 \mu} \evalrange{ \left( h^2 y - \frac{1}{{3}} y^3 \right) }{-h}{h} \\ &=\Delta z \frac{G h^3}{\mu} \left( 1 - \frac{1}{{3}} \right)  \\ &=\Delta z \frac{2 G h^3}{3 \mu}.\end{aligned}

Doubling the blood viscosity for our smoker, our respective fluxes are

\begin{aligned}\text{Flux}_{\text{smoker}} &= \Delta z \frac{2 G_{\text{smoker}} h^3}{3 (2 \mu)}  \\ \text{Flux}_{\text{non-smoker}} &= \Delta z \frac{2 G h^3}{3 \mu}.\end{aligned} \hspace{\stretch{1}}(2.12)

Demanding equality before and after smoking we find

\begin{aligned}G_{\text{smoker}} = 2 G.\end{aligned} \hspace{\stretch{1}}(2.14)

where G is the magnitude of the pressure gradient before the bad habits kicked in. The smoker’s poor little heart (soon to be a big overworked and weak heart) has to generate pressure gradients that are twice as big to get the same quantity of blood distributed through the body.

Problem Q2.


Consider steady simple shearing flow with no imposed pressure gradient (G = 0) of a two layer fluid with viscosity

\begin{aligned}\mu =\left\{\begin{array}{l l}\mu^{(1)} & \quad \mbox{latex -h < y < 0,$} \\ \mu^{(2)} & \quad \mbox{0 < y < h.}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(3.15)$

The boundary conditions are no-slip at the lower plate (y = -h) and at y = 0. The top plate is moving with a velocity -U at y = h and fluid is sticking to it. using the continuity of tangential (shear) stress at the interface (y = 0)

\item Derive the velocity profile of the two fluids.
\item Calculate the maximum speed.
\item Calculate the mean speed.
\item Calculate the flux (the volume flow rate.)
\item Calculate the tangential force (per unit width) F_x on the strip 0 \le x \le L of the wall y = -h.
\item Calculate the tangential force (per unit width) F_x^0 on the strip 0 \le x \le L at the interface y = 0 by the top fluid on the lower fluid.

Part 1. Derive the velocity profile of the two fluids.

With coordinates referring to figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig2}), our steady flow for layers 1 and 2 has the form

\caption{Two layer flow induced by moving wall.}


\begin{aligned}0 = - \frac{\partial {p}}{\partial {x}} + \mu^{(i)} \frac{\partial^2 {{u^{(i)}}}}{\partial {{y}}^2}\end{aligned} \hspace{\stretch{1}}(3.16a)

\begin{aligned}0 = - \frac{\partial {p}}{\partial {y}} \end{aligned} \hspace{\stretch{1}}(3.16b)

\begin{aligned}0 = - \frac{\partial {p}}{\partial {z}},\end{aligned} \hspace{\stretch{1}}(3.16c)


as we found in Q1. Only the boundary value conditions and the driving pressure are different here. In this problem and the next, we have constant pressure gradients dp/dx = -G to deal with, so we really have just the pair of equations

\begin{aligned}0 = G + \mu^{(i)} \frac{d^2 u^{(i)}(y)}{dy^2},\end{aligned} \hspace{\stretch{1}}(3.17)

to solve. For this Q2 problem we have G = 0, so the algebra to match our boundary value constraints becomes a bit easier. Our boundary value constraints are


\begin{aligned}u^{(1)}(-h) = 0 \\ \end{aligned} \hspace{\stretch{1}}(3.18a)

\begin{aligned}u^{(2)}(h) = -U \\ \end{aligned} \hspace{\stretch{1}}(3.18b)

\begin{aligned}u^{(1)}(0) = u^{(2)}(0),\end{aligned} \hspace{\stretch{1}}(3.18c)


plus one more to match the tangential components of the traction vector with respect to the normal \hat{\mathbf{n}} = (0, 1, 0). The components of that traction vector are

\begin{aligned}t_i &= \left( -p \delta_{ij} + 2 \mu e_{ij} \right) n_j \\ &= \left( -p \delta_{ij} + 2 \mu e_{ij} \right) \delta_{2j} \\ &= -p \delta_{i2} + 2 \mu e_{i2},\end{aligned}

but we are only interested in the horizontal component t_1 which is

\begin{aligned}t_1 &= -p \delta_{12} + 2 \mu e_{12} \\ &=2 \mu \frac{1}{{2}} \left( \frac{\partial {y}}{\partial {u}}+\not{{\frac{\partial {x}}{\partial {v}}}}\right).\end{aligned}

So the matching the tangential components of the traction vector at the interface gives us our last boundary value constraint

\begin{aligned}{\left.{{\mu^{(1)} \frac{\partial {y}}{\partial {u}}}}\right\vert}_{{y = 0}} = {\left.{{\mu^{(2)} \frac{\partial {y}}{\partial {u}}}}\right\vert}_{{y = 0}},\end{aligned} \hspace{\stretch{1}}(3.19)

and we are ready to do our remaining bits of algebra. We wish to solve the pair of equations

\begin{aligned}u^{(1)} &= A^{(1)} y + B^{(1)} \\ u^{(2)} &= A^{(2)} y + B^{(2)},\end{aligned} \hspace{\stretch{1}}(3.20)

for the four integration constants A^{(i)} and B^{(i)} using our boundary value constraints. The linear system to solve is

\begin{aligned}0 &= -A^{(1)} h + B^{(1)} \\ -U &= A^{(2)} h + B^{(2)} \\ B^{(1)} &= B^{(2)} \\ \mu^{(1)} A^{(1)} &= \mu^{(2)} A^{(2)}.\end{aligned} \hspace{\stretch{1}}(3.22)

With B = B^{(i)}, we have

\begin{aligned}0 &= -A^{(1)} h + B \\ -U &= \frac{\mu^{(1)}}{\mu^{(2)}} A^{(1)} h + B \end{aligned} \hspace{\stretch{1}}(3.26)

Subtracting these to solve for A^{(1)} we find

\begin{aligned}-U = h A^{(1)} \left( \frac{\mu^{(1)}}{\mu^{(2)}} + 1 \right).\end{aligned} \hspace{\stretch{1}}(3.28)

This gives us everything we need

\begin{aligned}A^{(1)} &= -\frac{U \mu^{(2)}}{h(\mu^{(1)} + \mu^{(2)})} \\ A^{(2)} &= -\frac{U \mu^{(1)}}{h(\mu^{(1)} + \mu^{(2)})} \\ B^{(1)} = B^{(2)} &= -\frac{U \mu^{(2)}}{\mu^{(1)} + \mu^{(2)}},\end{aligned} \hspace{\stretch{1}}(3.29)

Referring back to 3.20 our velocities are

\begin{aligned}\boxed{\begin{aligned}u^{(1)} &= -\frac{U \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \left( 1 + \frac{y}{h} \right) \\ u^{(2)} &= -\frac{U \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \left( 1 + \frac{\mu^{(1)}}{\mu^{(2)}} \frac{y}{h} \right) \end{aligned}}\end{aligned} \hspace{\stretch{1}}(3.32)

Checking, we see at a glance we see that we have u^{(2)}(h) = -U, u^{(1)}(-h) = 0, u^{(1)}(0) = u^{(2)}(0), and {\left.{{\mu^{(1)} du^{(1)}/dy}}\right\vert}_{{y=0}} = \mu^{(2)} {\left.{{du^{(2)}/dy}}\right\vert}_{{y=0}} as desired.

As an example, let’s add some numbers. With mercury and water in layers {(1)} and {(2)} respectively, we have

\begin{aligned}\mu^{(1)} &= 0.001526 \quad \text{Pa-s} \\ \mu^{(2)} &= 0.00089 \quad \text{Pa-s} \end{aligned} \hspace{\stretch{1}}(3.33)

so that our velocity is

\begin{aligned}u(y) = \left\{\begin{array}{l l}-1.37 U \left(1 + \frac{y}{h}\right) & \quad \mbox{latex y \in [-h, 0]$} \\ -1.37 U \left(1 + 1.7 \frac{y}{h}\right) & \quad \mbox{y \in [0, h]}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(3.35)$

This is plotted with h = U = 1 in figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig3})
\caption{Two layer shearing flow with water over mercury.}

Part 2. Calculate the maximum speed.

With u^{(1)}(-h) = 0, and u^{(1)} linearly decreasing, then u^{(2)} linearly decreasing further from the value at y = 0, it is clear that the maximum speed, no matter the viscosities of the fluids, is on the upper moving interface. This maximum takes the value {\left\lvert{u^{(2)}(h)}\right\rvert} = U.

Part 3. Calculate the mean speed.

As linear functions the average speeds of the respective fluids fall on the midpoints y = \pm h/2. These are

\begin{aligned}\left\langle{{u^{(1)}}}\right\rangle &= -\frac{U \mu^{(2)} }{ 2 (\mu^{(1)} + \mu^{(2)})} \\ \left\langle{{u^{(2)}}}\right\rangle &= -\frac{U \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \left( 1 + \frac{\mu^{(1)}}{2 \mu^{(2)}} \right) \end{aligned} \hspace{\stretch{1}}(3.36)

Averaging these two gives us the overall average, so we find

\begin{aligned}\left\langle{{u(y)}}\right\rangle = -\frac{U }{ 4 (\mu^{(1)} + \mu^{(2)})} \left( 3 \mu^{(2)} + \mu^{(1)} \right) \end{aligned} \hspace{\stretch{1}}(3.38)

Part 4. Calculate the flux

We can calculate the volume flux, much like the mass flux (although the mass flux seems a more sensible quantity to calculate). Looking at the rate of change of an element of fluid passing through the y-z plane we have

\begin{aligned}\frac{dV}{dt} = dy dz \mathbf{u} \cdot \hat{\mathbf{x}}\end{aligned} \hspace{\stretch{1}}(3.39)

Integrating over the total height, for a width \Delta z we have

\begin{aligned}\text{Volume Flux} &= \Delta z \int_{-h}^h u(y) dy \\ &= \Delta z 2h \left\langle{{u}}\right\rangle \\ \end{aligned}

So our volume flux through a width \Delta z is

\begin{aligned}\text{Volume Flux} =-\frac{2 h U \Delta z}{ 4 (\mu^{(1)} + \mu^{(2)})} \left( 3 \mu^{(2)} + \mu^{(1)} \right).\end{aligned} \hspace{\stretch{1}}(3.40)

Part 5. Tangential force on the wall.

From 3.32 the tangential components of our traction vectors are

\begin{aligned}t^{(1)} &= \mu^{(1)} \frac{d u^{(1)}}{dy} \\ &= -\frac{U \mu^{(1)} \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \frac{1}{{h}}\end{aligned}


\begin{aligned}t^{(2)} &= \mu^{(2)} \frac{d u^{(2)}}{dy} \\ &= -\frac{U \mu^{(1)} \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \frac{1}{{h}}\end{aligned}

We see that the tangential component of the traction vector is a constant throughout both fluids. Allowing this force to act on a length L of the lower wall, our force per unit width over that strip is just

\begin{aligned}F = -\frac{U \mu^{(1)} \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \frac{L}{h}.\end{aligned} \hspace{\stretch{1}}(3.41)

The negative value here makes sense since it is acting to push the fluid backwards in the direction of the upper wall motion.

Part 6. Tangential force on the lower fluid.

Due to constant nature of the tangential component of the traction vector shown above, the force per unit width of the upper fluid acting on the lower fluid, is also given by 3.41.

Problem Q3.


If on top of the problem described above a constant pressure gradient G = -dp/dx is applied between the boundaries y = \pm h, describe qualitatively what type of flow profile you would expect in the steady state. Draw the velocity profiles for two cases (i) \mu^{(1)} > \mu^{(2)} (ii) \mu^{(1)} < \mu^{(2)}. Explain your result.


We showed earlier that the Navier-Stokes equations for this Q3 case, where G is non-zero were given by 3.17, which restated is

\begin{aligned}0 = G + \mu^{(i)} \frac{d^2 u^{(i)}(y)}{dy^2}.\end{aligned} \hspace{\stretch{1}}(4.42)

Our solutions will now necessarily be parabolic, of the form

\begin{aligned}u^{(i)}(y) = -\frac{G}{2 \mu^{(i)}} y^2 + A^{(i)} y + B^{(i)},\end{aligned} \hspace{\stretch{1}}(4.43)

with the tangential traction vector components given by

\begin{aligned}t^{(i)} = -G y + A^{(i)} \mu^{(i)}\end{aligned} \hspace{\stretch{1}}(4.44)

The boundary value constants become a bit messier to solve for, and should we wish to do so we’d have to solve the system

\begin{aligned}0 &= -\frac{G}{2 \mu^{(1)}} h^2 - A^{(1)} h + B^{(1)} \\ -U &= -\frac{G}{2 \mu^{(2)}} h^2 + A^{(2)} h + B^{(2)} \\ B^{(1)} &= B^{(2)} \\ A^{(1)} \mu^{(1)} &= A^{(2)} \mu^{(2)} \end{aligned} \hspace{\stretch{1}}(4.45)

Without actually solving this system we should expect that our solution will have the form of our pure shear flow, with parabolas superimposed on these linear flows. For a higher viscosity bottom layer \mu^{(1)} > \mu^{(2)}, this should look something like figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig4}) whereas for the higher viscosity on the top, these would be roughly flipped as in figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig5}).

\caption{Superposition of constant pressure gradient and shear flow solutions (\mu^{(1)} > \mu^{(2)}).}

\caption{Superposition of constant pressure gradient and shear flow solutions (\mu^{(1)} < \mu^{(2)}).}

This superposition can be justified since we have no (\mathbf{u} \cdot \boldsymbol{\nabla})\mathbf{u} term in the Navier-Stokes equations for these systems.

Exact solutions.

The figures above are kind of rough. It’s not actually hard to solve the system above. After some simplification, we get

\begin{aligned}u^{(1)}(y) &= -\frac{\mu^{(2)} U -G h^2}{\mu^{(1)}+\mu^{(2)}}-\frac{y \left(G h^2 (\mu^{(2)} -\mu^{(1)}) +2 \mu^{(1)} \mu^{(2)} U \right)}{2 h \mu^{(1)} \left(\mu^{(1)}+\mu^{(2)}\right)}-\frac{G y^2}{2 \mu^{(1)}} \\ u^{(2)}(y) &= -\frac{\mu^{(2)} U -G h^2}{\mu^{(1)}+\mu^{(2)}}-\frac{y \left(G h^2 (\mu^{(2)} -\mu^{(1)}) +2 \mu^{(1)} \mu^{(2)} U \right)}{2 h \mu^{(2)} \left(\mu^{(1)}+\mu^{(2)}\right)}-\frac{G y^2}{2 \mu^{(2)}}.\end{aligned} \hspace{\stretch{1}}(4.49)

Should we wish a more exact plot for any specific values of the viscosities, we could plot exactly with software the vector field described by these velocities.


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: