A coupled implicit-explicit time integration method for compressible unsteady flows

TL;DR

This paper develops a stable coupled implicit-explicit time integration method for compressible unsteady flows, enabling efficient large eddy simulation and RANS computations in industrial applications.

Contribution

It introduces a novel coupling procedure for explicit and implicit schemes that maintains stability and spectral properties suitable for industrial flow simulations.

Findings

The alternative coupling procedure is stable and spectrally consistent.

Validation with standard test cases confirms effectiveness.

The method enables efficient LES/RANS coupled simulations.

Abstract

This paper addresses how two time integration schemes, the Heun's scheme for explicit time integration and the second-order Crank-Nicolson scheme for implicit time integration, can be coupled spatially. This coupling is the prerequisite to perform a coupled Large Eddy Simulation / Reynolds Averaged Navier-Stokes computation in an industrial context, using the implicit time procedure for the boundary layer (RANS) and the explicit time integration procedure in the LES region. The coupling procedure is designed in order to switch from explicit to implicit time integrations as fast as possible, while maintaining stability. After introducing the different schemes, the paper presents the initial coupling procedure adapted from a published reference and shows that it can amplify some numerical waves. An alternative procedure, studied in a coupled time/space framework, is shown to be stable and with spectral properties in agreement with the requirements of industrial applications. The coupling technique is validated with standard test cases, ranging from one-dimensional to three-dimensional flows.

Tables

Table 1. Table 1 : Flux formula depending on neighbour cells. The table is ”symmetrical”: the flux between two cells is independent of the direction of information propagation.

	Heun	Hybrid	IRK2
Heun	$F^{Heun}$	$F^{Heun}$	$\times$
Hybrid	$F^{Heun}$	$F^{Hybrid}$	$F^{IRK2}$
IRK2	$\times$	$F^{IRK2}$	$F^{IRK2}$

Table 2. Table 2 : Cell position, legend, and associated flux types. Note that left and right sides concern a 1D domain with a positive advection speed.

Legend	Flux type on left interface	Flux type on right interface	Cell position
Heun	explicit	explicit	$j=31$
Heun/Hyb	explicit	hybrid	$j=100$
Hyb/IRK2	hybrid	implicit	$j=104$
IRK2	implicit	implicit	$j=150$
IRK2/Hyb	implicit	hybrid	$j=196$
Hyb/Heun	hybrid	explicit	$j=200$

Table 3. Table 3 : Elapsed CPU time to reach 3 time periods with AION and IRK2 time integrators .

DOFs	${64}^2$	${128}^2$	${256}^2$
I/E Cells	$18\%$	$18\%$	$18\%$
AION/IRK2	$0.75$	$0.66$	$0.52$

Table 4. Table 4 : Total error slopes

Time integrator	$log\left(E_{total}\right)$ slope
Heun	1.87
AION	2.22
IRK2	2.22

Table 5. Table 5 : Total error slopes for l ​ o ​ g ​ ( h ) = [ − 2.4 , − 2.1 ] 𝑙 𝑜 𝑔 ℎ 2.4 2.1 log(h)=[-2.4,-2.1]

Time integrator	$log\left(E_{total}\right)$ slope
Heun	1.2
AION	1.84
IRK2	1.91

Table 6. Table 6 : Total error slopes for l ​ o ​ g ​ ( h ) = [ − 2.1 , − 1.8 ] 𝑙 𝑜 𝑔 ℎ 2.1 1.8 log(h)=[-2.1,-1.8]

Time integrator	$log\left(E_{total}\right)$ slope
Heun	2.54
AION	2.59
IRK2	2.54

Table 7. Table 7 : Elapsed CPU time to reach 1 time period with AION and IRK2 time integrators on a 260 2 superscript 260 2 260^{2} irregular mesh .

	CFL=0.1	CFL=0.4	CFL=0.9
Heun	1	$\times$	$\times$
AION	3.405	0.897	0.412
IRK2	3.64	0.948	0.434

Equations

$$\frac{\partial W}{\partial t}+\nabla \cdot F(W)+\nabla \cdot G(W,\nabla W)=0,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\iiint }_{{\Omega }_j}W\mathrm{d}\Omega =-{\iint }_{A_j}F(W)\cdot \vec{n}\mathrm{d}S-{\iint }_{A_j}G(W,\nabla W)\cdot \vec{n}\mathrm{d}S,$$

$${\overline{W}}_j=\frac{1}{|{\Omega }_j|}{\iiint }_{{\Omega }_j}W\mathrm{d}\Omega .$$

$$\frac{\mathrm{d}{\overline{W}}_j}{\mathrm{d}t}=R({\overline{W}}_j),$$

$$V^L=$$

$$V^R=$$

$$\{\begin{array}{rl}& \frac{\mathrm{d}{\overline{W}}_j(t)}{\mathrm{d}t}=R({\overline{W}}_j),\forall t\in R^{+},\\ & {\overline{W}}_j(0)=W_j^0,\forall j\in \{1,...,N\},\end{array}$$

$$\left\{\begin{array}[]{ll}\hbox{Predictor step:}&\displaystyle\widehat{W}=W^{n}+{\Delta t}R(W^{n}),\\ \vspace*{-3mm}&\\ \hbox{Corrector step:}&\displaystyle W^{n+1}=W^{n}+\frac{\Delta t}{2}\Big{(}R(W^{n})+R(\widehat{W})\Big{)}.\\ \end{array}\right.$$

$$\displaystyle W^{n+1}=W^{n}+\frac{\Delta t}{2}\Big{(}R(W^{n})+R(W^{n+1})\Big{)}.$$

$$f(X)=X-W^{n}-\frac{\Delta t}{2}\Big{(}R(X)+R(W^{n})\Big{)}.$$

$${f^{\prime }(X)}_{{|}_{X^p}}(X^{p+1}-X^p)=-f(X^p).$$

$$\{\begin{array}{rl}{X}^0& =W^n\\ (I-\frac{\Delta t}{2}J)\delta X^p& =-f(X^p)\\ X^{p+1}& =\delta X^p+X^p\end{array}$$

$$\left\{\begin{array}[]{ll}\hbox{Predictor step:}&\displaystyle\widehat{W}_{j}=W^{n}_{j}+\omega_{j}{\Delta t}R(W^{n}_{j})\\ \vspace*{-3mm}&\\ \hbox{Corrector step:}&\displaystyle W^{n+1}_{j}=W^{n}_{j}+\frac{\Delta t}{|\Omega_{j}|}(F_{j+\frac{1}{2}}^{Hybrid}-F_{j-\frac{1}{2}}^{Hybrid})\\ \end{array}\right.\text{for }0\leq\omega_{j}\leq 1.$$

$$\begin{array}[]{lrl}{\color[rgb]{0,0,0}{V}}^{L}=&\omega_{j}&\displaystyle\bigg{[}\frac{{\color[rgb]{0,0,0}{V}}^{n}_{j}+\widehat{{\color[rgb]{0,0,0}{V}}}_{j}}{2}+(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{j}^{n}\cdot\overrightarrow{C_{j}C_{f}}\bigg{]}+\\ &\vspace*{-3mm}\hfil\\ &(1-\omega_{j})&\displaystyle\bigg{[}{\color[rgb]{0,0,0}{V}}^{n+1}_{j}+(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{j}^{n+1}\cdot\overrightarrow{C_{j}C_{f}}-\frac{1-\omega_{j}}{2}(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})^{n+1}_{j}\bigg{]},\\ {\color[rgb]{0,0,0}{V}}^{R}=&\omega_{i}&\displaystyle\bigg{[}\frac{{\color[rgb]{0,0,0}{V}}^{n}_{i}+\widehat{{\color[rgb]{0,0,0}{V}}}_{i}}{2}+(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{i}^{n}\cdot\overrightarrow{C_{i}C_{f}}\bigg{]}+\\ &\vspace*{-3mm}\hfil\\ &(1-\omega_{i})&\displaystyle\bigg{[}{\color[rgb]{0,0,0}{V}}^{n+1}_{i}+(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{i}^{n+1}\cdot\overrightarrow{C_{i}C_{f}}-\frac{1-\omega_{i}}{2}(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})^{n+1}_{i}\bigg{]}.\end{array}$$

$$\Delta t=\underset{j\in \mathcal{D}}{\min }\hat{\nu }\frac{h_j}{\Vert \overrightarrow{v_j}\Vert +c_j},$$

$$\displaystyle\omega_{j}=\min\bigg{(}1,\frac{1}{\nu_{j}}\bigg{)},\hbox{ with }\nu_{j}=\frac{(\|\vec{v_{j}}\|+c_{j})\,\Delta t}{h_{j}}.$$

$$(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})^{n+1,s}_{j}=\left\{\begin{array}[]{ll}\displaystyle\max\left[0,\min\left(\min_{j}[\beta({\color[rgb]{0,0,0}{V}}_{j}^{n+1,s}-{\color[rgb]{0,0,0}{V}}_{j}^{n})+(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})^{n+1,s-1}_{j}],{\color[rgb]{0,0,0}{V}}_{j}^{n+1,s}-{\color[rgb]{0,0,0}{V}}_{j}^{n}\right)\right]&\\ &\hbox{ if }{\color[rgb]{0,0,0}{V}}_{j}^{n+1,s}-{\color[rgb]{0,0,0}{V}}_{j}^{n}\geq 0,\\ &\vspace*{-0mm}\\ \displaystyle\min\left[0,\max\left(\max_{j}[\beta({\color[rgb]{0,0,0}{V}}_{j}^{n+1,s}-{\color[rgb]{0,0,0}{V}}_{j}^{n})+(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})^{n+1,s-1}_{j}],{\color[rgb]{0,0,0}{V}}_{j}^{n+1,s}-{\color[rgb]{0,0,0}{V}}_{j}^{n}\right)\right]&\\ &\hbox{ if }{\color[rgb]{0,0,0}{V}}_{j}^{n+1,s}-{\color[rgb]{0,0,0}{V}}_{j}^{n}<0.\\ \end{array}\right.$$

$$\beta=\frac{2\big{(}1-\omega_{j}^{s}\nu_{j}^{s}(2-\omega_{j}^{s}\nu_{j}^{s})\big{)}}{\nu_{j}^{2}(1-\omega_{i}^{s})^{2}}$$

$$\left\{\begin{array}[]{ll}\hbox{Predictor step:}&\displaystyle\widehat{W_{j}}=W_{j}^{n}+\omega_{j}{\Delta t}R({W}_{j}^{n})\\ \vspace*{-3mm}&\\ \hbox{Corrector step:}&\displaystyle\left\{\begin{array}[]{ll}\displaystyle W^{n+1}_{j-1}=W^{n}_{j-1}+\frac{\Delta t}{2|\Omega_{j-1}|}(F_{j-\frac{1}{2}}^{n}+\widehat{F}_{j-\frac{1}{2}}-F_{j-\frac{3}{2}}^{n}-\widehat{F}_{j-\frac{3}{2}})\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j}=W^{n}_{j}+\frac{\Delta t}{|\Omega_{j}|}(F_{j+\frac{1}{2}}^{Hybrid}-\frac{F_{j-\frac{1}{2}}^{n}+\widehat{F}_{j-\frac{1}{2}}}{2})\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j+1}=W^{n}_{j+1}+\frac{\Delta t}{|\Omega_{j+1}|}(\frac{F_{j+\frac{3}{2}}^{n}+F_{j+\frac{3}{2}}^{n+1}}{2}-F_{j+\frac{1}{2}}^{Hybrid})\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j+2}=W^{n}_{j+2}+\frac{\Delta t}{2|\Omega_{j+2}|}(F_{j+\frac{5}{2}}^{n}+F_{j+\frac{5}{2}}^{n+1}-F_{j+\frac{3}{2}}^{n}-F_{j+\frac{3}{2}}^{n+1}).\end{array}\right.\end{array}\right.$$

$${(\nabla{\color[rgb]{0,0,0}{V}})}^{hyb}_{j}=\omega_{j}\bigg{(}\frac{(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})^{n}_{j}+(\widetilde{\widehat{\nabla{\color[rgb]{0,0,0}{V}}}})_{j}}{2}\bigg{)}+(1-\omega_{j})(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})^{n+1}_{j}\quad\text{for }0.5<\omega_{j}<1.$$

$$\left\{\begin{array}[]{ll}\displaystyle\frac{\partial y}{\partial t}+c\frac{\partial y}{\partial x}=0&\\ \vspace*{-3mm}\\ y(x,0)=\exp(ikx)&x\in[0,L]\\ \vspace*{-3mm}\\ y(0,t)=y(L,t)&t\in R^{+}\\ \end{array}\right.$$

$$y_j^{n+1}=G_j{y}_j^n,$$

$$\left\{\begin{array}[]{ll}\displaystyle\omega_{j}=\alpha\,\omega_{j-1}&\hbox{ for }\displaystyle\frac{N}{2}-50\leq j\leq\frac{N}{2}\\ \vspace*{-3mm}&\\ \displaystyle\omega_{j}=\frac{1}{\alpha}\,\omega_{j-1}&\hbox{ for }\displaystyle\frac{N}{2}+1\leq j\leq\frac{N}{2}+50\\ \vspace*{-3mm}&\\ \displaystyle\omega_{j}=1&\hbox{elsewhere.}\end{array}\right.$$

$$\left\{\begin{array}[]{ll}\hbox{Predictor step:}&\displaystyle\widehat{W_{j}}=W_{j}^{n}+{\Delta t}R(W_{j}^{n})\\ \vspace*{-3mm}&\\ \hbox{Corrector step:}&\displaystyle\left\{\begin{array}[]{ll}\displaystyle W^{n+1}_{j-1}=W^{n}_{j-1}+\frac{\Delta t}{2|\Omega_{j-1}|}(F_{j-\frac{1}{2}}^{n}+\widehat{F}_{j-\frac{1}{2}}-F_{j-\frac{3}{2}}^{n}-\widehat{F}_{j-\frac{3}{2}})\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j}=W^{n}_{j}+\frac{\Delta t}{|\Omega_{j}|}(F_{j+\frac{1}{2}}^{Hybrid}-\frac{F_{j-\frac{1}{2}}^{n}+\widehat{F}_{j-\frac{1}{2}}}{2})\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j+1}=W^{n}_{j+1}+\frac{\Delta t}{|\Omega_{j+1}|}(\frac{F_{j+\frac{3}{2}}^{n}+F_{j+\frac{3}{2}}^{n+1}}{2}-F_{j+\frac{1}{2}}^{Hybrid})\\ \vspace*{-3mm}&\\ \displaystyle W^{n+1}_{j+2}=W^{n}_{j+2}+\frac{\Delta t}{2|\Omega_{j+2}|}(F_{j+\frac{5}{2}}^{n}+F_{j+\frac{5}{2}}^{n+1}-F_{j+\frac{3}{2}}^{n}-F_{j+\frac{3}{2}}^{n+1}).\end{array}\right.\end{array}\right.$$

$$\begin{array}[]{ll}{\color[rgb]{0,0,0}{V}}^{L}=&\omega_{j}\displaystyle\bigg{[}\frac{{\color[rgb]{0,0,0}{V}}^{n}_{j}+\widehat{{\color[rgb]{0,0,0}{V}}}_{j}}{2}\bigg{]}+\frac{1}{2}(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{j}^{n}\cdot\overrightarrow{C_{j}C_{f}}+\bigg{(}\omega_{j}-\frac{1}{2}\bigg{)}(\widetilde{\widehat{\nabla{\color[rgb]{0,0,0}{V}}}})_{j}\cdot\overrightarrow{C_{j}C_{f}}+\\ &(1-\omega_{j})\displaystyle\bigg{[}{\color[rgb]{0,0,0}{V}}^{n+1}_{j}+(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{j}^{n+1}\cdot\overrightarrow{C_{j}C_{f}}-\frac{1-\omega_{j}}{2}(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})^{n+1}_{j}\bigg{]},\\ {\color[rgb]{0,0,0}{V}}^{R}=&\omega_{i}\displaystyle\bigg{[}\frac{{\color[rgb]{0,0,0}{V}}^{n}_{i}+\widehat{{\color[rgb]{0,0,0}{V}}}_{i}}{2}\bigg{]}+\frac{1}{2}(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{i}^{n}\cdot\overrightarrow{C_{i}C_{f}}+\bigg{(}\omega_{i}-\frac{1}{2}\bigg{)}(\widetilde{\widehat{\nabla{\color[rgb]{0,0,0}{V}}}})_{i}\cdot\overrightarrow{C_{i}C_{f}}+\\ &(1-\omega_{i})\displaystyle\bigg{[}{\color[rgb]{0,0,0}{V}}^{n+1}_{i}+(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})_{i}^{n+1}\cdot\overrightarrow{C_{i}C_{f}}-\frac{1-\omega_{i}}{2}(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})^{n+1}_{i}\bigg{]}.\\ \end{array}$$

$$\displaystyle{\color[rgb]{0,0,0}{V}}^{L}=$$

$$\displaystyle{\color[rgb]{0,0,0}{V}}^{R}=$$

$$c_j^N=\frac{arg(G_j)}{k\Delta t},$$

$$\displaystyle V^{gN}_{j}=\frac{d\big{(}arg(G_{j})\big{)}}{dk}\frac{c}{CFL\,\Delta x}.$$

$$y=\frac{1}{\sqrt{2\pi }\sigma }.\exp -\frac{(x-x_0{)}^2}{2{\sigma }^2}.$$

$$\Delta x_j=\alpha \Delta x_{j-1}$$

Full text

1 Introduction

Turbulence occurring in an unsteady compressible fluid flow is characterized by the cascade of vortices distributed over a large range of temporal and spatial scales  [1]. When the Reynolds number $Re$ is large (typically for $Re\simeq 10^{6}-10^{9}$), the resolution of the entire turbulence spectrum using Direct Numerical Simulation is still out of reach today for industrial configurations, even on the most powerful supercomputers. Currently, two major classes of approaches are mainly used to model turbulence effects in these situations.

The first technique consists in modeling all effects of the turbulence on the flow and leads to Reynolds-Averaged Navier-Stokes (RANS) equations. All turbulence effects are represented by their averaged effects on the mean flow, and turbulence fluctuations are accounted for through a turbulence model. The RANS approach is today the preferred technique for industrial applications, even if it is unable to represent unsteady effects of the turbulence. In particular, its accuracy strongly depends on the turbulence model used to represent the turbulence effects on the mean flow. The main interest of this approach remains its relatively low CPU cost for reasonable accuracy, especially for computing boundary layers. Implicit time integrators are generally considered in that case since they allow either fast convergence to steady-state solutions or large time steps for unsteady flows. Indeed, implicit schemes are generally chosen so to be unconditionally stable. These schemes generally need the linearization of the implicit residual, which makes the Jacobian of the residual appear. Finally, this leads to look for the solution increment as the solution of a linear system. Examples of implicit time integrators are the family of implicit Runge-Kutta methods (IRK) proposed by Kunstmann [2] and Butcher [3] or the Gear’s scheme [4]. Implicit methods have the drawback to require extra computational costs compared with explicit techniques, because a linear system must be solved at all time steps.

Another way to account for turbulence effects is Large Eddy Simulation (LES). Indeed, the principle of LES is to compute the largest turbulent scales and to model the smallest scales that cannot be captured by the mesh. Of course, the LES approach introduces additional terms to account for the intractable turbulence spectrum but the universal nature of the smallest scales makes the models more general than for the RANS approach. In practice, LES requires a very fine mesh to capture turbulence scales, which partly explains why the overall CPU cost is much higher than for RANS. This ratio of CPU cost is today typically about $10^{3}-10^{5}$, depending on the amount of modeled scales in LES, and therefore on the local mesh refinement. However, the intrinsically unsteady nature of LES makes it favorable for flows for which the hypothesis for averaging turbulence effects is not fully justified as for instance in wakes behind bodies at high incidence, or for simulating transition effects. Explicit time integration schemes are often used in LES to correctly capture the unsteady phenomenon. Indeed, schemes with good spectral properties are the pillars for capturing turbulence spectrum accurately and explicit schemes enable a simple and efficient control on spectral properties (dissipation and dispersion) and accuracy (order of accuracy) [5, 6]. However, explicit schemes suffer from small time steps that are limited by the Courant-Friedrichs-Lewy (CFL) number. Hence, computing boundary layers with a fine mesh at the wall with LES becomes very demanding from the CPU time point of view.

Runge-Kutta (RK) methods [7, 8] are examples of explicit time integrators. Some RK schemes were designed to respect good mathematical properties (low dissipation and low dispersion [9], Total Variation Bounded [10, 11, 12, 13], Total Variation Diminishing [14]). Standard Runge-Kutta methods are one-step and multi-stage schemes: they compute the solution at the next discrete time from the solution at the current time using additional intermediate solutions. An obvious way to extend the order of accuracy of the approach is to account for the solution history. As the number of stages increases non-linearly with the order of accuracy for standard RK methods, Williamson [15] designed low-storage RK schemes up to fourth order with only two stages in order to decrease memory cost. In the same way, linear multi-step methods, such as Adams-Bashforth methods [16], enable the computation of the new solution using several past solutions. For non-stiff problems, explicit Runge-Kutta or linear multi-step methods give accurate results. However, in case of stiffness, such standard methods can suffer from instability [17]. Therefore, for many stiff problems, the use of explicit methods requires very small time steps for the calculation to be feasible or stable.

The formula of Navier-Stokes equations for RANS or LES are very similar, and composed of a term with a time derivative, a convection term and a diffusion term. They both include turbulence effects by means of turbulent viscosity and diffusivity (defined using either an averaging or a filtering procedure). All the above comments on RANS and LES methods suggest that both methods may be coupled to take benefit from each of them in its best domain of confidence. The idea is to compute the turbulent boundary layer with the RANS approach, in a region where the LES can be less accurate because the mesh is generally not enough refined to capture turbulence effects. Far from the wall, it could be better to switch to LES in order to capture the largest turbulence vortices. Several ways to couple RANS and LES models can be found in the literature such as the Wall-Modeled LES [18] or the Detached Eddy Simulation [19, 20].

In this context, the present study concerns the design of a time integration procedure adapted to the spatial coupling of RANS and LES regions using implicit and explicit time integrators already included in the solver called FLUSEPA111FR4009261, 27/07/2013. This solver, developed by ArianeGroup, is used in the certification process of launchers during take-off and reentry, when complex physical phenomena occur. During the last years, the solver received attention and was the subject of many improvements [21, 22, 23, 24, 25, 26] .

The starting point of our study is the Heun’s [27] explicit time integration scheme available in this code. It has produced results in agreement with given LES industrial requirements and it is in addition made compatible with local mesh refinement, as in [28]. The remaining question is how to couple this second-order explicit scheme with a second-order implicit scheme, making a hybrid version of the underlying schemes that still remains second-order accurate.

Hybridization of explicit and implicit schemes is an important subject of research. In this context, IMplicit/EXplicit (IMEX) schemes consist in using a different time integration for the different terms in the original set of equations. Stiff terms are solved implicitly and regular terms explicitly. But since only one time integration term appears in the original set of equations, the IMEX approach is generally justified by the Strang’s splitting procedure [29]. The pillar of IMEX schemes lies on the separation of time integration schemes based on the different terms of the equations (see [30] and the analysis performed by Ascher et al. [31]). However, in our case, the coupling of RANS and LES leads to spatially couple explicit and implicit schemes as implicit schemes are suitable for RANS equations and explicit time integration for LES. To our knowledge, much less attention was paid on techniques to couple explicit and implicit schemes spatially. At this point, it must be highlighted that for practical applications, the monitoring of the convergence of the implicit formulation has a strong influence on the quality of the results , especially for LES. This is a reason why a full implicit formulation for both RANS and LES areas is not retained.

In 1986, Fryxell et al. [32] built an hybrid implicit-explicit method to perform one-dimensional Lagrangian hydrodynamics simulation based on the extension of a Godunov-type scheme to implicit regime. Dai and Woodward [33] proposed an iterative approach of this framework for Euler equations. Their scheme is second-order accurate in time and space but quite complicated to implement. This is especially a consequence of the requirement of solving the solution by an iterative process. In addition, the authors mention that the extension to multi-dimensional system is still an open question. Collins et al. [34] proposed the hybrid explicit-implicit scheme for Eulerian hydrodynamics, denoted CCG, that ensures the total variation diminishing (TVD) property at all CFL numbers for a linear equation. This hybrid scheme allows a second order of accuracy in space and time on pure explicit mesh cells. Unfortunately, the CCG scheme leads to incorrect numerical solutions and fails to maintain the TVD property in case of nonlinear equations. Men’shov and Nakamura [35] proposed an extension of this formulation to non-linear hyperbolic equation with Maximal Norm Diminishing (MND) property. There still remains a strong limitation: the switching procedure of CCG and Men’shov and Nakamura schemes (from explicit to implicit and vice versa) cannot maintain a second order of accuracy in space and time over the whole mesh. Timofeev and Norouzi [36] overcame this difficulty by blending second-order accurate schemes (in space and time) for explicit, implicit and hybrid cells using a smooth parameter to couple the considered time integrators. Their method is one of the most advanced ones since it maintains TVD property. In 2006, second order accurate explicit and implicit time integrators were hybridized by Tóth et al. [37] and coupled with an adaptive mesh refinement (AMR) strategy for three-dimensional applications in magnetohydrodynamics using structured grids. The procedure is quite complex to implement and needs to blend the first and second order fluxes on the interfaces.

In 2017, May and Berger [38] proposed a technique to switch between explicit and implicit approaches called flux bounding. The coupling of explicit and implicit schemes is performed by the coupling of explicit and implicit flux contributions on the interface in the context of finite volume discretization procedure. Unfortunately, it led to a result that is not second-order accurate in time and space for ”transient” cells for which residuals were computed using explicit and implicit flux contributions. Finally, it is important to notice that Persson [39] also provided a Runge-Kutta based IMEX scheme that allows to couple spatially explicit and implicit time integration schemes. However, this technique does not respect our requirement to focus on Heun’s scheme i.e. the explicit time integrator in FLUSEPA111FR4009261, 27/07/2013.

In this paper, a new way to couple spatially second-order accurate explicit and implicit time integrators for compressible unsteady flows is introduced. The goal is to propose a coupling technique that ends to a second-order accurate scheme over the whole computational domain. The remainder of this paper unfolds as follows. In Sec. 2, the standard form of the Navier-Stokes is recalled and the explicit and implicit schemes chosen for the present study are introduced. Our explicit and implicit basic schemes, Heun and implicit Runge Kutta schemes, are defined in Sec. 3. In Sec. 4, the method presented in [36] is recalled as it will serve as the basis of the new scheme. In Sec. 5, a first version of coupling Heun and implicit Runge Kutta (IRK2) time integrators is introduced with its dedicated space-time stability analysis. The first hybridized scheme is not acceptable and a second version of the coupling scheme is presented in Sec. 6. The space-time stability analysis is performed again. Before concluding, Sec. 7 is dedicated to the validation of the new coupling scheme for 1D, 2D and 3D test cases.

2 Discretization of the Navier-Stokes Equations

The Navier-Stokes conservation laws can be written in the following compact form,

[TABLE]

where $W$ is the vector of conservative variables and $\nabla W$ its gradient, $F$ and $G$ are respectively the flux density for convection and viscous effects. The computational domain $\Omega$ is split into $N$ non-overlapping rigid immobile cells $\Omega_{j}$ and Eq. (1) are integrated over every mesh cell. The Gauss relation ties the volume integrals of the divergence terms to the interface fluxes, leading to

[TABLE]

where $\Omega_{j}$ is the $j$-th control volume with its border $A_{j}$ and $\vec{n}$ is the outgoing unit normal vector. From now on, it is possible to define the averaged quantity of the conservative variables,

[TABLE]

This formula allows to rewrite the conservation laws (2) discretized with a finite-volume formulation in the following compact form,

[TABLE]

with $R(\overline{W}_{j})$ is the residual computed using the averaged quantities in cell $j$. Actually, the residual is defined as the sum of the flux over the whole boundary of a cell. Here, convection and diffusion fluxes are discretized by means of a $k$-exact formulation coupled with successive corrections [22, 40]. Indeed, any high-order scheme needs a polynomial representation of the unknowns locally around a cell of interest. Theoretically, the polynomial approximation can also be defined as the Taylor expansion of the local solution. Standard $k$-exact schemes need the definition of a large stencil in order to link Taylor expansion coefficients and averaged quantities (see [41] as an example of standard $k$-exact scheme). This large stencil is a drawback to cope with parallel computations since the number of data to exchange between domains depends on the partitioned (initial) mesh. The successive corrections were designed in order to avoid geometrical reconstructions for parallel computations. The number of ghost cells are limited to the minimum of one ghost cell per volume owing an interface on the boundary. Starting from the solution, the polynomial reconstruction is built using the truncation error of the Taylor expansion by increasing the order of the approximation. The starting point is the availability of the averaged quantities in the mesh cells and the definition of a discrete gradient scheme. The discrete gradient scheme applied to the averaged quantities leads to a given order of accuracy. A local analysis enables to estimate the truncation error performed on the gradient and by removing the last non-zero term in the truncation error, the gradient scheme accuracy is made more accurate by one order. This gradient is then usable to define a more accurate representation of the unknown through a local polynomial reconstruction. The process can then be repeated for the Hessian in order to increase accuracy on the computation of the successive derivatives [40, 22]. Of course, since the truncation error must be computed, it is necessary to locally exchange the set of local derivatives (first, the gradients, then the second-order derivatives, then the third-order derivatives and so on). The TVD property is ensured with a MUSCL reconstruction. For a second-order scheme, the flux on an interface is estimated using an affine extrapolation of the cell-centered unknowns onto the interface. If $f$ denotes a mesh interface with a unit normal vector $\vec{n}$ directed from the left cell $j$ to the right cell $i$, and if $C_{i}$, $C_{j}$ and $C_{f}$ respectively denote the cell centers and the interface center, the standard reconstruction on the vector of primitive variables $V$, is simply,

[TABLE]

where the gradient $\nabla V$ needs to be computed with the $k$-exact reconstruction. These gradients are also used to compute the viscous flux.

Finally, the finite-volume formulation of Eq. (4) can be expressed as the following Cauchy problem,

[TABLE]

where ${W_{j}^{0}}_{1\leq j\leq N}$ denotes the initial solution. For the sake of clarity, the averaging symbol will be dropped, and ${W}$ will represent the averaged quantities over the control volumes.

3 Heun and Crank-Nicolson/IRK2 Time Integration Methods

Heun and Crank-Nicolson/IRK2 are considered as the two basic time integrators of the present study from which a new method will be developed.

The Heun’s explicit scheme [27] is a second-order accurate predictor-corrector method. The state $W^{n+1}$ is first predicted with a forward Euler scheme, and then corrected by a standard trapezoidal rule,

[TABLE]

The Crank-Nicolson/IRK2 implicit scheme [42] is second-order accurate, and it can be expressed as,

[TABLE]

The residual being non linear in $W$, Eq. (8) is solved with an iterative Newton method. Starting from the iteration $n$ in Eq. (8), the next solution is defined as the root $X$ of $f(X)=0$ with

[TABLE]

At the $p$-th step of the Newton algorithm, the next iterative $X^{p+1}$ is such that

[TABLE]

The final algorithm is simply

[TABLE]

with $J=\frac{\partial R}{\partial W}$ the Jacobian matrix and $I$ the identity matrix. The linear system in the iterative loop is solved with a Krylov iterative method. The proposed Newton procedure is motivated by the strong non-linearities encountered during reentry and take-off. For small non-linear effects, the residual $R(W^{n+1})$ can be linearized around $W^{n}$.

The open question concerns how to spatially couple both second-order time integrator schemes, while maintaining the second order of accuracy and ensuring the overall stability, and taking into account the fact that both time integrators need two residuals computed at two different times. To this aim, the hybrid time integrator proposed by Timofeev and Norouzi [36] is first considered. In the following, all numerical schemes will be expressed in a one-dimensional domain without loss of generality, and for convective terms. Diffusive terms will be discussed if required.

4 Principle of the Hybrid Time Integration

The spatially-coupled explicit and implicit scheme proposed by Men’shov and Nakamura [35] and extended by Timofeev and Norouzi [36] is recalled hereafter to make an introduction to the hybrid time integration. The hybrid scheme [36] applied to the Cauchy problem of Eq. (4) has a predictor-corrector formulation:

[TABLE]

The main components of the scheme are presented in the following using a face $f$, with the face center $C_{f}$, separating two cells, $i$ and $j$, with respectively cell centers $C_{i}$ and $C_{j}$.

During the predictor step, the residual computation does not use any Riemann solver, and the predicted state $\widehat{W}_{j}$ is computed using a fully upwind scheme. So, the flux F({\color[rgb]{0,0,0}{V}}^{f}_{j}), computed from the vector of the primitive variables ${V}$ and involved in the computation of the residual $R({W^{n}_{j}})$ of the cell $j$ uses the reconstruction on the primitive state {\color[rgb]{0,0,0}{V}}^{f}_{j}={\color[rgb]{0,0,0}{V}}^{n}_{j}+(\widetilde{\nabla{\color[rgb]{0,0,0}{V}}})^{n}_{j}\cdot\overrightarrow{C_{j}C_{f}} [43, 44, 45]. This prediction step is not conservative since the flux F({\color[rgb]{0,0,0}{V}}^{f}_{j})\neq F({\color[rgb]{0,0,0}{V}}^{f}_{i}). Moreover, the gradient \widetilde{\nabla{\color[rgb]{0,0,0}{V}}} is computed from the standard gradient {\nabla{\color[rgb]{0,0,0}{V}}} and a slope limiter.

For any cell $j$, the hybrid flux $F^{Hybrid}$ is computed using a Riemann solver with the following reconstructed left and right primitive states on a face:

[TABLE]

The parameter $\omega_{j}$, subsequently named cell status, is defined for each cell $j$ such that $0\leq\omega_{j}\leq 1$. It is used to switch smoothly from a time-explicit ($\omega_{j}=1$) to a time-implicit ($\omega_{j}=0$) scheme. The cell status $\omega_{j}$ is locally adapted to the flow and the stability constraint. First, the global time step is chosen to be the maximum allowable time step on a part $\cal{D}$ of the whole computational domain without source of stiff terms:

[TABLE]

where $\hat{\nu}<1$, $h_{j}$ is a reference length scale, $\vec{v_{j}}$ the velocity vector, and $c_{j}$ the speed of sound in cell $j$. Then, the parameter $\omega_{j}$ is defined as:

[TABLE]

In Eq. (13), the term {(\widetilde{\Delta_{t}{\color[rgb]{0,0,0}{V}}})}^{n+1}_{j} comes from a generalization of the TVD property to the coupled space/time reconstruction. It can be seen as a time slope limiter. This term is given by a Newton’s algorithm which at any step $s$ reads:

[TABLE]

The parameter $\beta$ is computed from the local CFL number $\nu_{j}$ and $\omega_{j}$ at each step $s$:

[TABLE]

This method is second-order accurate in both space and time in fully explicit, fully implicit and explicit-implicit domains. Moreover, the scheme was proved to be TVD for 1D linear/non linear scalar conservation laws. It was applied successfully to problems including shocks [36]. This scheme will be referred as TN in the following.

It is clear that the underlying explicit (where $\omega_{j}=1$) and implicit (where $\omega_{j}=0$) schemes are not Heun’s and Crank-Nicolson/IRK2 time integrators. In section 5, a first straightforward attempt to couple these latter schemes is explained.

5 First Hybrid Coupled Scheme

This section presents the coupling of three schemes: Heun for the explicit part, IRK2 for the implicit part, and the hybrid formulation TN [36] for a smooth transition between the latter two schemes. The resulting numerical time integration scheme is defined to ensure a unique flux on each interface, leading to a conservative formulation.

An interface separates two cells that may have a different status given by the parameter $\omega_{j}$. The cell status is named Heun for $\omega_{j}=1$, IRK2 for $\omega_{j}=0$, and Hybrid otherwise, as shown by the 1D example on Fig. 1.

The flux formula used on an interface separating two cells depends on the status of these cells as shown in Tab. 1.

According to Tab. 1, for the 1D example on Fig. 1, the first version of the coupling scheme, referred as HCS1, is designed as,

[TABLE]

The hybrid flux $F^{Hybrid}$ is computed with the state reconstruction method introduced in Eq. (13). The scheme Eq. (18) ensures the unique definition of the interface flux, and then flux conservation.

The differences with the TN scheme are:

• •

In the predictor step, the residual is computed from the flux balance over all cell faces. The flux is computed using a linear extrapolation of the unknowns between cell centers and cell interfaces and accounting for a slope limiter. Since two different extrapolations are obtained on a face shared by two cells, the flux is computed by means of a (approximated) Riemann solver. As a consequence, $\widehat{W}_{j}$ has a physical meaning, something which was lost with TN scheme.

• •

The scheme recovers the Heun’s scheme in the explicit regime ($\omega_{j}=1$).

• •

The scheme recovers the IRK2 scheme in the implicit regime ($\omega_{j}=0$).

Remark 1: It is very important to make the transition between explicit and implicit time integration as fast as possible. This is the reason why in practice, the IRK2 formulation is applied when $\omega_{j}\leq 0.5$.

Remark 2: For the Navier-Stokes equations, the standard diffusion scheme available in the code was used for the implicit and explicit regions. The interface gradient is defined from the limited cell-centered gradients by means of an average taking into account local metrics. In the hybrid regions, the definition of the cell-centered gradient is modified as

[TABLE]

Every gradient in Eq. (19) is computed from the standard gradient and corrected using a slope limiter.

5.1 Space-Time von Neumann Analysis

The von Neumann analysis is a powerful tool to analyse the spectral properties of any scheme. This analysis is generally performed separately either for spatial or temporal schemes but it was shown in many recent papers (see [46, 6, 5, 47] for instance) that the true scheme behaviour is only obtained by a coupled analysis involving both space and time discretizations.

Starting from the one-dimensional linear advection equation coupled with a harmonic initial condition and with periodic boundary conditions,

[TABLE]

with $i^{2}=-1$ and $c>0$, the comparison between the exact theoretical solution and the numerical approximation gives information on both dissipation and dispersion of the fully discrete scheme. The fully discrete relation obtained from Eq. (20) using a second-order finite volume scheme reads:

[TABLE]

where $y_{j}^{n}$ represents the harmonic solution at discrete time $n$ and in the center of cell $j$. The complex coefficient $G_{j}$ represents the transfer function between two consecutive time solutions. The coefficient $G_{j}$ can be expressed as $G_{j}=|G_{j}|\,\exp(i\arg(G_{j}))$, which highlights the dissipation of the scheme related to $|G_{j}|$, and its dispersion related to $\exp(i\arg(G_{j}))$. In the following, let $\mu_{j}=|G_{j}|$ be the dissipation coefficient, and $\phi_{j}=arg(G_{j})/\hbox{CFL}$ the dispersion coefficient.

For standard finite element, difference or volume formulations, and focusing on a single scheme for space and time, the transfer coefficient does not depend on the space position. However, in our case, several time discretization schemes are blended, and the spectral properties of the global scheme will depend on the cell status $\omega_{j}$. Accounting for the cell status variations can be done through the change in the mesh spacing. But taking into account for a local change in space discretization in the spectral analysis is cumbersome (see for instance [48]). Once the mesh size $\Delta x$ is fixed, a manufactured choice of $\omega_{j}$ keeps the analysis quite simple to perform. The number of 1D segments is fixed to $N=300$ and the function $\omega_{j}$ is defined as

[TABLE]

with $\alpha=0.90$. The maximum CFL value for the Heun explicit scheme without amplification at any wavenumber is 0.6. Then, the spectral analysis will be performed for CFL=0.1 and CFL=0.6.

Figures 2 and 3 give an overview of the dissipation and the dispersion of the coupled time/space HCS1 scheme depending on the cell index $1\leq j\leq N$ and the non dimensionalized wave number $0\leq k\Delta x\leq\pi$ for CFL=0.1. They clearly show that some phenomena occur at the transitions between the different temporal schemes. As a consequence, the standard curves ($\mu_{j}$ and $\phi_{j}$ as a function of $k\Delta x$) are introduced in Figs. 4-6 for the specific cells at which there is a time integration transition. As indicated in Tab. 2, six specific discretization cells are analysed. In addition, Tab. 2 also gives the meaning of the legends for Figs. 4-6.

Figure 4 shows the dissipation coefficient. The focus for $k\Delta x\in[0,0.6]$ in Fig. 5 shows that the scheme amplifies some wavenumbers for Heun/Hybrid cells. The level of amplification is very low and one could assume that since the treatment is local (on one cell in 1D), the amplification could not lead to global instability and failure of the computation. Fig. 6 shows that cells with Hybrid/Heun and Hybrid/IRK2 fluxes produce the highest levels of dispersion, while dispersion is close to the reference Heun’s scheme for other choices.

The same analysis is performed at CFL=0.6, i.e. at the limit without amplification at any wavenumber for the Heun’s scheme. The dissipation and the dispersion curves are shown in Figs. 7-8. The dispersion curve for the Heun/Hybrid cell looks like the one for Heun’s scheme and presents a change of phase sign for $k\Delta x\simeq 1.8$. The dissipation curve presents again amplification but the level of amplification is much higher. Since the value of $\omega_{j}$ may be defined to follow some physical property, numerical amplification of a physical phenomena could follow the phenomena itself, leading to an unacceptable flow or simulation break. Notice that the discrete jumps on dispersion curves in Fig. 8 at $k\Delta x\simeq 1.2,1.6,1.8,2.2,2.4$, and $2.6$. These values of $k\Delta x$ are maxima of dispersion error, and correspond to values of wave number $k\Delta x$ when $\mathcal{R}e(G_{j})=0$.

The first coupled scheme HCS1 based on Heun and IRK2 schemes presents amplification. In order to remove it, a new way to couple the desired explicit and implicit schemes is proposed in Sec. 6.

6 Another tIme integratiON (AION) scheme

In order to get better spectral characteristics than the HCS1 scheme, a new hybrid scheme still based on both Heun and IRK2 schemes is designed. This new hybrid time integrator is called AION (Another tIme integratiON). It is applied to the Cauchy problem Eq. (4).

This hybrid scheme has a predictor-corrector formulation as TN’s and Heun’s schemes. According to Tab. 1, for the 1D example on Fig. 1, it is designed as:

[TABLE]

Remark: A situation not adressed in Tab. 1 nor in Fig. 1 must be described. For any hybrid cell (according to the value $\omega_{j}$) with hybrid flux for all the faces, the predictor state $\widehat{W}_{j}$ is time-integrated by accounting for $\omega_{j}$: $\widehat{W}_{j}=W_{j}^{n}+\omega_{j}{\Delta t}R(W_{j}^{n})$. For the mixed explicit / implicit regime, the hybrid flux differs due to the different reconstruction methods which is the key point to maintain stability without amplification. The demonstration of the stability without amplification is the purpose of the next section.

The reconstructed states are defined as:

[TABLE]

For $\omega_{j}=1$ and $\omega_{i}=1$, the following reconstruction is obtained:

[TABLE]

and the hybrid part leads to the Heun’s scheme for the one-dimension linear advection equation Eq. (20).

6.1 Space-Time Stability Analysis

The analysis introduced in Sec. 5.1 is now applied to the AION scheme. The spectral behaviour is shown in Figs. 9 and 10 for each cell index $1\leq j\leq N$ and wave number $k\Delta x\in[0,\pi]$. Using the isocontours of iso-$\mu_{j}$ and iso-$\phi_{j}$, it is clear that the spectral behaviour of the scheme depends on the cell index $j$ via the cell status $\omega_{j}$.

As before, the standard curves are now presented for the 6 values of $j$ introduced in Tab. 2 at CFL=0.1 and CFL=0.6. At CFL=0.1, the overall view of the dissipation coefficient in Fig. 11 shows that the AION scheme seems to have the same behaviour for any cell index. In addition, a zoom for $k\Delta x\in[0,0.6]$ does not show any amplification of waves and stability is kept for any kind of cell (Fig. 12). Finally, all schemes present the same kind of dispersion in Fig. 13. This shows that the limitation of the initial HCS1 scheme is now removed by our new coupling technique at CFL=0.1. In addition, as shown in Figs. 14-16, the same conclusion is obtained at CFL=0.6, which is the limit without amplification at any wavenumber for the pure explicit scheme.

For this example of linear advection, the proposed scheme is shown to be stable. The value of $\omega_{j}$ is set to vary between 0 and 1. For $\omega_{j}=1$, it was shown that the reconstructed states at the interface enable to recover the standard explicit Heun’s scheme for the linear equation. So, the idea is to keep a smooth transition from explicit to hybrid regions. However, it would be convenient to switch to implicit regions as fast as possible. In the next section, attention is paid on the transition between implicit and hybrid cells.

6.2 Analysis of the Hyb/IRK2 and IRK2/Hyb Transitions

This analysis is motivated by the desire to switch as fast as possible to the implicit formulation, while maintaining stability. The computational domain with 300 cells is split into two parts. Half of the domain is dedicated to the hybrid part of the AION scheme while the rest of the domain is dedicated to the IRK2 scheme. The spectral analysis is performed at the transition between the chosen schemes. Two configurations are analysed, one for the wave going from the hybrid domain to the implicit domain, and one for the wave going from the implicit domain to the hybrid domain.

Figs. 17 and 19 show the dissipation coefficient for CFL=0.5 and for CFL=0.6. For the transition between implicit and hybrid schemes. Amplification (positive dissipation) is represented by a grey area. Amplification occurs for the transition Hyb/IRK2. In the reverse mode, amplification does not occur, as highlighted by Figs. 18 and 20. Moreover, a stable formulation is obtained for $\omega_{j}\leq 0.6$, and this motivates our final choice to switch at $\omega_{j}=0.6$ between the hybrid flux computation and the implicit Runge-Kutta scheme. As a consequence and for the following applications, $\omega_{j}=1$ will lead to the Heun’s scheme, $\omega_{j}<0.6$ will lead to the IRK2 scheme. The new hybrid reconstruction method based on the reconstructed states Eq. (24) is applied for $0.6\leq\omega_{j}\leq 1$.

6.3 Analysis of $q-$waves

Focusing on the one-dimensional linear advection equation with a constant velocity $c$, several authors mentioned and explained the existence of numerical waves travelling in the direction opposite to $c$, introducing the notion of $q-$waves. Indeed, a $p$-wave is a wave travelling in the same direction as $c$, eventually not at the same speed. $q-$waves are travelling in the opposite direction. Existence of $q$-waves was shown to be independent of the stability analysis of numerical scheme and was highlighted for instance in the seminal work of Vichnevetsky, summarised in [49], by Poinsot and Veynante for combustion [50] or by Trefethen [51]. In this context, Sengupta et al. [46] linked the numerical group velocity with the existence of $q$-waves. Following the definition proposed in [46], new quantities related to the one-dimensional linear advection equation Eq. (20) are introduced. The numerical phase speed is defined as

[TABLE]

from which is deduced the numerical group velocity,

[TABLE]

Indeed, for a positive advection velocity $c>0$, $q$-waves are propagating upstream and their group velocity $V^{gN}$ is negative [46]. It is important to notice that the negative group velocity is only a necessary condition for apparition of $q$-waves. Indeed, observation of $q$-waves also depends upon the real and imaginary part of the amplification factor of the space-time discretization and with an excessive dissipation, filtering or damping, the $q$-waves can be removed very quickly. Even if it was mentioned in [46] that initial condition, grid resolution and multi-dimensional case have effects on $q$-waves, the analysis presented in the following is restricted to the one-dimensional case.

A mesh composed of 300 cells is considered with the manufactured choice of $\omega_{j}$ introduced in Sec. 5.1 and used again in Sec. 6.1. The numerical phase speed $c_{j}^{N}$ and the numerical group velocity $V^{gN}_{j}$ are computed from $G_{j}$ (Eq. 21) for any cell $j$ (and therefore for a given $\omega_{j}$). The analysis of $q-$waves is performed for the AION scheme and for CFL $\in[0,1]$. Attention is paid on the computation of the group velocity $V^{gN}_{j}$ for the different intersection of time integrators in AION scheme. In Figs. 21-24, dashed curves correspond to $|V^{gN}|/|c|>>1$, which represents discontinuity of dispersion $\phi_{j}$, and grey zones correspond to negative group velocity $V^{gN}_{j}<0$.

Remark: The choice of the CFL domain is dictated by the desire to compare several schemes. Of course, amplification always occurs for the Heun’s scheme at CFL$>0.6$ but the analysis can still be performed for CFL$>0.6$.

Concerning the areas of negative-$V^{gN}_{j}$ group velocity, it can be noticed that the zone is slightly more expanded for Heun/Hyb cells than for the pure Heun scheme, according to Figs. 21 and 22. Indeed, negative group velocities can be found up to CFL=0.57 at Heun/Hyb cells, but the full Heun’s scheme negative-group-velocity zone is present up to CFL=0.47. In addition, this negative group-velocity zone is larger for the other transitions according to Figs. 23 and 24. So, the presence of $q$-waves is possible in a larger domain of CFL for AION scheme than for the full Heun’s scheme. The last question concerns the capability of the AION scheme to damp $q-$waves. To answer it, it is mandatory to couple the observations on negative group velocity with the dissipation property of the AION scheme, as shown in Figs. 26-30.

According to Fig. 26, the CFL limit of the negative-$V_{gN}$ zone (CFL=0.47) due to Heun time integration corresponds to $\mu_{j}<0.2$. In the case of the time integrator transitions of AION scheme, the CFL limit of their negative-$V_{gN}$ zone corresponds to $\mu_{j}<0.3$ for cell Heun/Hyb (Fig. 27), $\mu_{j}<0.5$ for hyb/IRK2 and IRK2/Hyb cells (Fig. 28 29), and finally $\mu_{j}<0.2$ for cell Hyb/Heun cells (Fig. 30). Fortunately, it seems that the effect of extension of the $q$-waves zone due to our coupling scheme can be attenuated by dissipation. All this information allows us to conclude that the AION scheme seems to be as stable as the full Heun scheme on uniform mesh. Moreover, the possible additional $q$-waves that could propagate would be quite dissipated. The good numerical behaviour obtained through the proposed theoretical analysis must now be confirmed by numerical simulations, which is the goal of the next sections.

7 Validation

This section is devoted to the validation of the AION scheme using several test cases of increasing complexity, starting from 1D advection solution and Euler solutions to 3D Navier-Stokes computations.

7.1 Advection of a Gaussian hump

The first test case is the linear advection of a Gaussian hump (Eq. 1 with $F(W)=W$ and $G(W,\nabla W)=0$ in a periodic domain) and serves to assess both stability and accuracy of the AION, Heun and IRK2 schemes. Introducing the periodic domain $[0,1]$ and the parameters of the initial hump $x_{0}=0.5m$, $\sigma=0.1$, the analytic theoretical initial solution is simply:

[TABLE]

The initial solution must be recovered at any time $t=k\,s$ with $k\in\mathbb{N}$ and here, the analysis is performed at $t=10s$. The mesh is composed of $N=333$ cells and is irregular but manufactured using the following rule:

[TABLE]

with $\alpha=0.9$. Then the parameter $\omega_{j}$ is controlled as :

[TABLE]

Such a variation in mesh size enables a strong control of the maximum CFL value authorized for a computation using an explicit time integration. In addition, the coupling scheme is fully explicit in the regular part of the grid (ie. when $\Delta x_{j}=\Delta x_{max}$). The spatial scheme used is a second order spatial scheme for regular and irregular grids. The time step is fixed in order to obtain a CFL constant in regular grid ($=\hbox{CFL}_{min}$) and leads to higher values in the smaller cells of the irregular grid.

The first computations are performed using IRK2, Heun and AION scheme at $\hbox{CFL}_{min}=0.005$ in the explicit part. The largest CFL value in the implicit part is 1. Fig. 31 illustrates the fact that the AION scheme conserves the time accuracy of Heun’s scheme but the IRK2 scheme involves more dissipation. At this level, it must be mentionned that the Newton algorithm for the convergence of IRK2 and AION scheme is guaranteed: the final error (at any time step) is very low. This is to avoid any confidence in the results due to a bad convergence of the Newton algorithm.

Now, the same test case is performed with a higher $\hbox{CFL}_{min}=0.1$. In that case, the maximum CFL number in the implicit part reaches $20$, out of the stability region of Heun’s scheme. In this configuration, both AION and IRK2 schemes are stable. Here again, the IRK2 scheme is more dissipative than the AION scheme. Moreover, for the present case, the convergence of the implicit schemes (AION and IRK2) is performed until a low error on the solution. The number of Newton steps is much larger for the IRK2 scheme than for the AION scheme.

The last question of importance concerns the interest of using the AION scheme at moderate CFL numbers versus the IRK2 scheme at larger CFL value. The computation is now performed with the IRK2 scheme at an associated $\hbox{CFL}_{min}=1$ in the explicit part. At this CFL number, the AION computation cannot be performed since Heun’s scheme is now unstable. As expected (Fig. 33), the IRK2 scheme is usable in this case but the number of iterations to converge Newton’s algorithm has a strong influence on the solution accuracy. In our example, the reference solution needs 90 steps for Newton’s algorithm and the solution is still dissipated and dispersed. In order to try to represent an industrial computation, the same simulation is now performed but only 60 steps of Newton’s algorithm are allowed. The solution is more dispersive. For the sake of consiceness, the AION solution at $\hbox{CFL}_{min}=0.1$ is recalled. For the implicit part of the scheme, the convergence of Newton’s algorithm does not require more than 10 steps.

Indeed several articles mentionned that the time integration of LES and DNS simulation by implicit time integrators involve a damping of the high-frequency solution modes for large time steps [52, 53]. Furthermore the number of Newton subiteration is not easy to control in order to obtain convergence of the error. It remains more interesting to perform implicit time integration in stiff and/or stationnary zone of the numerical domain and perform explicit time integration elsewhere.

7.2 Sod’s tube

The Sod’s tube is an unsteady inviscid one-dimensional problem used to characterize the properties of both time integration and convection schemes. The computational domain of length $L_{x}=1~{}m$ is split into two parts separated by a membrane located initially at $x=0.5L_{x}$. The initial flow is defined by

[TABLE]

where $L$ refers to the left side and $R$ to the right side of the membrane. At $t=0$, the membrane blows up and waves are travelling inside the computational domain. At the final time $t=0.2s$, the solution is composed of a rarefaction wave, a contact discontinuity and a shock.

The Euler equations are solved using a 1-exact (second order) upwind scheme using the Roe approximate Riemann solver. In order to avoid spurious oscillation and to compute TVD solution, our limiting strategy is provided with the minmod slope limiter [54]. The computational domain is composed of $N=300$ cells located in regular parts with a uniform mesh size and in irregular parts with non-uniform mesh size such that

[TABLE]

with $\alpha=0.98$. Then the parameter $\omega_{j}$ is controlled as:

[TABLE]

This definition of $\omega_{j}$ leads to a fully explicit time integration on the regular part of the grid, where $\Delta x_{j}=\Delta x_{max}$. In the following, $\nu_{j}$ is defined by $\nu_{j}=\frac{\Delta t}{\Delta x_{j}}$. For $\nu_{min}=0.1$ in the explicit part, the maximum value of $\nu$ in the implicit part reaches $0.3$. Figures 34 and 35 show the density and velocity profiles, with zooms in the region of the rarefaction wave and near the contact discontinuity. The AION scheme appears to be as accurate as the second order time accurate standard schemes and in agreement with the solution obtained using the Euler explicit scheme.

A second set of computations is performed for $\nu_{min}=0.4$ in the explicit part, which leads to a maximum value of $\nu$ in the implicit part of $1.01$. The Heun’ scheme is unstable, but Figs 36 and 37 show that the AION and IRK2 schemes are still stable. Both solutions are very close. Paying attention to the velocity and density near the shock, it seems that the AION scheme slightly dissipates the overshoots given by the IRK2 implicit scheme.

From this case, it can be concluded that the reconstruction procedure to define the AION scheme is in agreement with the requirements to handle shock, contact discontinuity and rarefaction wave. Moreover, the AION scheme improves the explicit scheme with an enhanced stability and is more computationally efficient than the fully-implicit scheme. The next test case is dedicated to the analysis of the scheme accuracy.

7.3 Two-dimensional linear advection of an isentropic vortex

The advection of an isentropic vortex is a famous test case of the High Order Workshop. This problem is designed in order to test the solver capability to preserve vorticity in an unsteady simulation and this verification case is suitable for verifying total order of accuracy. The total error $E_{total}$ may be expressed as:

[TABLE]

with $(A,B)\in\mathbb{R}^{2}$. The total order of accuracy is estimated as $min(p,q)$. For self-consistence, the test case definition is summarized below. The computational domain is the square domain $[-\frac{L}{2},\frac{L}{2}]^{2}$ ($L=0.1$) with periodic boundary conditions. On a uniform flow of pressure $P_{\infty}$, temperature $T_{\infty}$ and Mach number $M_{\infty}$ is superimposed an isentropic vortex defined by its characteristic radius $R$ and strength $\beta$. The vortex is initialized in the center of the computational domain $(x_{c},y_{c})=(0,0)$. The initial state is defined by:

[TABLE]

with:

[TABLE]

with $R_{gas}=287.15\text{ }J/kg/K$ the gas constant and a constant ratio of specific heats $\gamma$ equal to $1.4$. The isentropic relation leads to the complete set of unknowns

[TABLE]

Here, the ”fast vortex” test case is considered and defined by

[TABLE]

The solution is time-marched during three rotations inside the periodic box. The computation is performed with a Successive-Correction 2-exact formulation for the spatial scheme (order three) and time integrated by Heun, IRK2 and AION schemes on the baseline Cartesian grids of $64^{2}$, $128^{2}$ and $256^{2}$ degrees of freedom (DOF). The time integration of theses computations was performed at CFL=0.1. For the AION scheme, it is mandatory to define the hybrid parameter $\omega_{j}$ ; it is configured according to following equation:

[TABLE]

Fig. 38 represents the distribution of $\omega_{j}$ on the $256^{2}$ mesh.

According to the pressure and velocity field obtained with Heun, AION and IRK2 time integrator on grid $256^{2}$, it appears that time integrators have slightly same properties of dissipation and dispersion. The AION scheme tends to less dissipate than other second-order time integrator (according to Figs. 40 and  42).

This feature is quite different regardless the ones found on one-dimensional test case with second-order space reconstruction, where the AION was more dissipative than other time integrators. Here the third-order accurate space reconstruction tends to improve the space-time spectral characteristics of the computation, time integrated with AION scheme. According to the computation of $L^{2}$-norm of the error using the two velocity-vector components (see Fig. 40), it appears, as expected, that the AION time integrator conserves second-order space-time accuracy (see Tab. 4). Nevertheless, according to the slope of $log(E_{total})$, it seems that the time integrator has a significant impact on the error computation (see Tabs. 5-6). Indeed the value of the slope is quite different between $log(h)=[-2.4,-2.1]$ and $log(h)=[-2.1,-1.8]$, particularly in case of the Heun scheme. The computational cost of the AION and the IRK2 time integrator at same CFL is also compared in Tab. 3. The ratio between cells with $\omega_{j}<1$ and cells with $\omega_{j}=1$ (I/E cells) is equal to $18\%$ for all grids. It appears that the AION scheme reaches an overall cost that is about nearly $50\%$ lower than the pure implicit IRK2 scheme for grid with $256^{2}$ DOFs at same CFL.

The convection of a vortex on an irregular grid of $260^{2}$ DOFs is also performed. The time integration is performed according to CFL condition of the biggest cells of the domain such as

[TABLE]

In this irregular domain, the size ratio between the largest and the smallest cell is equal to $11$. Hence, the constant time step of computation will be determined according to the CFL number of the largest cell, which leads to the CFL number $11$ time higher for the smallest cells than for the larger ones. The computation is performed at several CFL values, during one time period, in order to compare stability properties of the several time integrators. The parameter $\omega_{j}$ is configured such as

[TABLE]

The parameter $\omega_{j}$ is defined according to the size of cells in the domain (see Fig. 44 with grey zone corresponding to $\omega_{j}=1$).

In this case the proportion of cells with $\omega_{j}<1$ corresponding to hybrid and implicit cells is equal to $73\%$. The final computational times for the several time integrator and CFL numbers is provided in Tab. 7. The values are normalized with respect to computational cost of the Heun time integrator at CFL=0.1. According to Tab. 7, as expected, the explicit time integrator is unstable for high CFL number due to small cells. Then, the AION scheme, as expected, tends to have better stability properties than full Heun time integrator, due to implicit time treatment in the smallest cells. The computational gain compared to full implicit computation is equal to $7\%$, due to high percentage of hybrid and implicit cells.

In the following a focus is made on testing the hybrid treatment of the gradient for the diffusion scheme.

7.4 Three-dimensional Taylor Green Vortex

A three-dimensional laminar flow computation is now performed to assess the performance of the AION integration scheme for a viscous test case. The proposed test case is the Taylor Green vortex at Reynolds number ($Re=\frac{\rho_{\infty}U_{\infty}L}{\mu}$) equal to $1600$ and Mach number equal to 0.1. The initial solution is defined inside the periodic cubic domain $[-\pi L,\pi L]^{3}$ ($L=1$):

[TABLE]

This test case is defined in order to observe a kind of vortex cascade, as encountered in turbulent flow computations. For the computation, the fluid is assumed to be a compressible perfect gas with $\gamma=1.4$ and the Prandtl number ($Pr=\frac{\mu.c_{p}}{\kappa}$) is equal to $0.71$, with $c_{p}$ and $\kappa$, respectively, the heat capacities at constant pressure and the heat conductivity. The initial temperature is kept constant $T_{\infty}=\frac{P_{\infty}}{R_{gas}\rho_{\infty}}$ with $R_{gas}$ the perfect gas constant. From the local pressure and temperature is deduced the density $\rho=\frac{p}{R_{gas}T_{\infty}}$.

The simulation is performed until the physical time of $t=20t_{c}$ with $t_{c}=\frac{L}{U_{\infty}}$ the characteristic convective time is attained. The computation is performed with a 2-exact formulation for the spatial discretisation (three order accurate spatially) and time-integrated by Heun and AION schemes on a $256^{3}$ grid at $CFL=0.9$. The following outputs are performed:

• •

The temporal evolution of the kinetic energy integrated over the whole domain:

[TABLE]

• •

The temporal evolution of the kinetic energy dissipation rate:

[TABLE]

• •

The temporal evolution of the kinetic energy dissipation rate on the whole domain:

[TABLE]

with $w$ the vorticity.

These outputs are compared with the DNS results obtain on a $512^{3}$ grid with pseudo-spectral method, reference for the First International Workshop on High-Order CFD Methods held at the 50th AIAA Aerospace Meeting (https://www.grc.nasa.gov/hiocfd/). It appears that the AION scheme seems to be quite performing for this kind of viscous case. Indeed the AION scheme fits better the DNS results than Heun scheme for time range $t=[0t_{c},8t_{c}]$ and $t=[12t_{c},20t_{c}]$ according to the $E_{k}$ and $\varepsilon$ (see Fig. 45 - 46) when Heun scheme seems to be more dissipative in these time ranges. Nevertheless the Heun scheme seems to better reach the maximum of $\varepsilon$ than the AION scheme on the time range $t=[8t_{c},12t_{c}]$. On the contrary, the AION scheme gives better result for the enstrophy $E$, on the whole time range than Heun scheme (see Fig. 47). This results confirm the particularity that the Heun time integrator has a significant importance on accuracy of the numerical results. Indeed this conclusion was already observed on the slope of $log(E_{total})$ obtained in the previous case of the convected vortex. This comparison allows us to validate our hybridation of the viscous gradients presented previously.

8 Conclusion

RANS and LES equations have almost the same shape and their own pros and cons. In an industrial environment, it could be useful to perform LES far from the wall and RANS near the wall in order to account easily for the turbulence effects, while keeping an acceptable mesh size. Standard time integration schemes for LES and RANS equations do not follow the same constraints. For LES, an explicit time integration is chosen in order to control simply and efficiently spectral properties (dissipation and dispersion). For RANS equations, it is of paramount importance to reach the steady state as fast as possible, using an implicit time integration procedure.

Indeed, the current paper specifically addresses how explicit and implicit time integrators can be coupled spatially. Here, two standard time integration schemes (Heun and second order implicit Runge-Kutta schemes -IRK2-) are hybridized / blended using a transition function $\omega$, while keeping the standard expected properties (spectral behaviour). In the first part of the paper, a way to couple the proposed schemes, adapted from the literature, was first introduced but lead to instability for some wavenumbers and CFL values. A new alternative approach, named AION scheme, was proposed and designed in order to correct the unexpected behaviour of the first coupling procedure. The spectral analysis performed on the coupled space / AION schemes enabled us to check the stability of the coupling procedure. In addition, in order to minimize the CPU cost, transition to explicit and implicit schemes was performed by reducing the transition area, playing with the values of $\omega$. After the spectral analysis, attention was paid on simulations of increasing complexity, from the 1D Gaussian hump advection to 3D Taylor Green Vortex. Three aspects are of paramount importnace. First, the results with the new scheme were shown to have same or better quality than the standard basic schemes. Moreover, starting from the reference full-implicit time integration, the hybrid formulation enables to reduce the CPU cost. Finally, the AION formulation seems less dependent on the Newton’s algorithm convergence than the standard IRK2 scheme. Moreover, the convergence of the IRK2 in a fully unsteady simulation needs many steps, something which was not necessary with the proposed AION scheme.

We are working today on two new improvements. The first one deals with the use of an adaptive time-step procedure. In this case, explicit cells are separated into several classes depending of their local maximum stable time step. All the cells inside the same class are time-integrated using the same time step and for two adjacent classes, time step for the smallest cells is half the one for the largest cells. Attention must be paid on the computation of the flux for faces having the cell of one side in a given class and the cell on the other side in another one. This work is ongoing and will be summarized in a new paper. The last improvement concerns the application of the proposed technique to a coupled RANS/LES computation for demonstrating the capability in an industrial environment. We can assume that the corner stone will be the scaling of the hybrid parameter $\omega$ with the function that switches between RANS and LES models.