Step 1: Navier-Stokes solver for Lid-driven flow in rectangular cavity¶

This iPython notebook shows how to solve the lid-driven cavity problem with a Navier-Stokes solver. The main features are summarized as follows:

• Incompressible Navier-Stokes equations in two dimensions with primitive variables p, u and v and constant material properties

• Dirichlet boundary conditions for the velocities and Neumann boundary conditions for the pressure, except for a fixed point in the lower left (SW) corner

• Explicit projection method to decouple pressure and velocity

• Finite Difference discretization on a uniform orthogonal cartesian grid, first order explicit Euler time-stepping

• Pressure defined in cell centers and velocities defined in backward staggered grid

• Second derivatives discretized with centered differences, even for the non-linear terms in the momentum equation

• Viscous term in momentum equation is calculated implicitly

• Poisson equations for pressure and implicit viscous momentum term are solved with a direct solver or alternatively with an iterative solver and a multigrid preconditioner

Incompressible Navier-Stokes equations¶

The incompressible Navier-Stokes equations are given in vector notation with the velocity vector $u$, the time $t$, the density $\rho$, the pressure $p$, the kinematic viscosity $\nu$ and the Nabla operator $\nabla$.

$$\nabla\cdot\vec{u}=0$$ $$\frac{\partial \vec{u}}{\partial t}+(\vec{u}\cdot\nabla)\vec{u}=-\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{u}$$

Geometry and boundary conditions¶

The physical problem is a rectangular cavity with a moving lid as depicted in the figure below. The width and height of the cavity are $L_x=L_y=1\,\text{m}$. The top lid is a non-slip wall and it is moving with a velocity $u=1\,\frac{\text{m}}{\text{s}}$ in positive $x$-direction. All other boundarys are stationary non-slip walls with zero velocities. The pressure gradient is $\vec{n}\cdot\nabla p=0$ at all four boundarys as it is common for non-slip walls. Finally, the pressure is set to a constant value $p=0$ at the lower left corner to completely define the pressure boundary and to obtain a positive definite coefficient matrix.

Material properties, test cases and discretization¶

The density is set to a fixed value $\rho=1 \frac{\text{kg}}{\text{m}^3}$. The kinematic viscosity $\nu$ is varied to obtain different Reynolds numbers

$$Re=\frac{U L}{\nu},$$

where $U=u=1\,\frac{\text{m}}{\text{s}}$ is the characteristic velocity and $L=L_x=1\,\text{m}$ is the characteristic length.

Two different test cases and an appropriate discretization adjusted to the test case are defined as shown in the table below.

$Re$	$\nu$	$\Delta x$	$\Delta y$	grid	$\Delta t$	$t_{sim}$
20	0.05 Pa s	0.05 m	0.05 m	20x20	0.05 s	5 s
100	0.01 Pa s	0.02 m	0.02 m	50x50	0.02 s	20 s

The grid spacings $\Delta x$ and $\Delta y$ were adjusted to obtain a reasonable Peclet number for mass transport:

$$Pe=\frac{U \Delta x}{\nu} \le 2.$$

The time step $\Delta t$ was adjusted to obey the time step restriction for an explicit scheme declared by the Courant-Friedrichs-Levy number:

$$CFL = \frac{U \Delta t}{\Delta x} \le 1.$$

The simulation time $t_{sim}$ is the time after the simulation was found to be in a steady state.

Results¶

Pressure and velocity fields for Re = 20¶

First, a low Reynolds number $Re=20$ was investigated. The resulting pressure and velocity fields at a steady state after t = 5 s are illustrated in the figure below. To compare the results, a simulation of the same test case with a similar discretization was performed in OpenFOAM. The OpenFoam simulation data was imported and plotted with the Python script plotOpenFoamResults.py.

Pressure contours and velocity vectors of the lid-driven cavity with Re = 20 at steady state after time t = 5 s.

Python results	OpenFOAM results

Pressure and velocity fields for Re = 100¶

A higher Reynolds number $Re=100$ was investigated. The resulting pressure and velocity fields at a steady state after t = 5 s are again illustrated in the figure below together with OpenFOAM results.

Pressure contours and velocity vectors of the lid-driven cavity with Re = 100 at steady state after time t = 20 s.

Python results	OpenFOAM results

Velocity profiles for Re = 100¶

The results of the Python CFD code were also compared to literature data by Ghia et al. [Ghia1982]. Firstly, the $x$-velocity $u$ was evaluated on a vertical profile going through the center of the cavity. The results are plotted in the figure below.

Vertical profile of the $x$-velocity $u$ for the lid-driven cavity with Re = 100.

Secondly, the $y$-velocity $v$ was evaluated on a horizontal profile going through the center of the cavity. The results are plotted in the figure below.

Horizontal profile of the $y$-velocity $v$ for the lid-driven cavity with Re = 100.

Code documentation¶

TBD

References¶

[Ghia1982] U Ghia, K.N Ghia, C.T Shin. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics. Volume 48, Issue 3, 1982, Pages 387-411. https://doi.org/10.1016/0021-9991(82)90058-4.