$$\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 $$$$\frac{d^2U}{dY^2} =\frac{1}{\mu}\frac{dp}{dX} \tag{3.20}$$$$

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: $$$$\frac{d^2u}{dy^2}=-P \tag{3.21}$$$$ with corresponding boundary conditions: $$$$u=0 \, \text{for} \, y=0, \qquad u=1 \, \text{for} \, y=1 \tag{3.22}$$$$

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

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: $$$$\left.\begin{matrix} &u'(y)=u_1(y) \\ &u_1'(y)=-P \end{matrix}\right.\ \tag{3.24}$$$$

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

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
LNWDT=5; FNT=11
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT

N=200
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]

z0=np.zeros(2)

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

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)

plot(z[:,0],y,'-.')
legends.append('rk4: dp='+str(dpdx))

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