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).
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.
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.