$$ \newcommand{\D}{\displaystyle} \renewcommand{\eqref}[1]{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 [7])) 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

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
            y(np.array): radial coordinat
            t(float): time
            N(int): truncation integer. Truncate sumation after N elements

            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
            uold(array): solution from previous iteration
            D(float): Numerical diffusion number
            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[0] = 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[1] - y[0]
    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.