$$ \newcommand{\D}{\displaystyle} \renewcommand{\eqref}[1]{Eq.~(\ref{#1})} $$




3.1.1 Example: Couette-Poiseuille flow

In a Couette-Poiseuille flow, we consider fluid flow constrained between two walls, where the upper wall is moving at a prescribed velocity \( U_0 \) and at a distance \( Y=L \) from the bottom wall. Additionally, the flow is subjected to a prescribed pressure gradient \( \frac{\partial p}{\partial x} \).

Figure 32: An illustration of Couette flow driven by a pressure gradient and a moving upper wall.

We assume that the vertical velocity component \( V=0 \) and consequently we get from continuity: \( \frac{\partial U}{\partial X}=0 \) which implies that \( U=U(Y) \), i.e. the velocity in the streamwise X-direction depends on the cross-wise Y-direction only.

The equation of motion in the Y-direction simplifies to \( \frac{\partial p}{\partial Y}=-\rho g \), whereas the equation of motion in the X-direction has the form: $$ \begin{equation*} \rho U \cdot \frac{\partial U}{\partial X}=-\frac{\partial p}{\partial X} + \mu \left(\frac{\partial^2U}{\partial X^2}+ \frac{\partial^2U}{\partial Y^2}\right) \end{equation*} $$ which due to the above assumptions and ramifications reduces to $$ \begin{equation} \frac{d^2U}{dY^2} =\frac{1}{\mu}\frac{dp}{dX} \tag{3.20} \end{equation} $$

with no-slip boundary conditions: \( U(0)=0,\ U(L)=U_0 \). To render equation (3.20) on a more appropriate and generic form we introduce dimensionless variables: \( u=\frac{U}{U_0},\ y=\frac{Y}{L},\ P=-\frac{1}{U_0}(\frac{dp}{dX})\frac{L^2}{\mu} \), which yields: $$ \begin{equation} \frac{d^2u}{dy^2}=-P \tag{3.21} \end{equation} $$ with corresponding boundary conditions: $$ \begin{equation} u=0 \, \text{for} \, y=0, \qquad u=1 \, \text{for} \, y=1 \tag{3.22} \end{equation} $$

An analytical solution of equation (3.21) with the corresponding boundary conditions (3.22) may be found to be: $$ \begin{equation} u=y \cdot \left[ 1+\frac{P}{2}(1-y) \right] \tag{3.23} \end{equation} $$

Observe that for \( P\leq -2 \) we will get negative velocities for some values of \( y \). In Figure 33 velocity profiles are illustrated for a range of non-dimensional pressure gradients P.

Figure 33: Velocity profiles for Couette-Poiseuille flow with various pressure gradients.

To solve (3.21) numerically, we represent is as a system of equations: $$ \begin{equation} \left.\begin{matrix} &u'(y)=u_1(y) \\ &u_1'(y)=-P \end{matrix}\right.\ \tag{3.24} \end{equation} $$

with corresponding boundary conditions: $$ \begin{equation} u(0)=0,\ u(1)=1 \tag{3.25} \end{equation} $$

To solve this boundary value problem with a shooting method, we must find \( s=u'(0)=u_1(0) \) such that the boundary condition \( u(1)=1 \) is satisfied. We can express this condition by: $$ \begin{align*} \phi(s)=u(1;s)-1, \qquad \text{ such that } \qquad \phi(s)=0 & \qquad \text{when} \qquad s=s^* \end{align*} $$

We guess two values \( s^0 \) and \( s^1 \) and compute the correct \( s \) by linear interpolation due to the linearity of system of ODEs (3.24). For the linear interpolation see equation (3.14).

The shooting method is implemented in the python-code Poiseuille_shoot.py and results are computed and plotted for a range of non-dimensional pressure gradients and along with the analytical solution.

# src-ch2/Couette_Poiseuille_shoot.py;ODEschemes.py @ git@lrhgit/tkt4140/src/src-ch2/ODEschemes.py;

from ODEschemes import euler, heun, rk4
import numpy as np
from matplotlib.pyplot import *

# change some default values to make plots more readable 
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT

L = 1.0
y = np.linspace(0,L,N+1)

def f(z, t):
    """RHS for Couette-Posieulle flow"""
    zout = np.zeros_like(z)
    zout[:] = [z[1], -dpdx]
    return zout 

def u_a(y,dpdx):
    return y*(1.0 + dpdx*(1.0-y)/2.0);

beta=1.0 # Boundary value at y = L

# Guessed values
s=[1.0, 1.5]


dpdx_list=[-5.0, -2.5, -1.0, 0.0, 1.0,2.5, 5.0]

for dpdx in dpdx_list:
    phi = []
    for svalue in s:
        z0[1] = svalue
        z = rk4(f, z0, y)
        phi.append(z[-1,0] - beta)
    # Compute correct initial guess 
    s_star = (s[0]*phi[1]-s[1]*phi[0])/(phi[1]-phi[0])
    z0[1] = s_star
    # Solve the initial value problem which is a solution to the boundary value problem
    z = rk4(f, z0, y)

    legends.append('rk4: dp='+str(dpdx))
    # Plot the analytical solution
    plot(u_a(y, dpdx),y,':')
    legends.append('exa: dp='+str(dpdx))

# Add the labels
legend(legends,loc='best',frameon=False) # Add the legends