Simulating a Fluid with Chorin’s Method

This article describes the numerical method implemented in this project: the equations, Chorin’s splitting scheme, the staggered grid, the discretisation of each term, the boundary conditions, the solution of the pressure Poisson equation, and the choice of the time step. The discretisation follows the book by Griebel, Dornseifer, Neunhoeffer, Numerical Simulation in Fluid Dynamics (hereafter “Griebel”); deviations from it are noted explicitly.

1. Equations

The flow of a viscous incompressible fluid is described by the Navier–Stokes equations:

where u = (u, v) is the velocity, p the pressure, ρ the density, and ν the kinematic viscosity. There are no body forces. Bold u denotes the velocity vector, plain u its x-component.

Differential operators. In two-dimensional Cartesian coordinates the “nabla” operator is the vector of partial derivatives

and all operators in the equations are built from it.

• Divergence ∇·u — the dot product of ∇ with a vector; the result is a scalar:

It is the relative rate of change of the volume of a fluid parcel (the balance of inflow and outflow through a small volume). The condition ∇·u = 0 expresses incompressibility: as much fluid flows into any volume as flows out of it.

• Gradient ∇p — nabla applied to a scalar; the result is a vector:

It points in the direction of the fastest pressure increase, so the force −∇p (per unit volume) pushes the fluid from high-pressure regions towards low-pressure ones.

• Convective derivative (u·∇)u. First, u and ∇ are combined into a scalar operator — the derivative along the flow:

which is then applied to each velocity component; the result is a vector:

Together with ∂u/∂t it gives the acceleration of a fluid parcel, Du/Dt = ∂u/∂t + (u·∇)u: the parcel’s velocity changes not only with time at a fixed point, but also because the parcel is carried to a point where the velocity is different. This is the only nonlinear term in the equations.

• Laplacian Δ = ∇·∇ — the divergence of the gradient, a second-order scalar operator:

It is applied to a vector component-wise:

The Laplacian is proportional to the deviation of the value at a point from the average over its neighbourhood, so the term νΔu describes viscous diffusion: friction evens out the velocity of neighbouring fluid layers.

Component form. Substituting these expressions gives a system of three scalar equations for the three unknowns u, v, p:

The convective term (u·∇)u can be written in two forms. For the u component, expand the derivatives of the products:

Since ∇·u = 0, the second term vanishes, and the advective form (the first brace) coincides with the conservative (divergence) form, which is the one used in the code:

or, in vector notation, (u·∇)u = ∇·(u ⊗ u).

A distinctive feature of the incompressible formulation: pressure has no evolution equation of its own. It acts as a Lagrange multiplier that at every instant “adjusts itself” so that the velocity field stays divergence-free. Chorin’s method exploits exactly this property.

2. Chorin’s Scheme (Projection Method)

2.1. Idea

The method was proposed by A. Chorin (1968) and independently by R. Temam (1969). The time step tⁿ → tⁿ⁺¹ = tⁿ + Δt is split into two parts:

1. Predictor: move the fluid without pressure. This gives an intermediate field u* which in general violates incompressibility: ∇·u* ≠ 0.
2. Corrector: pressure fixes the velocity. Choose the pressure p so that its gradient corrects u* exactly as much as needed to satisfy the incompressibility equation:

Where fluid has “accumulated”, the pressure comes out high. Its gradient pushes the excess to the neighbours, and in the end ∇·uⁿ⁺¹ = 0.

Why such a correction always exists and is unique. This follows from the Helmholtz–Hodge decomposition (statement and proof — [7, § 1.3]; survey — [8]). Any vector field can be split, and in a unique way, into two parts:

The first part describes how the fluid swirls and flows without being compressed. The second describes only the “excess” compression and expansion. The corrector simply discards the second part.

2.2. Three Stages

Stage 1. Predictor — an explicit Euler step for convection and diffusion, without pressure:

Stage 2. Pressure Poisson equation. We seek the pressure pⁿ⁺¹ whose gradient will make the velocity divergence-free in stage 3. Stage 3 has not been performed yet, but we already know its formula. So we write two requirements on uⁿ⁺¹:

Derivation. Take the divergence of the first equality and substitute the second:

Since ∇·∇p = Δp (the divergence of the gradient is the Laplacian), we obtain

The velocity uⁿ⁺¹ has dropped out of the equation. What remains is a single scalar equation for a single unknown function p, with a known right-hand side: u* was already computed in stage 1.

Stage 3. Corrector:

after which the boundary values are set again.

3. Staggered (MAC) Grid

The domain is divided into n_x × n_y cells with spacings Δx = L/n_x, Δy = H/n_y. The unknowns are stored at different points of a cell:

• p — at the cell centre;
• u — at the centres of the vertical (left/right) faces;
• v — at the centres of the horizontal (bottom/top) faces.

v(j+1,i)
     +-------o-------+
     |               |
u(j,i) o     * p     o u(j,i+1)
     |               |
     +-------o-------+
          v(j,i)

This arrangement (“marker-and-cell”, Harlow & Welch, 1965) has three important properties:

1. The pressure difference between neighbouring cells falls exactly on the face where the corresponding velocity component lives — the pressure gradient is computed without interpolation.
2. The velocity divergence in a cell is the difference of fluxes through its faces, again without interpolation.
3. There are no “checkerboard” pressure oscillations, which are typical of collocated grids.

4. Spatial Discretisation

Below we use “physical” half-integer indices: u_(i+1/2, j) is the velocity on the right face of cell (i, j), v_(i, j+1/2) on the top face, p_(i,j) at the centre.

4.1. Diffusion

Five-point Laplacian, second-order accurate:

and similarly for v.

4.2. Convection: a Blend of Central and Donor-Cell Schemes

Each convective derivative is written as a difference of fluxes through the faces of a control volume around the velocity node, for example

The flux through a face depends on the transporting velocity a on that face and on the values of the transported quantity φ_L, φ_R on either side of it:

• γ = 0 — central scheme (second order, but unstable/oscillatory at large grid Reynolds numbers);
• γ = 1 — donor-cell scheme (upwind, first order): for a > 0 the flux is aφ_L, for a < 0 it is aφ_R, i.e. the quantity is taken “from the donor cell” upstream;
• intermediate γ give a linear blend. The upwind contribution is equivalent to an artificial viscosity ~γ|a|Δx/2.

The transporting velocities are the arithmetic mean of the two nearest nodes. For the u equation:

and symmetrically for the v equation. This is exactly Griebel’s scheme (formulas 3.19–3.20 of the book), merely written in terms of a flux function.

The conservative form guarantees that the flux leaving one control volume enters the neighbouring one — momentum is neither created nor lost by the scheme inside the domain.

Get Simon’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

Why the upwind contribution is needed. A pure central scheme produces unphysical sawtooth oscillations at high Reynolds numbers.

4.3. Divergence and Gradient

Divergence of u* at the cell centre:

Pressure gradient on the faces:

5. Boundary Conditions

This section describes boundary conditions on a MAC grid in general form: walls, inflow, outflow, symmetry, periodicity.

5.1. How the Boundary Works on a MAC Grid

The domain boundary runs along cell faces. On each side one velocity component is normal to the boundary and the other is tangential, and their nodes are placed differently:

• normal component — the nodes lie exactly on the boundary, and the value is set directly;
• tangential component — there are no nodes on the boundary. The nearest interior node is half a step away from the boundary, and beyond the boundary lies a ghost node. The condition is imposed through the ghost value.

tangential (u at the bottom):          normal (v at the bottom):

 y = +Δy/2  ── o  u_in  (interior)      y = +Δy  ── o  v_in  (interior)
 y = 0   ══════════════  boundary       y = 0   ═══ o ═══  v_b  (on boundary)
 y = −Δy/2  ── o  u_gh  (ghost)

Any velocity condition reduces to two elementary techniques.

Prescribed value (Dirichlet), φ = φ_b on the boundary. For the normal component simply φ_(on boundary) = φ_b. For the tangential component the boundary value is obtained by linear interpolation between the ghost and interior nodes, so we require their mean to equal φ_b:

Simply writing φ_gh = φ_b is wrong: the boundary value would then be (φ_b + φ_in)/2 ≠ φ_b.

Prescribed normal derivative (Neumann), ∂φ/∂n = 0. The value on the far side of the boundary is copied from the near side: φ_gh = φ_in for the tangential component and φ_(on boundary) = φ_in for the normal one. A nonzero derivative g gives a shift: φ_gh = φ_in + g·h, where h is the spacing across the boundary.

5.2. Types of Velocity Conditions

Denote the normal component by u_n and the tangential one by u_t. The standard set:

Notes:

• No-slip and moving wall. A viscous fluid neither slides along a wall nor passes through it. For a stationary wall the ghost value is the mirror image with a minus sign: the straight line through u_gh and u_in crosses zero exactly at the wall.
• Free-slip. No friction against the wall, only impermeability. The same formulas describe an axis or plane of symmetry: there is neither flux nor shear across it.
• Inflow. The profile U_in is specified by the user: constant, parabolic, etc.
• Outflow. The flow “leaves as it is”: normal derivatives are set to zero. This is an approximate condition. The outflow boundary is placed far from the region of interest so that it does not distort the solution.
• Periodicity. The domain is “glued”: the right boundary is the same face as the left one. The normal-component nodes on the boundary become ordinary unknowns (shared by both sides), and the ghost columns of the tangential component are filled with values from the other side: for periodicity in x, v(:,1) = v(:,end-1), v(:,end) = v(:,2).

The same rules apply on the faces of internal obstacles: obstacle cells are excluded from the computation, and their faces are walls.

5.3. Pressure

Pressure has no physical boundary condition of its own. It follows from what the corrector uⁿ⁺¹ = u* − (Δt/ρ)∇p does to the normal velocity on the boundary. There are two options.

Neumann, ∂p/∂n = 0. Used where the normal velocity on the boundary is prescribed: walls, inflow, free-slip, and also outflow with a copied u_n. The normal velocity already has the required value and must not change during the correction. From u_nⁿ⁺¹ = u_n* − (Δt/ρ)·∂p/∂n and u_nⁿ⁺¹ = u_n* = u_(n,b) it follows that

Discretely: the divergence of a boundary cell includes the prescribed flux through the boundary, and the pressure gradient on the boundary face is not computed. This is equivalent to a ghost pressure p₀ = p₁, for which the term (p₀ − p₁)/h² in the Laplacian vanishes. Hence the diagonal of the one-dimensional matrix has −1 instead of −2 at such an end.

This condition is not physical: the true normal pressure gradient at a wall equals ρν·∂²u_n/∂n² and is nonzero. This produces a numerical boundary layer in the pressure error of thickness O(√(νΔt)).

Dirichlet, p = p_b. Used where the pressure is prescribed and the normal velocity is not known in advance: outflow to the atmosphere, flow driven by a pressure difference. On such a boundary the normal velocity is not imposed but found by the corrector. Via the ghost cell p₀ = 2p_b − p₁:

i.e. −3 on the diagonal, and the known term 2p_b/h² is moved to the right-hand side. The gradient on the boundary face, (p₁ − p_b)/(h/2), enters the corrector.

Periodicity. The matrix “wraps around”: the first and last cells in the periodic direction become neighbours. The one-dimensional matrix is circulant, with −2 along the whole diagonal and ones in the corners.

If at least one part of the boundary carries a Dirichlet condition, the pressure level is determined and the matrix is nonsingular. If the conditions are Neumann or periodic everywhere, the pressure is determined only up to a constant, and its level must be fixed.

6. Solving the Linear System for Pressure

The Poisson equation from stage 2, written in every grid cell, is a system of linear equations Ap = b with N = n_x·n_y unknowns. The system has the following properties:

• Sparsity. Each row has at most five nonzero entries: a cell is coupled only to its four neighbours (five-point Laplacian, section 4.1).
• Constant matrix, changing right-hand side. A depends only on the grid, while b = (ρ/Δt)∇·u* changes every step. So the system has to be solved anew at every time step, and this is usually the most expensive part of the computation.
• Singularity. With Neumann conditions on all walls (section 5.3) the pressure is determined only up to a constant, and the matrix is singular. The pressure level has to be fixed, but only its gradient is physically meaningful.

Such systems are solved by iterative or direct methods. Both approaches are discussed in detail, using the heat equation as an example, in articles [9], [10]. The project uses successive over-relaxation for the following reasons:

• Memory. Iterative methods store only the pressure array itself (plus an array of old values for the Jacobi method) — O(N). The matrix need not be assembled: the formula “cell value in terms of its neighbours” is enough. A direct method stores the factors L, U, and elimination fills in the zeros inside the matrix band: in 2D this means O(N log N) nonzeros, in 3D O(N^(4/3)), and the factorisation costs O(N²) operations.
• Simplicity. The method hardly changes when going from 2D to 3D or to domains with obstacles: only the cell’s neighbours and the boundary formulas change.
• Accuracy. There is no point in solving the system more accurately than the scheme error ~h², and iterations can be stopped at the required tolerance. A good initial guess from the previous step makes this cheap.

The price is slow convergence: the number of iterations grows as O(N_x) even with the optimal ω, and the convergence rate is very sensitive to the choice of ω. In addition, after the corrector the velocity is divergence-free only up to the tolerance ε, not up to round-off. For more on SOR, the choice of ω, and the stopping criterion, see the articles “Solving the Steady-State Heat Equation” [9, § 3.3], [10, § 3.3].

6.1. Successive Over-Relaxation (SOR)

The discrete Poisson equation in cell (i, j) is the five-point Laplacian from section 4.1 with right-hand side b_(i,j) = (ρ/Δt)(∇·u*)_(i,j). Solve it for the central value:

The cells are visited one after another, and each immediately uses the already updated values of its neighbours (the Gauss–Seidel method). The correction is then “amplified” with the relaxation parameter ω:

All values are overwritten “in place” in a single array.

Specifics for pressure:

• Boundary condition ∂p/∂n = 0 is imposed through a ghost layer: before each iteration the values in the wall-adjacent cells are copied into the ghost cells (section 5.3).
• Initial guess — the pressure from the previous time step. The pressure field changes little from step to step, so few iterations are needed.
• Stopping criterion — the relative residual ‖b − Ap‖ / ‖b‖ < ε.
• Pressure level. With Neumann conditions on all walls the pressure is determined only up to a constant. SOR converges on such a system too, provided the right-hand side sums to zero. Therefore its mean is subtracted from b (a correction at round-off level), and the mean of the computed p is subtracted from p. This does not affect the velocity: a constant shift does not change the gradient.

6.2. Alternative: Direct Methods

Direct methods for this system are described in the same articles [9, § 3.2], [10, § 3.2]: sparse LU factorisation, Cholesky factorisation, block tridiagonal (Thomas) algorithm, fast direct solvers.

For a two-dimensional cavity of moderate size, sparse LU with a single factorisation is faster than SOR. Iterative methods become preferable as the grid grows and in 3D, where fill-in of the factors eats up memory.

7. Time Step

The predictor is explicit, so the step is limited by two stability conditions.

Diffusion (explicit scheme for the heat equation on a five-point stencil):

Convection (the Courant–Friedrichs–Lewy condition in two-dimensional form):

The time step is therefore:

where cfl < 1.

When estimating |u|_max, |v|_max, ghost values are excluded: in the ghost row above the lid, u = 2U_lid − u_inner can reach almost 2U_lid, which would halve the step for no physical reason.

Donor-cell weight. At the same time we compute

which corresponds to Griebel’s stability condition for the hybrid scheme, γ ≥ max(|u|Δt/Δx, |v|Δt/Δy). The smallest admissible value is taken, so the scheme stays as close as possible to the central one (less dissipative). The upper bound of one is a safeguard: when the CFL condition holds, it never activates.

The MATLAB implementation of the lid-driven cavity flow problem can be found here: https://github.com/Semen1996/incompressible-flow-chorin-method

References

1. Griebel M., Dornseifer T., Neunhoeffer T. Numerical Simulation in Fluid Dynamics: A Practical Introduction. SIAM, 1998.
2. Chorin A. J. Numerical solution of the Navier–Stokes equations. Mathematics of Computation, 22 (1968), 745–762.
3. Temam R. Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires. Arch. Rational Mech. Anal., 33 (1969), 377–385.
4. Harlow F. H., Welch J. E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics of Fluids, 8 (1965), 2182–2189.
5. Ghia U., Ghia K. N., Shin C. T. High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method. J. Comput. Phys., 48 (1982), 387–411.
6. Guermond J.-L., Minev P., Shen J. An overview of projection methods for incompressible flows. Comput. Methods Appl. Mech. Engrg., 195 (2006), 6011–6045.
7. Chorin A. J., Marsden J. E. A Mathematical Introduction to Fluid Mechanics. 3rd ed., Springer, 1993. — § 1.3: the Helmholtz–Hodge decomposition and its connection to pressure and the projection method.
8. Bhatia H., Norgard G., Pascucci V., Bremer P.-T. The Helmholtz–Hodge decomposition — a survey. IEEE Trans. Visualization and Computer Graphics, 19(8) (2013), 1386–1404.
9. Solving the Steady-State Heat Equation. Part 1. The One-Dimensional Case. https://github.com/Semen1996/steady-state-heat-1d
10. Solving the Steady-State Heat Equation. Part 2. The Two-Dimensional Case. https://github.com/Semen1996/steady-state-heat-2d