Crank-Nicolson time scheme

Properties

• Implicit
• Second order in time (when the off-centring coefficient is 1, after startup)
• Transient
• Off-centring often required for stability on complex flows
• Boundedness not guaranteed

Discretisation

Consider a generic transport equation for a scalar $\phi$:

$$\ddt{\phi} + \div \left(\vec U \phi \right) - \laplacian \left(\Gamma \phi \right) = S$$

Crank-Nicolson evaluates terms at the mid-point of the time step. OpenFOAM achieves this by splitting each contribution into an implicit part that enters the linear system (fvm:: terms) and an explicit part that is carried forward from previous time levels. For the temporal derivative alone, the scheme stores an auxiliary field $\mathrm{ddt0}(\phi)$ on the mesh registry.

On a static mesh, with uniform time step $\Delta t$ and off-centring coefficient $\alpha \in [0, 1]$, the temporal rate after startup is represented as:

$$\ddt{\phi} \approx \frac{1 + \alpha}{\Delta t}\left(\phi - \old{\phi}\right) - \mathrm{ddt0}(\phi)$$

The field $\mathrm{ddt0}(\phi)$ is updated each time step from the previous time levels:

$$\mathrm{ddt0}(\phi) = \frac{1 + \alpha}{\Delta t} \left(\old{\phi} - \oldold{\phi}\right) - \alpha\, \mathrm{ddt0}_{o}(\phi)$$

When $\alpha = 1$ the scheme is fully centred; when $\alpha = 0$ it reduces to the Euler implicit scheme.

Implicit and explicit contributions

When fvm::ddt(phi) is used, the temporal term contributes to the linear system for $\phi$ as:

Contribution	Discrete form
Implicit (matrix diagonal)	$\displaystyle\frac{1 + \alpha}{\Delta t}\, V\,\phi$
Explicit (source)	$\displaystyle\left[\frac{1 + \alpha}{\Delta t}\,\old{\phi} + \mathrm{ddt0}(\phi)\right] V$

where $V$ is the cell volume. The corresponding explicit evaluation from fvc::ddt(phi) is the difference between these two parts.

The implicit weighting of off-centred terms is related to the off-centring coefficient by:

$$c_{\mathrm{CN}} = \frac{1}{1 + \alpha}$$

On the first time step the scheme uses the Euler implicit coefficient $1/\Delta t$ to avoid start-up problems; from the second step onward the $(1 + \alpha)/\Delta t$ formulation applies.

Usage

The scheme is specified using:

ddtSchemes
{
    default         CrankNicolson <coeff>
    ddt(phi)        CrankNicolson <coeff>;
}

The mandatory <coeff> is the off-centring coefficient $\alpha$ in the range $[0, 1]$:

• 0: equivalent to Euler implicit
• 1: fully centred Crank-Nicolson
• Values such as 0.9 are commonly used as a compromise between accuracy and robustness

The coefficient may also be specified as a Function1 to ramp from Euler to Crank-Nicolson over an initial period, for example:

ddtSchemes
{
    default         CrankNicolson
    ocCoeff
    {
        type        scale;
        scale       linearRamp;
        duration    0.01;
        value       0.9;
    };
}

Compared with the backward scheme, Crank-Nicolson requires only the current and old-time field values (not old-old-time values in the matrix assembly), but is generally less robust. The backward scheme provides formal second-order accuracy with stronger stability characteristics.

Further information

Source code

• Foam::fv::CrankNicolsonDdtScheme

Reference

• Crank and Nicolson 11