A spectral deferred correction method for incompressible flow with variable viscosity

by   Jörg Stiller, et al.
TU Dresden

This paper presents a semi-implicit spectral deferred correction (SDC) method for incompressible Navier-Stokes problems with variable viscosity and time-dependent boundary conditions.The proposed method integrates elements of velocity- and pressure-correction schemes, which yields a simpler pressure handling in comparison to the SDPC method of Minion Saye (J. Comput. Phys. 375: 797-822, 2018). Combined with the discontinuous Galerkin spectral-element method for spatial discretization it can in theory reach arbitrary order of accuracy in time and space. Numerical experiments in three space dimensions demonstrate up to order 12 in time and 17 in space for constant, spatiotemporally varying as well as solution-dependent viscosity. The phenomenon of order reduction also reported by Minion Saye is observed in the case of time-dependent boundary conditions, where it manifests in terms of a slower convergence of the correction sweeps.



There are no comments yet.


page 1

page 2

page 3

page 4


Improved a priori error estimates for a discontinuous Galerkin pressure correction scheme for the Navier-Stokes equations

The pressure correction scheme is combined with interior penalty discont...

Assessment of high-order IMEX methods for incompressible flow

This paper investigates the competitiveness of semi-implicit Runge-Kutta...

A divergence-free finite element method for the Stokes problem with boundary correction

This paper constructs and analyzes a boundary correction finite element ...

A pressure-correction and bound-preserving discretization of the phase-field method for variable density two-phase flows

In this paper, we present an efficient numerical algorithm for solving t...

Generalized Internal Boundaries (GIB)

Representing large-scale motions and topological changes in the finite v...

A frequency-dependent p-adaptive technique for spectral methods

When using spectral methods, a question arises as how to determine the e...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

High-order discretization methods are gaining interest in fluid mechanics [76, 84, 15, 36], solid mechanics [30, 85], electrodynamics [22, 21]

and other areas of computational science that are governed by partial differential equations

[65]. This development is driven by the expectation of achieving a superior algorithmic efficiency which enables high-fidelity simulations at a scale beyond the reach of low-order methods. In simulations of processes evolving in space-time, the accuracy of spatial and temporal approximations must be tuned to each other. The natural and only scalable way for achieving this ist to match the convergence rates in space and time. However, examining recent work on high-order methods in computational fluid dynamics (CFD) reveals that the spatial order ranges typically from 4 to 16, whereas the temporal order rarely exceeds 3. This discrepancy constitutes the principal motivation of the present work which adopts the spectral deferred correction (SDC) method to reach arbitrary temporal convergence rates for incompressible flows with variable viscosity.

Before reviewing the state of the art in high-order time integration the implications of high-order space discretization shall be briefly recapitulated. In CFD, element based Galerkin methods with piecewise polynomial expansions represent the prevalent approach apart from spectral, finite-difference and isogeometric methods [32, 54, 16, 47]. The numerical properties of these methods have been thoroughly studied at the hand of convection, diffusion and wave problems as well as combinations thereof. Ainsworth and Wajid [2] and Gassner and Kopriva [40] showed that high-order spectral element approximations achieve far lower dispersion and dissipation errors than second-order finite-element or finite-volume methods using a comparable mesh spacing. This property yields a tremendous advantage in marginally or underresolved simulations of turbulent flows [10]. On the other hand, the condition of the discrete operators worsens when increasing the degree of the expansion basis. Denoting the element size with and the polynomial degree with

, the largest eigenvalues grow asymptotically as

for the convection operator and for the diffusion operator, see [17, Ch. 7.3]

. For explicit time integration schemes these estimates imply stability restrictions of the form

and , respectively. This corresponds to a reduction of the admissible time step by a factor of for convection and for diffusion in comparison to low-order finite-element or finite-volume methods. These stability issues lead to a preference of implicit time integration methods, especially for diffusion. With convection, however, nonlinearity may impair implicit methods and render semi-implicit or even fully explicit approaches attractive.

Time integration methods applied in computational fluid dynamics cover a wide range of approaches, including multistep, Runge-Kutta (RK), Rosenbrock, extrapolation, deferred correction and variational methods, see e.g. [46, 45, 42, 86, 28]. Early work on high-order space discretization for incompressible flows advocated semi-implicit multistep methods [73, 55] which are based on the projection method introduced by Chorin [23]. These methods were investigated and generalized in numerous follow-up studies, and gained considerable popularity, mainly because of their simplicity and low cost per time step [44, 43, 63, 37, 58, 35]. Their advantage is offset, however, by aggravating stability restrictions of the explicit part and the loss of A-stability of the implicit part for convergence orders greater than two (second Dahlquist barrier [45]). Accordingly, the vast majority of studies based on multistep methods uses order two in time, whereas the spatial order ranges from 4 to well above 10. As a notable exception Klein et al. [57] used a fully implicit SIMPLE method based on backward differentiation formulas (BDF) up to order 4.

In contrast to linear multistep methods, implicit Runge-Kutta methods can be constructed to reach high convergence orders along with excellent stability properties. For flow problems, however, convection or solution dependent viscosity yield nonlinear equations which need to be linearized and solved on every stage. In order to constrain complexity, all studies known to the author used only diagonally implicit Runge-Kutta (DIRK) methods. Uranga et al. [87] applied a third order DIRK for large-eddy simulations of transitional flow wings in conjunction with Newton’s method and preconditioned conjugate gradients for solving the nonlinear equation systems on each stage. Rosenbrock-type methods are build on the Jacobian of the right-hand side (RHS) and, thus, achieve linearization more directly. John et al. [52] compared Rosenbrock methods of order 3 for incompressible Navier-Stokes problems with fractional step methods and showed their competitiveness, especially in terms of accuracy and robustness. More recently, Bassi et al. [9] and Noventa et al. [72] applied Rosenbrock methods of orders up to 6 combined with discontinuous Galerkin methods in space to compressible and incompressible flows past airfoils and other configurations. In [72] the authors mention the possibility of order reduction with time-dependent boundary conditions, but did not include corresponding convergence studies. As an alternative to the fully implicit approach, implicit-explicit (IMEX) methods combine implicit RK for the stiff (and often linear) part of the RHS with explicit RK schemes for the nonlinear part [8, 56, 20, 13]. So far, however, IMEX RK methods appear to be rarely used in flow simulations featuring high-order spatial approximations [74, 48, 39]. Moreover, the increasing number of matching conditions [56] complicates the construction of higher order IMEX RK and methods of order higher than five are not available to the knowledge of the author. A further issue arising in the construction of high-order RK methods is the order reduction phenomenon, which can be triggered, e.g., by stiff source terms or time-dependent boundary conditions [18, 45, 56]. Although several approaches have been proposed to mitigate this issue, there exists no general solution for complex equations such as Navier-Stokes problems [5, 6, 75]. One possibility to mitigate order reduction is to use RK methods possessing a high stage order [14, 45]. Unfortunately, DIRK and Rosenbrock methods are limited to stage order two by construction.

Extrapolation and deferred correction methods employ low-order time-integration schemes within an iterative framework to achieve convergence of (in principle) arbitrary high order [29, 33, 59, 68]. In recent years, the spectral deferred correction (SDC) method gained attention [24, 25, 26, 4, 27, 79, 88, 80, 66, 19]. It was proposed by Dutt et al. [33]

for solving the Cauchy problem for ordinary differential equations and extended by

Kress and Gustafsson [59]

to initial boundary value problems. The basic idea is to convert the differential evolution problem into a Picard integral equation which is solved by a deferred correction procedure, driven by a lower order marching scheme. In this procedure, the lower order scheme sweeps repeatedly through subintervals defined by a set of collocation points, which also serve for Langrange interpolation and integration. Choosing these points from a Gauss-type quadrature yields an SDC method that converges toward the solution of the corresponding implicit Gauss collocation method

[45]. Applying a first order corrector such as implicit Euler, each sweep ideally elevates the order by one, until reaching the maximum depending on the chosen set of collocation points. Like other high-order methods SDC is also susceptible to order reduction. However, in contrast to RK methods, the phenomenon manifests rather in a slower convergence toward the collocation solution, whereas the converged solution still attains the optimal order. Although several approaches have been proposed for accelerating convergence, the general solution to this problem remains an open issue [60, 61, 24, 25, 88]. In spite of its capability to reach arbitrary high orders and straightforward extension to parallel-in-time methods boosting the efficiency on high-performance computers [70, 12], SDC has been rarely applied to fluid dynamics problems. Moreover, most studies are confined to simple configurations with constant properties and periodic boundaries [67, 69, 4, 79]. Only recently, Minion and Saye [66] proposed an SDC method for computing incompressible flows with time-dependent BC.

Variational methods resemble SDC in harnessing piecewise polynomial expansions in the time direction. Unlike the latter, they achieve discretization by application of a variational principle, predominantly the discontinuous Galerkin (DG) method. Recent applications of DG in time include incompressible flows [84, 1], elasticity [85] and hyperbolic conservation laws [38]. Although inherently implicit, the method allows to incorporate semi-implicit strategies similar to SDC as proposed e.g. in [84]. Like SDC, DG methods can be based on Lagrange polynomials constructed from Gauss points. However, assuming that points are used, DG methods generally converge with order , whereas SDC methods reach order with Gauss-Lobatto-Legendre (GLL) and with Gauss-Legendre points [19].

The goal of this study is to develop an SDC method for incompressible Navier-Stokes problems which, in combination with the DG spectral element method for spatial discretization, is capable to reach arbitrary high order in time and space. As the backbone of the new method, a semi-implicit correction scheme is devised which yields a simpler and more robust pressure handling than the SDPC method proposed by Minion and Saye [66]. Moreover, the present approach allows for variable, solution-dependent viscosity as well as time-dependent boundary conditions.

The remainder of the paper is organized as follows: Section 2 summarizes the incompressible Navier-Stokes equations with variable viscosity. Section 3 reviews the spectral deferred correction method and extends it to the flow problem. Section 4 presents the spatial discretization followed by a compilation and discussion of numerical results in section 5. Section 6 concludes the paper.

2. The Navier-Stokes equations with variable viscosity

This paper considers incompressible flows with constant density and variable viscosity in a simply connected domain . The velocity satisfies the momentum (Navier-Stokes) and continuity equations


in , where represents the pressure and


the viscous stress tensor, both divided by density;

is the kinematic viscosity and an explicitly defined forcing term. The flow problem is closed be stating initial and boundary conditions


For continuity, must be divergence free and satisfy the compatibility condition


Assuming a constant viscosity and using the identity the viscous term in the momentum equation (1) can be rewritten in the following forms

(native) (7a)
(laplacian) (7b)

For variable they can be generalized to


where corresponds to the native, the laplacian, and

the rotational form, respectively. These forms are equivalent when applied to solenoidal vector fields, but not with approximate solutions that are not divergence free. This needs to be considered in the discrete case.

3. Spectral deferred correction method

3.1. General approach

This section briefly reviews the SDC method based on the model problem


with , and initial condition .

3.1.1. Preliminaries

As a prerequisite for developing the method, the time interval is divided into subintervals such that and . Further, let denote the length of the -th subinterval, the discrete solution at time and the -th approximation to . In addition to this, the intermediate times serve as collocation points for Lagrange interpolation and as quadrature points for numerical integration. Depending on the underlying quadrature rule, one or both endpoints may be dropped, see e.g. [62, 19]. The following description is based on the Gauss-Lobatto-Legendre (GLL) rule and, hence, includes both endpoints for interpolation and integration. Accordingly, represent the GLL points scaled to , the discrete solution vector, , and the corresponding Langrange interpolant at time .

3.1.2. Predictor

The initial approximation is obtained by performing a predictor sweep of the form


where and is an approximation to . Using, for example, a combination of forward and backward Euler rules based on the decomposition yields an IMEX Euler predictor with


Alternatively, the predictor can be constructed from higher order time integration schemes such as RK or multistep methods [60]. This approach may give an advantage by providing more accurate start values, but will not be investigated in frame of the present study.

3.1.3. Corrector

The goal of the corrector is to remove the error from a given approximation . For deriving the correction equation, the error function is defined as


Further, the residual function is introduced by


Differentiating and subtracting (12) and (13) yields


This equation can be rearranged using (9) and (12) to give the error equation


which is supplemented with the initial condition .

The error equation is solved numerically by means of a time integration scheme sweeping through the subintervals. For example, application of IMEX Euler yields


for , where represents the approximation to . Substituting (13) for , eliminating by means of (12) and defining the new approximate solution by finally gives the update equation


The last term in (17) is usually approximated by replacing the integrand by its Lagrange interpolant, i.e.


with . The approximate integral can be expressed in terms of a quadrature formula,


where equals the integral of the -th interpolation polynomial over subinterval , normalized with . As a consequence, the sum recovers the weights of the underlying quadrature rule and, hence, converges to the solution of the related collocation method.

The sketched SDC method attains order for a single interval and when applied to a sequence of multiple intervals. Using a first order corrector as sketched above, every sweep, ideally, elevates the order by one, until reaching the maximum order [33]. However, stiff terms and boundary conditions may affect convergence such that more iterations are required to attain the optimal order. For details and possible remedies the reader is referred to [49, 88, 66].

3.2. Application to incompressible flow

3.2.1. Considerations

The generalization of the SDC approach to the incompressible Navier-Stokes problem (15) follows a similar approach as outlined above. Starting from an identical partitioning of a given time interval, the semi-discrete solution is denoted by for the velocity and for the pressure. Before proceeding it is important to note the following differences between the flow problem and the model problem (9): Although the momentum balance (1) resembles an evolution equation for the velocity, it involves an additional variable in terms of the pressure. Moreover, the velocity is required to satisfy the continuity equation (2), which lacks a time derivative and, hence, looks like an algebraic constraint from perspective of time integration. Finally, the flow equations are subject to boundary conditions (5) that may depend on time themselves.

3.2.2. Predictor

The complex nature of the flow problem complicates the construction of the predictor (10). Instead of defining the operator directly it is more appropriate to derive its structure from a single time step across some subinterval . In analogy to the model problem, on could use the IMEX Euler method for incompressible flow, i.e.


While this scheme looks reasonably simple and elegant, it yields a coupled system for and , which renders the solution costly, especially in view of the pertinent stability restrictions and low accuracy. Therefore, it seems attractive turning to projection schemes that decouple continuity from the momentum balance. These schemes employ some approximation of pressure in the (incomplete) momentum step and achieve continuity by performing a separate projection step [41, 43, 66]. Moreover, the latter yields a Poisson equation for correcting the pressure.

Depending on the order of the substeps two classes of projection schemes can be distinguished: pressure-correction and velocity-correction methods. Pressure-correction methods were introduced by Chorin [23]. They first solve an implicit diffusion problem for each velocity component, including approximations for convection and pressure terms, and then project the provisional velocity to a divergence-free field. In comparison to IMEX Euler, the splitting leads to an additional error caused by the violation of tangential boundary conditions in the projection step. However, several approaches exist for controlling the splitting error and retaining first order convergence [43]. Minion and Saye [66] investigated different variants of the pressure-correction scheme as a basis for their SDPC method.

Velocity-correction methods were introduced by Orszag et al. [73] and further extended, e.g. in [55, 44]. As a common feature, these methods perform the projection step before solving the diffusion problem. With constant viscosity this approach preserves continuity. However, similar to the pressure-correction method, it introduces a splitting error due to inaccurate pressure boundary conditions in the projection step. Using the rotational form of the velocity-correction method mitigates this error and recovers optimal convergence [44].

Several authors adapted the pressure-correction method to simulate flows with variable viscosity [34, 71, 31], whereas the author is not aware of corresponding extensions of the velocity-correction method. On the other hand, the pressure-correction method implies a rather complicated handling of the pressure when applied as a base method for SDC [66]. The method proposed in the following combines the advantages of both approaches. It starts with a velocity-correction step and concludes with a projection like the pressure-correction method.

For stating the base time-integration method, the momentum equation is rewritten in the form






The predictor is then defined as follows:


for . The first three substeps comprise a velocity-correction method. In particular, (30) represents an incomplete Euler step using the forward rule for convection and diffusion, backward rule for the forcing term and skipping the pressure part. It is followed by the first projection (31) and the viscous correction (32). The latter computes a new semi-implicit approximation of the diffusion term , while the explicit guess and the viscous divergence contribution introduced in (30) are dropped. For constant viscosity (3032) reproduce with the standard and with the rotational velocity-correction method as defined in [44]. In the case of variable viscosity, the diffusion step produces a divergence error of the order . This error is removed by the final projection step (33). Alternatively, it can be tolerated as a part of the overall discretization error, which will be considered as an option in the numerical experiments.

In contrast to pressure-correction, the proposed method requires no initial approximation of the pressure. The intermediate pressure and the final pressure are obtained each by solving a Poisson problem which follows from the corresponding projection step. For example, taking the divergence and, respectively, the normal projection of the first equation in (31) leads to


Similarly, (33) yields a Poisson equation for with homogeneous Neumann conditions.

3.2.3. Corrector

Apart from the additional low- and high-order contributions, the corrector resembles the predictor. The substeps for sweep read




represents the subinterval integral similar to (19) with the pressure part yet to be defined. In contrast to the predictor, the corrector exploits the previous approximation of to provide a better starting value for the implicit diffusion term in the extrapolation step (35). The pressures computed in the projection steps (36) and (38), respectively, are generally no approximations of and, hence, marked by a tilde. Adding and rearranging the equations for , , and gives


and leads to the following ansatz for the low-order contribution


Substituting the high- and low-order contributions (39, 41) in (40) and considering the converged case yields


Comparing this result to the corresponding collocation formulation implies


As a consequence it may be conjectured that the choice of and has no influence on the corrector and, hence, does not affect the convergence of the SDC method. This conjecture is confirmed by preliminary studies exploring several approaches, including the evaluation of using the recomputed pressure obtained from


Consequently, all studies in this work were performed with the simplest choice, and . Finally it is noted that the high-order contribution (39) includes the viscous divergence contribution . While this term vanishes for the exact solution, it improves stability and accuracy with non-solenoidal approximations.

4. Spatial discretization

The semi-discrete SDC formulation developed in the previous section is discretized in space using the discontinuous Galerkin spectral element method (DG-SEM) with nodal base functions [47]. The method is based largely on the approach of Fehn et al. [35], but also includes measures to attain pressure robustness proposed in [53, 77, 58, 3].

Note that the following description is constrained to cuboidal domains. This restriction serves only for convenience and can be lifted easily without affecting the proposed SDC method.

4.1. Preliminaries

First, the computational domain is decomposed into rectangular hexahedral elements to obtain the discrete domain


Let denote the set of all interior (including periodic) faces in and the set of boundary faces. The union of these sets defines the skeleton . For any interior face there exist two adjoining elements and with unit normal vectors and , respectively. The standard average and jump operators for element-wise continuous functions of any dimension are defined as


where are the traces of the function from within . A further jump operator, involving the inner product with the normal vectors, is introduced for vector or higher rank tensor functions:


These definitions are extended to boundary faces by assuming and providing exterior values on depending on boundary conditions [47]. Hereafter, the index is dropped to indicate a quantity that is defined on any face or a set of faces.

Let denote the tensor-product space of all polynomials on with degree less or equal in each direction. Glueing all element spaces together yields the global space of element-wise polynomial, discontinuous functions


This allows to define the ansatz spaces


for velocity and pressure, respectively. Note that is required for inf-sup stability, see e.g. [11, 51]. In the following the degree is set to for velocity and for pressure.

4.2. DG-SEM formulation

4.2.1. Gradient and divergence functionals

Consider the test functions and . The gradient functional of a scalar is given by




define the divergence functionals for any vector and second rank tensor , respectively. Boundary conditions are considered by providing proper exterior values, as will be detailed below. Based on these functionals the discrete gradient and divergence are introduced such that


4.2.2. Time derivative

The time derivative comprises the discrete counterparts of the convection, diffusion, pressure and forcing terms introduced in (23 – 25). For the convection term application of the local Lax-Friedrichs flux leads to


where . The diffusive and pressure terms are based on the divergence and gradient functionals, i.e.


for all . Finally, the forcing term follows from


which is equivalent to element-wise projection.

4.2.3. Laplace and viscous diffusion operators

The Laplacian occurring in the pressure equations such as (34) and the viscous diffusion operator in (32) and (37) are discretized using the symmetric interior penalty (SIP) method [7].

Application to the pressure Laplacian yields


for . The penalty parameter is defined as with


where is the mesh spacing normal to the face and a constant parameter [82]. On Neumann boundaries the face averages and jumps are given by


Note that the latter implies .

For the diffusion operator with constant the SIP method gives


where and . Dirichlet boundary conditions are weakly imposed by setting


The first of these relations is equivalent to and thus also defines the jump . Remarkably, does not couple across such that the corresponding diffusion problems can be solved component by component as long as and the RHS do not depend on the solution.

4.2.4. Divergence/mass-flux stabilization.

The discrete projection steps are augmented with the penalty functional


This functional was introduced by Joshi et al. [53] in frame of a post-processing technique for stabilizing pressure-correction methods for incompressible inviscid flow. It has no counterpart in the differential formulation, but vanishes for continuous, element-wise divergence-free vector fields . In the general case, the first part of penalizes the divergence of within elements and the second part jumps of the normal flux across faces. Akbas et al. [3] recently proved that both parts are required for pressure robustness.

The divergence penalty functional (68) has been applied with projection methods as well as coupled methods [53, 58, 3, 36]. In these studies, various expressions have been proposed for the stabilization parameters and . The present work follows [3] by setting


where is a positive constant, the mesh spacing in the normal direction and a reference value of viscosity.

4.2.5. Predictor

Starting from at time the predictor sweeps through all subintervals , performing the following steps:

  1. Extrapolation:

  2. First projection: Determine pressure by solving


    with Neumann boundary conditions


    Subsequently compute the intermediate velocity such that


    with homogeneous Neumann conditions (constant extrapolation) on .

  3. Diffusion: Find such that for all


    with Dirichlet conditions

  4. Final projection: Solve


    for and, subsequently,


    to obtain the velocity . For both problems, (76) as well as (77), homogeneous Neumann conditions are imposed.

4.2.6. Corrector

The DG-SEM formulation of the corrector resembles that of the predictor. It is summarized below, skipping specifications of function spaces and homogeneous boundary conditions for brevity.

For sweep through all subintervals and perform the following steps:

  1. [start=0]

  2. Low- and high-order contributions:

  3. Extrapolation: