1. Introduction
Highorder 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 highfidelity simulations at a scale beyond the reach of loworder methods. In simulations of processes evolving in spacetime, 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 highorder 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 highorder time integration the implications of highorder space discretization shall be briefly recapitulated. In CFD, element based Galerkin methods with piecewise polynomial expansions represent the prevalent approach apart from spectral, finitedifference 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 highorder spectral element approximations achieve far lower dispersion and dissipation errors than secondorder finiteelement or finitevolume 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 loworder finiteelement or finitevolume 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 semiimplicit or even fully explicit approaches attractive.Time integration methods applied in computational fluid dynamics cover a wide range of approaches, including multistep, RungeKutta (RK), Rosenbrock, extrapolation, deferred correction and variational methods, see e.g. [46, 45, 42, 86, 28]. Early work on highorder space discretization for incompressible flows advocated semiimplicit multistep methods [73, 55] which are based on the projection method introduced by Chorin [23]. These methods were investigated and generalized in numerous followup 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 Astability 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 RungeKutta 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 RungeKutta (DIRK) methods. Uranga et al. [87] applied a third order DIRK for largeeddy simulations of transitional flow wings in conjunction with Newton’s method and preconditioned conjugate gradients for solving the nonlinear equation systems on each stage. Rosenbrocktype methods are build on the Jacobian of the righthand side (RHS) and, thus, achieve linearization more directly. John et al. [52] compared Rosenbrock methods of order 3 for incompressible NavierStokes 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 timedependent boundary conditions, but did not include corresponding convergence studies. As an alternative to the fully implicit approach, implicitexplicit (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 highorder 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 highorder RK methods is the order reduction phenomenon, which can be triggered, e.g., by stiff source terms or timedependent 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 NavierStokes 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 loworder timeintegration 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 Gausstype 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 highorder 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 parallelintime methods boosting the efficiency on highperformance 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 timedependent 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 semiimplicit 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 GaussLobattoLegendre (GLL) and with GaussLegendre points [19].
The goal of this study is to develop an SDC method for incompressible NavierStokes 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 semiimplicit 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, solutiondependent viscosity as well as timedependent boundary conditions.
The remainder of the paper is organized as follows: Section 2 summarizes the incompressible NavierStokes 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 NavierStokes 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 (NavierStokes) and continuity equations
(1)  
(2) 
in , where represents the pressure and
(3) 
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(4)  
(5) 
For continuity, must be divergence free and satisfy the compatibility condition
(6) 
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)  
(7c) 
For variable they can be generalized to
(8) 
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
(9) 
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 GaussLobattoLegendre (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
(10) 
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
(11) 
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
(12) 
Further, the residual function is introduced by
(13) 
Differentiating and subtracting (12) and (13) yields
(14) 
This equation can be rearranged using (9) and (12) to give the error equation
(15) 
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
(16) 
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
(17) 
The last term in (17) is usually approximated by replacing the integrand by its Lagrange interpolant, i.e.
(18) 
with . The approximate integral can be expressed in terms of a quadrature formula,
(19) 
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 NavierStokes problem (1–5) follows a similar approach as outlined above. Starting from an identical partitioning of a given time interval, the semidiscrete 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.
(20)  
(21) 
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: pressurecorrection and velocitycorrection methods. Pressurecorrection 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 divergencefree 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 pressurecorrection scheme as a basis for their SDPC method.
Velocitycorrection 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 pressurecorrection method, it introduces a splitting error due to inaccurate pressure boundary conditions in the projection step. Using the rotational form of the velocitycorrection method mitigates this error and recovers optimal convergence [44].
Several authors adapted the pressurecorrection method to simulate flows with variable viscosity [34, 71, 31], whereas the author is not aware of corresponding extensions of the velocitycorrection method. On the other hand, the pressurecorrection 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 velocitycorrection step and concludes with a projection like the pressurecorrection method.
For stating the base timeintegration method, the momentum equation is rewritten in the form
(22) 
where
(23) 
with
(24)  
(25)  
(26)  
(27)  
(28)  
(29) 
The predictor is then defined as follows:
(30)  
(31)  
(32)  
(33) 
for . The first three substeps comprise a velocitycorrection 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 semiimplicit approximation of the diffusion term , while the explicit guess and the viscous divergence contribution introduced in (30) are dropped. For constant viscosity (30 – 32) reproduce with the standard and with the rotational velocitycorrection 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 pressurecorrection, 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
(34) 
Similarly, (33) yields a Poisson equation for with homogeneous Neumann conditions.
3.2.3. Corrector
Apart from the additional low and highorder contributions, the corrector resembles the predictor. The substeps for sweep read
(35)  
(36)  
(37)  
(38) 
where
(39) 
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
(40) 
and leads to the following ansatz for the loworder contribution
(41) 
Substituting the high and loworder contributions (39, 41) in (40) and considering the converged case yields
(42) 
Comparing this result to the corresponding collocation formulation implies
(43) 
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
(44)  
(45) 
Consequently, all studies in this work were performed with the simplest choice, and . Finally it is noted that the highorder contribution (39) includes the viscous divergence contribution . While this term vanishes for the exact solution, it improves stability and accuracy with nonsolenoidal approximations.
4. Spatial discretization
The semidiscrete SDC formulation developed in the previous section is discretized in space using the discontinuous Galerkin spectral element method (DGSEM) 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
(46) 
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 elementwise continuous functions of any dimension are defined as
(47)  
(48) 
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:
(49) 
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 tensorproduct space of all polynomials on with degree less or equal in each direction. Glueing all element spaces together yields the global space of elementwise polynomial, discontinuous functions
(50) 
This allows to define the ansatz spaces
(51) 
for velocity and pressure, respectively. Note that is required for infsup stability, see e.g. [11, 51]. In the following the degree is set to for velocity and for pressure.
4.2. DGSEM formulation
4.2.1. Gradient and divergence functionals
Consider the test functions and . The gradient functional of a scalar is given by
(52) 
Further,
(53)  
(54) 
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
(55)  
(56) 
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 LaxFriedrichs flux leads to
(57)  
where . The diffusive and pressure terms are based on the divergence and gradient functionals, i.e.
(58)  
(59)  
(60)  
(61) 
for all . Finally, the forcing term follows from
(62) 
which is equivalent to elementwise 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
(63)  
for . The penalty parameter is defined as with
(64) 
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
(65) 
Note that the latter implies .
For the diffusion operator with constant the SIP method gives
(66)  
where and . Dirichlet boundary conditions are weakly imposed by setting
(67) 
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/massflux stabilization.
The discrete projection steps are augmented with the penalty functional
(68) 
This functional was introduced by Joshi et al. [53] in frame of a postprocessing technique for stabilizing pressurecorrection methods for incompressible inviscid flow. It has no counterpart in the differential formulation, but vanishes for continuous, elementwise divergencefree 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
(69) 
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:

Extrapolation:
(70) 
First projection: Determine pressure by solving
(71) with Neumann boundary conditions
(72) Subsequently compute the intermediate velocity such that
(73) with homogeneous Neumann conditions (constant extrapolation) on .

Diffusion: Find such that for all
(74) with Dirichlet conditions
(75)
4.2.6. Corrector
The DGSEM 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:

[start=0]

Low and highorder contributions:
(78) (79) 
Extrapolation:
(80)
Comments
There are no comments yet.