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
#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:
y(np.array): radial coordinat
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[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.