Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Channel flow with step pressure gradient

Posted by peeterjoot on April 1, 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 and statement.

This is problem 2.5 from [1].

Viscous fluid is at rest in a two-dimensional channel between stationary rigid walls with y = \pm h. For t \ge 0 a constant pressure gradient P = -dp/dx is imposed. Show that u(y, t) satifies

\begin{aligned}\frac{\partial {u}}{\partial {t}} = \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \frac{P}{\rho},\end{aligned} \hspace{\stretch{1}}(1.1)

and give suitable initial and boundary conditions. Find u(y, t) in form of a Fourier series, and show that the flow approximates to steady channel flow when t \gg h^2/\nu.

Solution

With only horizontal components to the flow, the Navier-Stokes equations for incompressible flow are

\begin{subequations}

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

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

\end{subequations}

Substitution of 2.2b into 2.2a gives us

\begin{aligned}\frac{\partial {u}}{\partial {t}} = -\frac{1}{{\rho}} \frac{\partial {p}}{\partial {x}} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.3)

Our equation to solve is therefore

\begin{aligned}\boxed{\frac{\partial {u}}{\partial {t}} = \Theta(t) \frac{P}{\rho} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.}\end{aligned} \hspace{\stretch{1}}(2.4)

This equation, for t < 0, allows for solutions

\begin{aligned}u = A y + B\end{aligned} \hspace{\stretch{1}}(2.5)

but the problem states that the fluid is at rest initially, so we don’t really have to solve anything (i.e. A = B = 0).

The no-slip conditions introduce boundary value conditions u(\pm h, t) = 0.

For t \ge 0 we have

\begin{aligned}\frac{\partial {u}}{\partial {t}} = \frac{P}{\rho} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.6)

If we attempt separation of variables with u(y, t) = Y(y) T(t), our equation takes the form

\begin{aligned}T' Y = \frac{P}{\rho} + \nu T Y''.\end{aligned} \hspace{\stretch{1}}(2.7)

We see that the non-homogenous term prevents successful application of separation of variables. Let’s modify our problem by attempting to recast our equation into a homogenous form by adding a particular solution for the steady state flow problem. That problem was the solution of

\begin{aligned}\frac{\partial^2 {{}}}{\partial {{y}}^2} u_s(y, 0) = -\frac{P}{\rho \nu}\end{aligned} \hspace{\stretch{1}}(2.8)

which has solution

\begin{aligned}u_s(y, 0) = \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right) + A y + B.\end{aligned} \hspace{\stretch{1}}(2.9)

The freedom to incoporate an h^2 constant into the equation as an integration constant has been employed, knowing that it will kill the y^2 contributions at y = \pm h to make the boundary condition matching easier. Our no-slip conditions give us

\begin{aligned}0 &= A h + B \\ 0 &= -A h + B.\end{aligned} \hspace{\stretch{1}}(2.10)

Adding this we have 2 B = 0, and subtracting gives us 2 A h = 0, so a specific solution that matches our required boundary value (and initial value) conditions is just the steady state channel flow solution we are familiar with

\begin{aligned}u_s(y, 0) = \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.12)

Let’s now assume that our general solution has the form

\begin{aligned}u(y, t) = u_H(y, t) + u_s(y, 0).\end{aligned} \hspace{\stretch{1}}(2.13)

Applying the Navier-Stokes equation to this gives us

\begin{aligned}\frac{\partial {u_H}}{\partial {t}} = \frac{P}{\rho} + \nu \frac{\partial^2 {{u_H}}}{\partial {{y}}^2} + \nu \frac{\partial^2 {{u_s}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.14)

But from 2.8, we see that all we have left is a homogenous problem in u_H

\begin{aligned}\frac{\partial {u_H}}{\partial {t}} = \nu \frac{\partial^2 {{u_H}}}{\partial {{y}}^2},\end{aligned} \hspace{\stretch{1}}(2.15)

where our boundary value conditions are now given by

\begin{aligned}0 &= u_H(\pm h, t) + u_s(\pm h) \\ &= u_H(\pm h, t),\end{aligned}

and

\begin{aligned}0 &= u(y, 0) \\ &= u_H(y, 0) + \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right),\end{aligned}

or
\begin{subequations}

\begin{aligned}u_H(\pm h, t) = 0\end{aligned} \hspace{\stretch{1}}(2.16a)

\begin{aligned}u_H(y, 0) = -\frac{P}{2 \rho \nu} \left( h^2 - y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.16b)

\end{subequations}

Now we can apply separation of variables with u_H = T(t) Y(y), yielding

\begin{aligned}T' Y = \nu T Y'',\end{aligned} \hspace{\stretch{1}}(2.17)

or

\begin{aligned}\frac{T'}{T} = \nu \frac{Y''}{Y} = \text{constant} = - \nu \alpha^2.\end{aligned} \hspace{\stretch{1}}(2.18)

Here a positive constant \nu \alpha^2 has been used assuming that we want a solution that is damped with time.

Our solutions are

\begin{aligned}T &\propto e^{- \nu \alpha^2 t} \\ Y &= A \sin \alpha y + B \cos\alpha y,\end{aligned} \hspace{\stretch{1}}(2.19)

or

\begin{aligned}u_H(y, t) = \sum_\alpha e^{-\alpha^2 \nu t} \left( A_\alpha \sin \alpha y + B_\alpha \cos\alpha y \right).\end{aligned} \hspace{\stretch{1}}(2.21)

We have constraints on \alpha due to our boundary value conditions. For our sin terms to be solutions we require

\begin{aligned}\sin (\alpha (\pm h)) = \sin n \pi\end{aligned} \hspace{\stretch{1}}(2.22)

and for our cosine terms to be solutions we require

\begin{aligned}\cos (\alpha (\pm h)) = \cos \left( \frac{\pi}{2} + n \pi \right)\end{aligned} \hspace{\stretch{1}}(2.23)

\begin{aligned}\alpha &= \frac{2 n \pi}{2 h} \\ \alpha &= \frac{2 n + 1 \pi}{2 h}\end{aligned} \hspace{\stretch{1}}(2.24)

respectively.

Our homogeneous solution therefore takes the form

\begin{aligned}u_H(y, t) = C_0 + \sum_{m > 0} C_m e^{ -(m \pi/2h)^2 \nu t } \left\{\begin{array}{l l}\sin \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{latex m$ even} \\ \cos \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{m odd} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.26)$

Our undetermined constants should be provided by the boundary value constraint at t = 0 2.16b, leaving us to solve the Fourier problem

\begin{aligned}-\frac{P}{2 \mu} \left( h^2 - y^2 \right)=\sum_{m \ge 0} C_m \left\{\begin{array}{l l}\sin \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{latex m$ even} \\ \cos \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{m odd} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.27)$

Multiplying by a sine and integrating will clearly give zero (even times odd function over a symmetric interval). Let’s see if there’s any scaling required to select out the C_m term

\begin{aligned}\int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) dy&=\frac{2h}{\pi} \int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) \pi dy/2h \\ &=\frac{2h}{\pi} \int_{-\pi/2}^{\pi/2}\cos m x \cos n x dx \\ &=\frac{h}{\pi} \int_{-\pi/2}^{\pi/2}\left( \cos( (m - n) \pi/2 ) +\cos( (m + n) \pi/2 ) \right) dx\end{aligned}

Note that since m and n must be odd, m \pm n = 2 c for some integer c, so this integral is zero unless m = n (consider m = 2 a + 1, n = 2 b + 1). For the m = n term we have

\begin{aligned}\int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) dy&=\frac{h}{\pi} \int_{-\pi/2}^{\pi/2}\left( 1+\cos( m \pi ) \right) dx \\ &=h\end{aligned}

Therefore, our constants C_m (for odd m) are given by

\begin{aligned}C_m &= -\frac{P h }{2 \mu} \int_{-h}^h\left( 1 - \left( \frac{y}{h}\right)^2 \right)\cos \left( \frac{ m \pi y }{2 h} \right) dy \\ &=-\frac{P h^2 }{2 \mu} \int_{-1}^1\left( 1 - x^2 \right)\cos \left( \frac{ m \pi x }{2} \right) x \\ \end{aligned}

With m = 2 n + 1, we have

\begin{aligned}C_{2 n + 1} = -\frac{16 P h^2 (-1)^n}{\mu \pi^3 (2 n + 1)^3}.\end{aligned} \hspace{\stretch{1}}(2.28)

For that calculation see: phy454/channelFlowWithStepPressureGradient.cdf.

Our complete solution is

\begin{aligned}u(y, t) =\frac{P h^2}{2 \mu} \left( 1 - \left( \frac{y}{h} \right)^2 \right)- \frac{16 P h^2 }{\mu \pi^3 }\sum_{n = 0}^\infty \frac{(-1)^n}{(2 n + 1)^3}e^{ -((2 n + 1) \pi/2h)^2 \nu t } \cos \left( \frac{ (2 n + 1) \pi y }{2 h} \right) .\end{aligned} \hspace{\stretch{1}}(2.29)

The largest of the damped exponentials above is the n = 0 term which is

\begin{aligned}e^{ - \pi^2 \nu t /h^2 },\end{aligned} \hspace{\stretch{1}}(2.30)

so if \nu t >> h^2 these terms all die off, leaving us with just the steady state.

Rather remarkably, this Fourier series is actually a very good fit even after only a single term. Using the viscosity and density of water, h = 1 \text{cm}, and P = 3 \times \mu_{\text{water}} \times (2 \text{cm}/{s})/ h^2 (parameterizing the pressure gradient by the average velocity it will induce), a plot of the parabola that we are fitting to and the difference of that from the first Fourier term is shown in figure (1).

Figure 1: Parabolic channel flow steady state, and difference from first Fourier term.

The higher order corrections are even smaller. Even the first order deviations from the parabola that we are fitting to is a correction on the scale of 1/100 of the height of the parabola. This is illustrated in figure (2) where the magnitude of the first 5 deviations from the steady state are plotted.

Figure 2: Difference from the steady state for the first five Fourier terms.

An animation of the time evolution above can be found in figure (3). If this animation is unavailable, it can also be found at http://youtu.be/0vZuv9HBtmo.

Figure 3: Time evolution of channel flow velocity profile after turning on a constant pressure gradient.

It’s also interesting to look at the very earliest part of the time evolution. Observe in figure (4)¬† (or http://youtu.be/dDkx8iLwOew) the oscillatory phenomina. Could some of that be due to not running with enough Fourier terms in this early part of the evolution when more terms are probably significant?

Figure 4: Early time evolution of channel flow velocity profile after turning on a constant pressure gradient.

References

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

Advertisements

Leave a Reply

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

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 )

Google+ photo

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

Connecting to %s

 
%d bloggers like this: