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

# 7.2 Confined, unsteady Couette flow

The classical version of unsteady Couette flow with $$b=\infty$$ has been presented as Stokes first problem and discussed in section (3.3 Notes on similarity solutions).

Figure 89: Confined, unsteady Couette flow or channel flow. The channel with is $$b$$ In this section we will look at the problem of unsteady Couette flow, confined by two walls, which has the following governing equation: $$\begin{equation} \tag{7.6} \frac{\partial U}{\partial \tau}=\nu \frac{\partial^2U}{\partial Y^2},\ 0 < Y < b \end{equation}$$

with the following boundary conditions, representing the presence of the two walls or channel if you like: $$\begin{equation} \tag{7.7} \left.\begin{matrix} U(0,\ \tau)=&U_0\\ U(b,\ \tau)=&0 \end{matrix}\right\}=\tau \geq 0 \end{equation}$$

Further, the parabolic problem also needs initial conditions to be solved and we assume: $$\begin{equation} \tag{7.8} U(Y,\ \tau)=0,\ \tau < 0 \end{equation}$$

In the section 3.3 Notes on similarity solutions we have presented several ways to render (7.6)) dimensionless, and for the current problem we introduce the following dimensionless variables: $$\begin{equation} \tag{7.9} y=\frac{Y}{b},\ u=\frac{U}{U_0},\ t= \frac{\tau\nu}{b^2} \end{equation}$$

which allow (7.6) to be written: $$\begin{equation} \tag{7.10} \frac{\partial u}{\partial t}=\frac{\partial^2u}{\partial y^2},\ 0 < y < 1 \end{equation}$$

with the corresponding boundary conditions: $$\begin{equation} \tag{7.11} \left.\begin{matrix} u(0,\ t)=&1\\ u(1,\ t)=&0 \end{matrix}\right\},\ t \geq 0 \end{equation}$$

and initial conditions: $$\begin{equation} \tag{7.12} u(y,\ t)=0,\ t < 0 \end{equation}$$

As for the problem in the section 3.3 Notes on similarity solutions, this present example may also be formulated as a heat conduction problem. The problem of confined, unsteady Couette flow has the following analytical solution: $$\begin{equation} \tag{7.13} u(y,t)=1-y-\frac{2}{\pi}\cdot \sum^{\infty}_{n=1}\frac{1}{n}\exp[-(n\pi)^2t]\sin(n\pi y) \end{equation}$$

A detailed derivation of (7.13) may be found in appendix G.6 of Numeriske beregninger.

We discretize (7.10) by a forward difference for the time $$t$$: $$\begin{equation} \tag{7.14} \frac{\partial u}{\partial t}\bigg|^n_j\approx \frac{u^{n+1}_j-u_j^n}{\Delta t} \end{equation}$$

and central differences for the spatial coordinate $$y$$: $$\begin{equation} \tag{7.15} \frac{\partial^2u}{\partial y^2} \approx \frac{u^n_{j+1}-2u^n_j+u^n_{j-1}}{(\Delta y)^2} \end{equation}$$ where: $$\begin{equation*} t_n=n\cdot \Delta t,\ n=0,1,2,\dots,\ y_j=j\cdot \Delta y,\ j=0,1,2,\dots \end{equation*}$$

Substitution of (7.14) in (7.10) results in the following difference equation: $$\begin{equation} \tag{7.16} u_j^{n+1}= D\, (u_{j+1}^n+u^n_{j-1})+(1-2D) \, u_j^n \end{equation}$$

where $$D$$ is a dimensionless group, commonly denoted the diffusion number or the Fourier number in heat conduction. See section 13.1 in )) for a discussion of (7.16). $$\begin{equation} \tag{7.17} D=\frac{\Delta t}{(\Delta y)^2}=\nu \frac{\Delta \tau}{(\Delta Y)^2} \end{equation}$$

A scheme with Forward differences for the Time (evolutionary) variable and Central differences for the Space variable, is commonly referred to as a FTCS (Forward Time Central Space). For a FTCS-scheme we normally mean a scheme which is first order in $$t$$ and second order in $$y$$. Another common name for this scheme is the Euler-scheme.

In the section 7.6 Truncation error, consistency and convergence we show that for $$D=1/6$$, the scheme (7.16)) is second order in $$t$$ and fourth order in $$y$$. Further, in fluid mechanics it is customary to write $$u^n_j$$ rather than $$u_{j,n}$$, such that index for the evolutionary variable has a super index. We seek to adopt this convention in general.

Figure 90: Numerical stencil of the FTCS scheme. In Figure 90 we try to illustrate that the scheme is explicit, meaning that the unknown value at time $$n+1$$ can be found explicitly from the formula without having to solve an equation system. In other words, the unknown value at time $$n+1$$ is not implicitly dependent on other values at other spatial locations at time $$n+1$$.

The above example is implemented in the code below. Download the code and experiment using different diffusion numbers. The FTCS-scheme is explicit and thus has a stability constraint. We will look further into stability in the next sections, but as we will see in the animations displayed below, the stability limit in this example is $$D = \frac{1}{2}$$.

# src-ch5/couette_FTCS.py;Visualization.py @ git@lrhgit/tkt4140/src/src-ch5/Visualization.py;

import matplotlib; matplotlib.use('Qt5Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()

import numpy as np
from math import exp, sin, pi

def analyticSolution(y, t, N=100):

""" Method that calculates the analytical solution to the differential equation:
du/dt = d^2(u)/dx^2 , u = u(y,t), 0 < y < 1
Boundary conditions: u(0, t) = 1, u(1, t) = 0
Initial condition: u(t, 0) = 0 t<0,  u(t, 0) = 1 t>0

Args:
t(float): time
N(int): truncation integer. Truncate sumation after N elements

Returns:
w(float): velocity, us - ur
"""
sumValue = 0
for n in range(1,N+1):
temp = np.exp(-t*(n*np.pi)**2)*np.sin(n*np.pi*y)/n
sumValue += temp
u = 1 - y - (2/pi)*sumValue
return u

def solveNextTimestepFTCS(Uold, D, U_b=1, U_t=0):
""" Method that solves the transient couetteflow using the FTCS-scheme..
At time t=t0 the plate starts moving at y=0
The method solves only for the next time-step.
The Governing equation is:

du/dt = d^2(u)/dx^2 , u = u(y,t), 0 < y < 1

Boundary conditions: u(0, t) = 1, u(1, t) = 0

Initial condition: u(t, 0) = 0 t<0,  u(t, 0) = 1 t>0

Args:
uold(array): solution from previous iteration
D(float): Numerical diffusion number

Returns:
unew(array): solution at time t^n+1
"""
Unew = np.zeros_like(Uold)

Uold_plus = Uold[2:]
Uold_minus = Uold[:-2]
Uold_mid = Uold[1:-1]

Unew[1:-1] = D*(Uold_plus + Uold_minus) + (1 - 2*D)*Uold_mid
Unew = U_b
Unew[-1] = U_t

return Unew

if __name__ == '__main__':

import numpy as np
from Visualization import createAnimation

D = 0.505 # numerical diffusion number

N = 20
y = np.linspace(0, 1, N + 1)
h = y - y
dt = D*h**2
T = 0.4 # simulation time
time = np.arange(0, T + dt, dt)

# Spatial BC
U_bottom = 1.0 # Must be 1 for analytical solution
U_top = 0.0 # Must be 0 for analytical solution

# solution matrices:
U = np.zeros((len(time), N + 1))
U[0, 0] = U_bottom   # no slip condition at the plate boundary
U[0,-1] = U_top
#
Uanalytic = np.zeros((len(time), N + 1))
Uanalytic[0, 0] = U[0,0]

for n, t in enumerate(time[1:]):

Uold = U[n, :]

U[n + 1, :] = solveNextTimestepFTCS(Uold, D, U_b=U_bottom, U_t=U_top)

Uanalytic[n + 1, :] = analyticSolution(y,t)

U_Visualization = np.zeros((1, len(time), N + 1))
U_Visualization[0, :, :] = U

createAnimation(U_Visualization, Uanalytic, ["FTCS"], y, time, symmetric=False)


In the following two animations we show numerical solutions obtained with the explicit FTCS scheme, as well as with two implicit schemes (Crank-Nicolson and Laasonen schemes), along with the analytical solution. The first animation shows results for $$D=0.5$$, for which the FTCS scheme is stable. In the second animation we use $$D=0.504$$, a value for which the FTCS scheme's stability limit is exceeded. This observation is confirmed by the oscillatory character of the numerical solution delivered by this scheme. .

Animation of numerical results obtained using three numerical schemes as well as the analytical solution using $$D=0.5$$.

Animation of numerical results obtained using three numerical schemes as well as the analytical solution using $$D=0.504$$. The FTCS displays an unstable solution.