Efficient high order accurate staggered semi-implicit discontinuous Galerkin methods for natural convection problems

TL;DR

This paper introduces high-order staggered semi-implicit discontinuous Galerkin methods for natural convection, combining explicit and semi-Lagrangian schemes to improve efficiency and stability in simulating buoyancy-driven flows.

Contribution

It develops a novel class of high-order staggered semi-implicit DG schemes for natural convection, extending previous work to include buoyancy effects and compressible flows.

Findings

Schemes show good agreement with benchmark data.

Semi-Lagrangian approach enhances computational efficiency.

Method effectively handles both incompressible and compressible flows.

Abstract

We propose a new family of high order staggered semi-implicit discontinuous Galerkin (DG) methods for the simulation of natural convection problems. Assuming small temperature fluctuations, the Boussinesq approximation is valid and the flow can simply be modeled by the incompressible Navier-Stokes equations coupled with a transport equation for the temperature and a buoyancy source term in the momentum equation. Our numerical scheme is developed starting from the work presented in [TD14], in which the spatial domain is discretized using a face-based staggered unstructured mesh. For the computation of the advection and diffusion terms, two different algorithms are presented: i) a purely Eulerian explicit upwind-type scheme and ii) a semi-Lagrangian approach. The first methodology leads to a conservative scheme whose major drawback is the time step restriction imposed by the CFL stability condition. On the contrary, computational efficiency can be notably improved relying on a semi-Lagrangian approach. This method leads to an unconditionally stable scheme if the diffusive terms are discretized implicitly. Once the advection and diffusion contributions have been computed, the pressure Poisson equation is solved and the velocity is updated. As a second model for the computation of buoyancy-driven flows, we also consider the full compressible Navier-Stokes equations. The staggered semi-implicit DG method first proposed in [TD17] for all Mach number flows is properly extended to account for the gravity source terms arising in the momentum and energy conservation laws. The validity and the robustness of our novel class of staggered semi-implicit DG schemes is assessed at the aid of several classical benchmark problems, showing in all cases a good agreement with available numerical reference data. Finally, a detailed comparison between the incompressible and the compressible solver is presented.

Tables

Table 1. Table 1: Taylor-Green vortex. Numerical convergence results for the velocity vector field.

$N_e$	$\nu =0$				$\nu =0.1$
	$(p,p_{\gamma })=(1,1)$		$(p,p_{\gamma })=(2,2)$		$(p,p_{\gamma })=(1,1)$		$(p,p_{\gamma })=(2,2)$
	$ϵ\left(𝐯\right)$	$𝒪$	$ϵ\left(𝐯\right)$	$𝒪$	$ϵ\left(𝐯\right)$	$𝒪$	$ϵ\left(𝐯\right)$	$𝒪$
$44$	$2.81E-01$		$6.72E-02$		$2.62E-01$		$4.45E-02$
$176$	$9.56E-02$	$1.6$	$1.33E-02$	$2.3$	$7.43E-02$	$1.8$	$5.02E-03$	$3.1$
$396$	$5.32E-02$	$1.4$	$5.25E-03$	$2.3$	$3.34E-02$	$2.0$	$1.43E-03$	$3.1$
$704$	$3.54E-02$	$1.4$	$2.65E-03$	$2.4$	$1.86E-02$	$2.0$	$5.91E-04$	$3.1$
$1100$	$2.59E-02$	$1.4$	$1.54E-03$	$2.4$	$1.18E-02$	$2.1$	$2.97E-04$	$3.1$
$1584$	$1.99E-02$	$1.4$	$9.80E-04$	$2.5$	$8.15E-03$	$2.0$	$1.70E-04$	$3.1$
$2156$	$1.59E-02$	$1.5$	$6.66E-04$	$2.5$	$5.96E-03$	$2.0$	$1.06E-04$	$3.0$

Table 2. Table 2: Differentially heated cavity test. Operating conditions.

Rayleigh	$𝐠$	${\rho }_0$	$\alpha$	$\beta$	$\nu$
${10}^3$	${(0,-9.81)}^T$	$1$	$6.836\cdot {10}^{-3}$	$3.4112\cdot {10}^{-3}$	$4.874\cdot {10}^{-3}$
${10}^4$			$2.171\cdot {10}^{-3}$		$1.541\cdot {10}^{-3}$
${10}^5$			$6.836\cdot {10}^{-4}$		$4.874\cdot {10}^{-4}$
${10}^6$			$2.171\cdot {10}^{-4}$		$1.541\cdot {10}^{-4}$
${10}^7$			$6.836\cdot {10}^{-5}$		$4.874\cdot {10}^{-5}$

Table 3. Table 3: Differentially heated cavity (DHC) test. Comparison of the Nusselt number obtained with the Eulerian (Eu) and Eulerian-Lagrangian (EL) advection scheme with available reference solutions from the literature.

Ra	STIN2D Eu.		STIN2D EL		Ref. [96]	Ref. [97]	Ref. [98]	Ref. [99]	Ref. [6]
	${\Gamma }_h$	${\Gamma }_c$	${\Gamma }_h$	${\Gamma }_c$	Average	Average	Average	Average	Average
${10}^3$	$1.1199$	$1.1199$	$1.1212$	$1.1212$	$1.118$	$1.117$	$1.1149$	$1.117$	$1.112$
${10}^4$	$2.2471$	$2.2471$	$2.2486$	$2.2486$	$2.243$	$2.243$	$2.2593$	$2.254$	$2.198$
${10}^5$	$4.5283$	$4.5302$	$4.6605$	$4.5302$	$4.519$	$4.521$	$4.4832$	$4.598$	$4.465$
${10}^6$	$8.8651$	$8.8655$	$8.8592$	$8.8595$	$8.800$	$8.806$	$8.8811$	$8.976$	$8.783$
${10}^7$	$16.8394$	$16.8428$	$16.7684$	$16.7714$	$--$	$16.40$	$16.3869$	$16.656$	$16.46$

Table 4. Table 4: Differentially heated cavity (DHC) test. Comparison of the normalized maximum vertical velocity at y = 0.5 𝑦 0.5 y=0.5 obtained with the Eulerian (Eu) and Eulerian-Lagrangian (EL) advection schemes with available reference solutions from the literature.

Ra	STIN2D Eu.	STIN2D EL	Ref. [96]	Ref. [97]	Ref. [98]	Ref. [99]
${10}^3$	$3.7227$	$3.7324$	$3.679$	$3.692$	$3.6962$	$3.686$
${10}^4$	$19.6342$	$19.6372$	$19.51$	$19.63$	$19.6177$	$19.79$
${10}^5$	$68.9334$	$68.9526$	$68.22$	$68.85$	$68.6920$	$70.63$
${10}^6$	$220.603$	$220.212$	$216.75$	$221.6$	$220.8331$	$227.11$
${10}^7$	$701.051$	$700.679$	$--$	$702.3$	$703.2526$	$714.48$

Table 5. Table 5: Differentially heated cavity (DHC) test. Comparison of the normalized maximum horizontal velocity at x = 0.5 𝑥 0.5 x=0.5 obtained with the Eulerian (Eu) and Eulerian-Lagrangian (EL) advection schemes with available reference solutions from the literature.

Ra	STIN2D Eu.	STIN2D EL	Ref. [96]	Ref. [98]	Ref. [99]
${10}^3$	$3.6736$	$3.6828$	$3.634$	$3.6493$	$3.489$
${10}^4$	$16.1856$	$16.1877$	$16.2$	$16.1798$	$16.12$
${10}^5$	$34.8869$	$34.896$	$34.81$	$34.7741$	$33.39$
${10}^6$	$64.8511$	$64.9075$	$65.33$	$64.6912$	$65.40$
${10}^7$	$148.681$	$149.561$	$--$	$145.2666$	$143.56$

Equations

$$\beta ({\theta }_{\mathrm{b}}-{\theta }_0)\ll 1,$$

$$\nabla \cdot \mathbf{v}=0,$$

$$\frac{\partial \mathbf{v}}{\partial t}+\nabla \cdot {\mathbf{F}}_c^{\mathbf{v}}+\nabla p=\nabla \cdot \left(\nu \nabla \mathbf{v}\right)+\left(1-\beta \delta \theta \right)\mathbf{g},$$

$$\frac{\partial \theta }{\partial t}+\nabla \cdot {\mathbf{F}}_c^{\theta }=\nabla \cdot \left(\alpha \nabla \theta \right),$$

$${\mathbf{F}}^{\mathbf{v}}:={\mathbf{F}}_c^{\mathbf{v}}-\nu \nabla \mathbf{v},{\mathbf{L}}^{\mathbf{v}}\left(\mathbf{v}\right):=\frac{\partial }{\partial t}\mathbf{v}+\nabla \cdot {\mathbf{F}}^{\mathbf{v}},$$

$${\mathbf{F}}^{\theta }:={\mathbf{F}}_c^{\theta }-\nabla \cdot \alpha \nabla \theta ,{\mathbf{L}}^{\theta }\left(\theta ,\mathbf{v}\right):=\frac{\partial }{\partial t}\theta +\nabla \cdot {\mathbf{F}}^{\theta }.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot {\mathbf{F}}_c^{\rho \mathbf{v}}+\nabla P=\nabla \cdot \mathbf{\sigma }+\rho \mathbf{g},$$

$$\frac{\partial \rho E}{\partial t}+\nabla \cdot {\mathbf{F}}_c^{\rho E}=\nabla \cdot \left(\mathbf{\sigma }\mathbf{v}+\kappa \nabla \theta \right)+\rho \mathbf{g}\cdot \mathbf{v},$$

$${\mathbf{F}}^{\rho \mathbf{v}}:={\mathbf{F}}_c^{\rho \mathbf{v}}-\mathbf{\sigma },{\mathbf{L}}^{\rho \mathbf{v}}\left(\rho \mathbf{v}\right):=\frac{\partial }{\partial t}\mathbf{v}+\nabla \cdot {\mathbf{F}}^{\rho \mathbf{v}}.$$

$$P=\rho R\theta ,e=c_v\theta ,\text{hence}P=(\gamma -1)\rho e,$$

$$p_i(\mathbf{x},t)$$

$${\mathbf{v}}_j(\mathbf{x},t)$$

$${\theta }_i(\mathbf{x},t)$$

$${\int }_{{\mathbf{T}}_i}{\varphi }_{\kappa }^{(i)}\nabla \cdot \mathbf{v}d\mathbf{x}=0,$$

$${\int }_{{\mathbf{R}}_j}{\psi }_k^{(j)}\left(\frac{\partial \mathbf{v}}{\partial t}+\nabla \cdot {\mathbf{F}}^{\mathbf{v}}\right)d\mathbf{x}+{\int }_{{\mathbf{R}}_j}{\psi }_k^{(j)}\nabla pd\mathbf{x}={\int }_{{\mathbf{R}}_j}{\psi }_k^{(j)}\left(1-\beta \delta \theta \right)\mathbf{g}d\mathbf{x},$$

$${\int }_{{\mathbf{T}}_i}{\varphi }_{\kappa }^{(i)}\frac{\partial \theta }{\partial t}d\mathbf{x}+{\int }_{{\mathbf{T}}_i}{\varphi }_{\kappa }^{(i)}\nabla \cdot {\mathbf{F}}^{\theta }d\mathbf{x}=0.$$

$${\oint }_{\partial {\mathbf{T}}_i}{\varphi }_{\kappa }^{(i)}\mathbf{v}\cdot {\vec{n}}_{i}dS-{\int }_{{\mathbf{T}}_i}\nabla {\varphi }_{\kappa }^{(i)}\cdot \mathbf{v}d\mathbf{x}=0,$$

$$\sum _{j\in S_i}\left({\int }_{{\Gamma }_j}{\varphi }_{\kappa }^{(i)}{\mathbf{v}}_j\cdot {\vec{n}}_{ij}dS-{\int }_{{\mathbf{T}}_{i,j}}\nabla {\varphi }_{\kappa }^{(i)}\cdot {\mathbf{v}}_{j}d\mathbf{x}\right)=0,$$

$${\int }_{{\mathbf{R}}_j}{\psi }_k^{(j)}\left(\frac{\partial {\mathbf{v}}_j}{\partial t}+\nabla \cdot {\mathbf{F}}_j^{\mathbf{v}}\right)d\mathbf{x}+{\int }_{{\mathbf{T}}_{\ell (j),j}}{\psi }_k^{(j)}\nabla p_{\ell (j)}d\mathbf{x}+{\int }_{{\mathbf{T}}_{r(j),j}}{\psi }_k^{(j)}\nabla p_{r(j)}d\mathbf{x}+{\int }_{{\Gamma }_j}{\varphi }_k^{(j)}\left(p_{r(j)}-p_{\ell (j)}\right){\vec{n}}_{j}dS$$

$$={\int }_{{\mathbf{R}}_j}{\psi }_k^{(j)}\mathbf{g}d\mathbf{x}-{\int }_{{\mathbf{T}}_{\ell (j),j}}{\psi }_k^{(j)}\beta \delta {\theta }_{\ell (j)}\mathbf{g}d\mathbf{x}-{\int }_{{\mathbf{T}}_{r(j),j}}{\psi }_k^{(j)}\beta \delta {\theta }_{r(j)}\mathbf{g}d\mathbf{x},$$

$${\int }_{{\mathbf{T}}_i}{\varphi }_{\kappa }^{(i)}\frac{\partial {\theta }_i}{\partial t}d\mathbf{x}+{\int }_{{\mathbf{T}}_i}{\varphi }_{\kappa }^{(i)}\nabla \cdot {\mathbf{F}}_i^{\theta }d\mathbf{x}=0,$$

$$\sum _{j\in S_i}\left({\int }_{{\Gamma }_j}{\varphi }_{\kappa }^{(i)}{\varphi }_l^{(j)}{\vec{n}}_{ij}dS\cdot {\hat{\mathbf{v}}}_{l,j}^{n+1}-{\int }_{{\mathbf{T}}_{i,j}}\nabla {\varphi }_{\kappa }^{(i)}{\varphi }_l^{(j)}d\mathbf{x}\cdot {\hat{\mathbf{v}}}_{l,j}^{n+1}\right)=0,$$

$${\int }_{{\mathbf{R}}_j}{\psi }_k^{(j)}{\psi }_l^{(j)}d\mathbf{x}{\mathbf{L}}_h^{\mathbf{v}}({\hat{\mathbf{v}}}_{l,j})+{\int }_{{\mathbf{T}}_{\ell (j),j}}{\psi }_k^{(j)}\nabla {\varphi }_l^{(\ell (j))}d\mathbf{x}{\bar{p}}_{l,\ell (j)}+{\int }_{{\mathbf{T}}_{r(j),j}}{\psi }_k^{(j)}\nabla {\varphi }_l^{(r(j))}d\mathbf{x}{\bar{p}}_{l,r(j)}+{\int }_{{\Gamma }_j}{\psi }_k^{(j)}{\varphi }_l^{(r(j))}{\vec{n}}_{j}dS{\bar{p}}_{l,r(j)}$$

$$-{\int }_{{\Gamma }_j}{\psi }_k^{(j)}{\varphi }_l^{(\ell (j))}{\vec{n}}_{j}dS{\bar{p}}_{l,\ell (j)}={\int }_{{\mathbf{R}}_j}{\psi }_k^{(j)}\mathbf{g}d\mathbf{x}-\beta {\int }_{{\mathbf{T}}_{r(j),j}}{\psi }_k^{(j)}{\varphi }_l^{(r(j))}d\mathbf{x}\delta {\bar{\theta }}_{l,r(j)}\mathbf{g}-\beta {\int }_{{\mathbf{T}}_{\ell (j),j}}{\psi }_k^{(j)}{\varphi }_l^{(\ell (j))}d\mathbf{x}\delta {\bar{\theta }}_{l,\ell (j)}\mathbf{g},$$

$${\int }_{{\mathbf{T}}_i}{\varphi }_{\kappa }^{(i)}{\varphi }_l^{(i)}d\mathbf{x}{\mathbf{L}}_h^{\theta }({\bar{\theta }}_{l,i},{\bar{\mathbf{v}}}_{l,i})=0,$$

$$\sum _{j\in S_i}{\mathcal{D}}_{i,j}{\hat{\mathbf{v}}}_j=0,$$

$${\hat{\mathbf{M}}}_j{\mathbf{L}}_h^{\mathbf{v}}({\hat{\mathbf{v}}}_j)+{\mathcal{R}}_j{\bar{\mathbf{p}}}_{r(j)}-{\mathcal{L}}_j{\bar{\mathbf{p}}}_{\ell (j)}={\mathcal{G}}_j-{\mathcal{R}}_j^{\theta }\delta {\bar{\mathbf{\theta }}}_{r(j)}-{\mathcal{L}}_j^{\theta }\delta {\bar{\mathbf{\theta }}}_{\ell (j)},$$

$${\bar{\mathbf{M}}}_i{\mathbf{L}}_h^{\theta }({\bar{\mathbf{\theta }}}_i,{\bar{\mathbf{v}}}_i)=0,$$

$${\mathcal{D}}_{i,j}$$

$${\hat{\mathbf{M}}}_j$$

$${\mathcal{R}}_j$$

$${\mathcal{L}}_j$$

$${\mathcal{G}}_j$$

$${\mathcal{R}}_j^{\theta }$$

$${\mathcal{L}}_j^{\theta }$$

$${\bar{\mathbf{M}}}_i$$

$${\mathcal{Q}}_{i,j}={\int }_{{\mathbf{T}}_{i,j}}{\psi }_k^{(j)}\nabla {\varphi }_l^{(i)}d\mathbf{x}-{\int }_{{\Gamma }_j}{\psi }_k^{(j)}{\varphi }_l^{(i)}{\sigma }_{i,j}{\vec{n}}_{j}dS,$$

$${\sigma }_{i,j}=\frac{r(j)-2i+\ell (j)}{r(j)-\ell (j)}.$$

$${\hat{\mathbf{M}}}_j{\mathbf{L}}_h^{\mathbf{v}}({\hat{\mathbf{v}}}_j)+{\mathcal{Q}}_{r(j),j}{\bar{\mathbf{p}}}_{r(j)}+{\mathcal{Q}}_{\ell (j),j}{\bar{\mathbf{p}}}_{\ell (j)}={\mathcal{G}}_j-{\mathcal{R}}_j^{\theta }\delta {\bar{\mathbf{\theta }}}_{r(j)}-{\mathcal{L}}_j^{\theta }\delta {\bar{\mathbf{\theta }}}_{\ell (j)}.$$

$$\sum _{j\in S_i}{\mathcal{D}}_{i,j}{\hat{\mathbf{v}}}_j^{n+1}=0,$$

$$\displaystyle\frac{1}{\Delta t}\hat{\bm{M}}_{j}\left(\hat{\mathbf{v}}_{j}^{n+1}-\mathbf{F}_{t}^{\mathbf{v}}\left(\mathbf{v}^{n}\right)\right)+\bm{\mathcal{Q}}_{r(j),j}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\bar{\mathbf{p}}_{r(j)}^{n+\Theta}}+\bm{\mathcal{Q}}_{\ell(j),j}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\bar{\mathbf{p}}_{\ell(j)}^{n+\Theta}}=\bm{\mathcal{G}}_{j}-\bm{\mathcal{R}}_{j}^{\theta}\,\delta\bar{\bm{\theta}}_{r(j)}^{n+1}-\bm{\mathcal{L}}_{j}^{\theta}\,\delta\bar{\bm{\theta}}_{\ell(j)}^{n+1},$$

$${\bar{\mathbf{M}}}_i{\bar{\mathbf{\theta }}}_i^{n+1}-{\bar{\mathbf{M}}}_i{\bar{\mathbf{F}}}_t^{\theta }\left({\mathbf{\theta }}^n,{\mathbf{v}}^n\right)=0,$$

Full text

††journal: Computers & fluids

1 Introduction

Natural convection problems play an important role in computational fluid dynamics. They appear in numerous engineering applications and natural phenomena ranging from the design of cooling devices in industrial processes, electronics, building isolation or solar energy collectors, to the simulation of atmospheric flows. In the last decades, the scientific community has put a lot of efforts into the study of these phenomena, see e.g. [5, 6, 7, 8, 9, 10] for a non-exhaustive overview. Nowadays, the main challenge is to develop efficient high order numerical methods which are able to capture even small scale structures of the flow, avoiding the use of RANS turbulence models (see [11, 12]). In this paper, we propose a novel family of high order accurate staggered semi-implicit discontinuous Galerkin (DG) methods, which extends the works presented in [1, 3, 4] appropriately to deal also with gravity driven flows.

Depending on the magnitude of the temperature perturbation and on the importance of density changes, natural convection problems are usually divided into two main groups. If the Mach number and the temperature fluctuations are small, the incompressible Navier-Stokes equations under the usual Boussinesq assumption can be applied. Otherwise, the full compressible Navier-Stokes equations have to be employed. In the following, we will mainly focus on the first case. However, in this paper also the compressible model will be considered, thus allowing for a direct comparison of the results obtained using the two different systems of governing partial differential equations. Therefore, we will be able to further validate the applicability of the Boussinesq approach for the flow regimes we are interested in.

In the literature there are numerous approaches that have been proposed for the solution of the Navier-Stokes equations, such as finite difference methods [13, 14, 15, 16] or continuous finite element schemes [17, 18, 19, 20, 21, 22, 23]. Nevertheless, the construction of high order numerical methods, and especially of high order discontinuous Galerkin (DG) finite element schemes, is still a very active research field, which has started with the pioneering works of Bassi and Rebay [24], and Baumann and Oden [25, 26]. Later, several high order DG methods for the incompressible and compressible Navier-Stokes equations have been proposed, see for example [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]. We also would like to mention recent works on semi-implicit DG schemes that can be found in [40, 41, 42, 43, 44], to which our approach is indirectly related.

The algorithm proposed in this article makes use of the novel family of staggered semi-implicit DG schemes that has been introduced in [1, 2, 3] for the incompressible Navier-Stokes equations in two and three space dimensions and which was later also extended to the full compressible regime in [4], following the work outlined in [45, 46, 47]. These arbitrary high order accurate DG schemes are constructed on staggered unstructured meshes. The pressure, the density and the energy are defined on the triangular or tetrahedral primal grid, whereas the velocity is computed on a face-based staggered dual mesh. While the use of staggered grids is a very common practice in the finite difference and finite volume framework (see e.g. [13, 14, 48]), its use is not so widespread in the context of high order DG schemes. The first staggered DG methods, which adopted a vertex-based dual grid, have been proposed in [49, 50]. Other recent high order staggered DG algorithms that rely on an edge-based dual grid have been advanced in [51, 52]. For high order staggered semi-implicit discontinuous Galerkin schemes on uniform and adaptive Cartesian meshes, see [53, 54].

Focusing on the incompressible model, we propose two different approaches for the computation of the nonlinear convective terms. On the one hand, we consider the methodology already introduced in [3]. There, the convective subsystem for the velocity is solved considering the Rusanov flux function for an explicit upwind-type discretization of the nonlinear convective terms. Instead, the viscous terms are discretized implicitly, making again use of the dual mesh in order to obtain the discrete gradients, without needing any additional numerical flux function for the viscous terms. One of the major drawbacks of this approach is its high computational cost coming from the small time step dictated by the CFL stability condition due to the explicit discretization of the convective terms. Moreover, to avoid spurious oscillations, a limiter should be used (see [4]). As an alternative option, which is at the same time able to deal with large gradients and substantially reduces the computational cost, we propose the use of an * Eulerian-Lagrangian* approach, recently forwarded in [55] also in the context of high order in space staggered DG schemes. The trajectory of the flow particles is followed backward in time by integrating the associated trajectory equations at the aid of a high order Taylor series expansion, where time derivatives are replaced by spatial derivatives using the Cauchy-Kovalevskaya procedure, similar to the ADER approach of Toro and Titarev [56, 57, 58]. The high order spatial discretization of the DG scheme is then employed to obtain a high order time integration for each point needed to solve numerically the advection part of the governing equations. For further information on efficient semi-Lagrangian and Eulerian-Lagrangian schemes we refer the reader to [59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69].

The use of the Boussinesq assumption yields the coupling of the incompressible Navier-Stokes equations with an additional conservation equation for the temperature. The computation of the related advection and diffusion terms is performed similarly to what is done for the velocity in the momentum equation. Nevertheless, let us remark that the temperature is defined on the primal mesh so that interpolation from one mesh to the other is avoided in the fully Eulerian approach. Once the new temperature is known, the gravity source term in the momentum equation can be evaluated. Finally, the pressure Poisson equation is solved and the velocity at the new time step is computed.

Regarding the compressible Navier-Stokes equations, we extend the numerical scheme introduced in [4] to consider the additional gravity terms. To this end, two new terms are included in the pressure system which has been obtained by formal substitution of the discrete momentum equation into the discrete energy equation. The first gravity term, coming from the momentum equation, is computed jointly with the convective and viscous terms of the momentum equation at the beginning of each time step. Meanwhile, the gravity term embedded in the energy equation is computed at each Picard iteration using the updated values of the linear momentum density.

The rest of the paper is organized as follows. In Section 2 we recall the incompressible and compressible Navier-Stokes equations. For the incompressible model, the Boussinesq assumption is made to account for fluid flow with small temperature variations under gravity effects. Concerning the compressible model, we consider the full Navier-Stokes equations including the conservation law for the total energy density and assuming here the equation of state for an ideal gas. Section 3 is devoted to the description of the semi-implicit staggered DG method used to solve the incompressible model in two and three space dimensions. We start by recalling some basic definitions about the usage of staggered meshes and the polynomial spaces which are employed. Next, we derive the numerical method considering two different frameworks for the discretization of convective and diffusive terms, namely an Eulerian and an Eulerian-Lagrangian approach. The extension of the algorithm to the compressible case is described in Section 4. Several benchmarks are presented in Section 6, aiming at assessing the validity, efficiency and the robustness of our novel numerical schemes. The main pros and cons of the Eulerian and the Eulerian-Lagrangian approaches are analyzed as well. Finally, we compare the results obtained with the incompressible solver against those computed with the compressible solver in the low Mach number regime.

2 Governing equations

As already mentioned, natural convection problems may be studied using two different models: the incompressible and the compressible Navier-Stokes equations, both including proper gravitational terms. The choice of the model usually depends on specific features of the flow under consideration, like the magnitude of the temperature fluctuations or the importance of capturing density variations. Here, we are mainly interested in small temperature changes, so that we will firt focus on the incompressible case. Later, even the full compressible model will be studied and numerical results will be compared considering both approaches.

2.1 Incompressible Navier-Stokes equations

Let us consider the laminar flow of a single phase Newtonian fluid without neither radiation nor chemical reactions. Let furthermore $\beta$ be the thermal expansion coefficient of the fluid, $\theta_{0}$ the reference temperature of the flow and $\theta_{\mathrm{b}}-\theta_{0}$ the maximum temperature fluctuation. Under the assumption

[TABLE]

the Boussinesq approximation for buoyancy-driven flows holds. Therefore, natural convection problems with small temperature gradients may be analyzed by solving the system of incompressible Navier-Stokes equations coupled with an energy conservation equation through an additional source term in the momentum equation. The governing PDE system reads as follows:

[TABLE]

where $\mathbf{v}$ is the velocity field; $p=P/\rho$ indicates the normalized fluid pressure; $P$ is the physical pressure and $\rho$ is the fluid density, which according to the Boussinesq approximation is assumed to be constant everywhere apart from the gravity source term in the momentum equation; $\nu=\mu/\rho$ is the kinematic viscosity coefficient; $\mathbf{F}^{\mathbf{v}}_{c}=\mathbf{v}\otimes\mathbf{v}$ is the flux tensor of the nonlinear convective terms; $\delta\theta:=\theta-\theta_{0}$ is the temperature difference; $\mathbf{g}$ is the gravity acceleration; $\theta$ is the temperature; $\mathbf{F}_{c}^{\theta}=\mathbf{v}\theta$ is the flux tensor of the nonlinear convective terms of the energy equation; $\alpha=\kappa/(\rho c_{P})$ is the thermal diffusivity, which depends on the thermal conductivity $\kappa$, the density $\rho$ and the heat capacity $c_{P}$. Let us also introduce the following notation:

[TABLE]

2.2 Compressible Navier-Stokes equations

Large temperature fluctuations may produce substantial changes in the density of a fluid. As a consequence, the incompressible model (2)-(4) is no longer valid for large temperature fluctuations. In this case, one must use the full compressible Navier-Stokes equations under gravitational effects, that read

[TABLE]

where $\mathbf{F}^{\rho\mathbf{v}}_{c}=\rho\mathbf{v}\otimes\mathbf{v}$ is the convective flux for the momentum equation; $\bm{\sigma}=\mu(\nabla{\mathbf{v}}+\nabla{\mathbf{v}}^{\top})-\frac{2}{3}\mu\left(\nabla\cdot{\mathbf{v}}\right)\mathbf{I}$ is the viscous stress tensor; $\rho E=\rho{e}+\rho{k}$ is the total energy density; $\rho{k}=\frac{1}{2}\rho{\mathbf{v}}^{2}$ is the kinetic energy density; $e=e(P,\rho)$ represents the specific internal energy per unit mass and is given by the equation of state (EOS) as a function of the pressure $P$ and the density $\rho$; $H=e+\frac{P}{\rho}$ denotes the specific enthalpy; $\mathbf{F}^{\rho E}_{c}=\mathbf{v}\left(\rho E+P\right)=\rho\mathbf{v}k+\rho\mathbf{v}H$ is the convective flux for the energy conservation equation; $\kappa$ is the thermal conductivity coefficient. Let us also define the following operators:

[TABLE]

In our approach we assume that we are dealing with an ideal gas, so that the thermal and caloric equation of state (EOS) that are needed to close the above system are given by

[TABLE]

with $R=c_{P}+c_{v}$ denoting the specific gas constant; $c_{P}$ and $c_{v}$ representing the heat capacities at constant pressure and at constant volume, respectively, and $\gamma=c_{P}/c_{v}$ being the usual ratio of specific heats.

3 Numerical method for the incompressible model

System (2)-(4) will be solved starting by the staggered semi-implicit discontinuous Galerkin scheme detailed in [1, 2, 3]. Here, we recall the main ingredients of the algorithm, while for an exhaustive description the reader is referred to the aforementioned references.

3.1 Staggered unstructured mesh

The computational domain is discretized using a face-based staggered unstructured meshes, as adopted in [1, 2, 3, 4]. In what follows, we briefly summarize the grid construction and the main notation for the two dimensional triangular grid. After that, the primal and dual spatial elements are extended to the three dimensional case.

The spatial computational domain $\Omega\subset\mathbb{R}^{2}$ is covered with a set of $N_{e}$ non-overlapping triangular elements $\mathbf{T}_{i}$ with $i=1\ldots N_{e}$. By denoting with $N_{d}$ the total number of edges, the $j-$th edge will be called $\Gamma_{j}$. $\mathcal{B}(\Omega)$ refers to the set of indices $j$ corresponding to boundary edges. The three edges of each triangle $\mathbf{T}_{i}$ constitute the set $S_{i}$ defined by $S_{i}=\{j\in[1,N_{d}]\,\,|\,\,\Gamma_{j}\mbox{ is an edge of }\mathbf{T}_{i}\}$. For every $j\in[1\ldots N_{d}]-\mathcal{B}(\Omega)$ there exist two triangles $i_{1}$ and $i_{2}$ that share $\Gamma_{j}$. We assign arbitrarily a left and a right triangle called $\ell(j)$ and $r(j)$, respectively. The standard positive direction is assumed to be from left to right. Let $\vec{n}_{j}$ denote the unit normal vector defined on the edge $j$ and oriented with respect to the positive direction. For every triangular element $i$ and edge $j\in S_{i}$, the neighbor triangle of element $\mathbf{T}_{i}$ that share the edge $\Gamma_{j}$ is denoted by $\mathbb{\wp}(i,j)$.

For every $j\in[1,N_{d}]-\mathcal{B}(\Omega)$ the quadrilateral element associated to $j$ is called $\mathbf{R}_{j}$ and it is defined, in general, by the two centers of gravity of $\ell(j)$ and $r(j)$ and the two terminal nodes of $\Gamma_{j}$, see also [70, 71, 72, 52]. We denote by $\mathbf{T}_{i,j}=\mathbf{R}_{j}\cap\mathbf{T}_{i}$ the intersection element for every $i$ and $j\in S_{i}$. Figure 1 summarizes the notation, the primal triangular mesh and the dual quadrilateral grid.

According to [2], we will call the main grid, or primal grid, the mesh of triangular elements $\{\mathbf{T}_{i}\}_{i\in[1,N_{e}]}$, whereas the quadrilateral grid $\{\mathbf{R}_{j}\}_{j\in[1,N_{d}]}$ is addressed as the dual grid.

The definitions given above are then readily extended to three space dimensions within a domain $\Omega\subset\mathbb{R}^{3}$.

An example of the resulting main and dual grids in three space dimensions is reported in Figure 2. The main grid consists of tetrahedral simplex elements, and the face-based dual elements contain the three vertices of the common triangular face of two tetrahedra (a left and a right one), and the two barycenters of the two tetrahedra that share the same face. Therefore, in three space dimensions the dual grid consists of non-standard five-point hexahedral elements. The same face-based staggered dual mesh has also been used in [71, 73, 74, 75].

3.2 Basis functions

The basis functions are chosen according to [1, 3], thus, in the two dimensional case we first construct the polynomial basis up to a generic polynomial of degree $p$ on some triangular and quadrilateral reference elements with local coordinates $\xi$ and $\eta$. The reference triangle is considered to be $T_{std}=\{(\xi,\eta)\in\mathbb{R}^{2}\,\,|\,\,0\leq\xi\leq 1,\,\,0\leq\eta\leq 1-\xi\}$ and the reference quadrilateral element is defined as $R_{std}=[0,1]^{2}$. Then, the standard nodal approach of conforming continuous finite elements yields $N_{\phi}=\frac{(p+1)(p+2)}{2}$ basis functions denoted with $\{\phi_{k}\}_{k\in[1,N_{\phi}]}$ on $T_{std}$, and $N_{\psi}=(p+1)^{2}$ basis functions referred to as $\{\psi_{k}\}_{k\in[1,N_{\psi}]}$ on $R_{std}$. The transformation between the reference coordinates $\bm{\xi}=(\xi,\eta)$ and the physical coordinates $\mathbf{x}=(x,y)$ is performed by the maps $T_{i}:\mathbf{T}_{i}\longrightarrow T_{std}$ for every $i=1\ldots N_{e}$ and $T_{j}:\mathbf{R}_{j}\longrightarrow R_{std}$ for every $j=1\ldots N_{d}$, with the associated inverse relations, namely $T_{i}^{-1}:\mathbf{T}_{i}\longleftarrow T_{std}$ and $T_{j}^{-1}:\mathbf{R}_{j}\longleftarrow R_{std}$, respectively.

For the three-dimensional tetrahedra, we use again the standard nodal basis functions of conforming finite elements based on the reference element $T_{std}=\{(\xi,\eta,\zeta)\in\mathbb{R}^{3}\,\,|\,\,0\leq\xi\leq 1,\,\,0\leq\eta\leq 1-\xi,\,\,0\leq\zeta\leq 1-\xi-\eta\}$ and then we define a map to connect the reference space, $\bm{\xi}=(\xi,\eta,\zeta)$, to the physical space, $\mathbf{x}=(x,y,z)$, and vice-versa. Unfortunately, the non-standard five-point hexahedral elements of the dual mesh entail the definition of a polynomial basis directly in the physical space using a simple modal basis based on rescaled Taylor monomials, such as the ones proposed in [3]. We thus obtain $N_{\phi}=N_{\psi}=\frac{1}{6}(p+1)(p+2)(p+3)$ basis functions per element for both the main grid and the dual mesh.

3.3 Staggered semi-implicit DG scheme

The discrete pressure and temperature are defined on the main grid, that is $p_{i}(\mathbf{x},t)=p_{h}(\mathbf{x},t)|_{\mathbf{T}_{i}}$ and $\theta_{i}(\mathbf{x},t)=\theta_{h}(\mathbf{x},t)|_{\mathbf{T}_{i}}$, while the discrete velocity is defined on the dual grid, namely ${\mathbf{v}}_{j}(\mathbf{x},t)={\mathbf{v}}_{h}(\mathbf{x},t)|_{\mathbf{R}_{j}}$.

Therefore, the numerical solution of (2)-(4) is given within each spatial element by

[TABLE]

In the rest of the paper, we use the convention that variables indexed by $j$ are defined on the dual grid, while the index $i$ is used for the quantities which refer to the main grid. Furthermore, we will use the hat symbol to denote degrees of freedom on the dual grid, whereas bars are used to distinguish degrees of freedom on the primal mesh. The set of variables on the main grid will be denoted by $\{Q_{i}\}_{i\in[1,N_{e}]}$, while $\{Q_{j}\}_{j\in[1,N_{d}]}$ corresponds to the unknowns defined on the dual grid. The vector of basis functions $\mathbf{\phi}^{(i)}(\mathbf{x})$ is generated via the map $T_{i}^{-1}$ from $\mathbf{\phi}(\bm{\xi})$ on $T_{std}$. The vector $\mathbf{\psi}^{(j)}(\mathbf{x})$ is generated from $\mathbf{\psi}(\bm{\xi})$ on $R_{std}$ through the mapping $T_{j}^{-1}$ in the two dimensional case, and it is directly defined in the physical space for each element in the three dimensional case, see [3].

Multiplying equations (2) and (4) by $\phi_{k}$, equation (3) by $\psi_{\kappa}$ and integrating on the related control volumes $\mathbf{T}_{i}$ and $\mathbf{R}_{j}$, respectively, we obtain the weak formulation of the incompressible model (2)-(4) for every $i=1\ldots N_{e}$, $\kappa=1\ldots N_{\phi}$, $j=1\ldots N_{d}$ and $k=1\ldots N_{\psi}$:

[TABLE]

Integration by parts of (15) leads to

[TABLE]

with $\vec{n}_{i}$ indicating the outward pointing unit normal vector. Next, taking into account the discontinuities of $p_{h}$ and $\theta_{h}$ on $\Gamma_{j}$ and $\mathbf{v}_{h}$ on the edges of $\mathbf{R}_{j}$, the weak form becomes

[TABLE]

where $\vec{n}_{ij}=\vec{n}_{i}|_{\Gamma_{j}}$. Note that the pressure has a discontinuity along $\Gamma_{j}$ inside the dual element $\mathbf{R}_{j}$ and hence the pressure gradient in (16) needs to be interpreted in the sense of distributions, as in path-conservative finite volume schemes [76, 77]. This leads to the jump terms present in (20), see [2]. Alternatively, the same jump term can be produced also via forward and backward integration by parts, see e.g. the well-known work of Bassi and Rebay [24]. Using definitions (12)-(14), we rewrite the above equations as

[TABLE]

where the Einstein summation convention over repeated indexes holds. Moreover, $\mathbf{L}_{h}^{\mathbf{v}}$ and $\mathbf{L}^{\theta}_{h}$ represent appropriate discretizations of the operators $\mathbf{L}^{\mathbf{v}}$ and $\mathbf{L}^{\theta}$, which will be defined later.

For every $i$ and $j$ equations (22)-(24) can be written in a compact matrix form as

[TABLE]

with the matrix definitions

[TABLE]

The action of the matrices $\bm{\mathcal{L}}_{j}$ and $\bm{\mathcal{R}}_{j}$ can be generalized by introducing the new matrix $\bm{\mathcal{Q}}_{i,j}$, defined as

[TABLE]

with the sign function $\sigma_{i,j}$ given by

[TABLE]

In this way, $\bm{\mathcal{Q}}_{\ell(j),j}=-\bm{\mathcal{L}}_{j}$ and $\bm{\mathcal{Q}}_{r(j),j}=\bm{\mathcal{R}}_{j}$, hence the momentum equation (26) writes

[TABLE]

Time discretization using the theta ($\Theta$) method leads to

[TABLE]

where $\bar{\mathbf{p}}^{n+\Theta}=\Theta\bar{\mathbf{p}}^{n+1}+(1-\Theta)\bar{\mathbf{p}}^{n}$, with $\Theta$ an implicitness factor to be taken in the range $\Theta\in\left[\frac{1}{2},1\right]$, and $\mathbf{F}_{t}^{\mathbf{v}}\left(\mathbf{v}^{n}\right)$ and $\bar{\mathbf{F}}_{t}^{\theta}\left({\bm{\theta}^{n}},{\mathbf{v}^{n}}\right)$ are proper discretizations of the convective and diffusive terms of the momentum and energy equations, respectively, which read

[TABLE]

Two different approaches are considered in this work to obtain a numerical approximation of the operators (42)-(43): i) a fully Eulerian discretization and ii) an Eulerian-Lagrangian method. The Eulerian scheme provides a fully conservative formulation, contrarily to the Eulerian-Lagrangian approach [55]. However, the latter scheme would be unconditionally stable, so that no CFL restrictions need to be considered for the time step, thus substantially reducing the computational cost of the algorithm.

3.4 Pressure system

Formal substitution of equation (40) into (39), i.e. making use of the Schur complement, leads to a linear system in which the pressure is the only unknown, that is a scalar quantity:

[TABLE]

The above system is symmetric and in general positive semi-definite, see [2, 4], so that an efficient conjugate gradient method can be applied in order to obtain the new pressure $\bar{\mathbf{p}}^{n+1}$. The previous equation (44) coupled with (40)-(41) constitute the system of equations to be solved.

3.5 Nonlinear advection and diffusion

In the framework of semi-implicit schemes [78, 79, 80, 81, 74, 47], the nonlinear convective terms are typically discretized explicitly, while the pressure terms are treated implicitly. For more details on how these numerical methods are related to flux-vector splitting schemes, see [46].

Following [3], we exploit the advantages of using staggered grids to develop a suitable discretization of the nonlinear advection and diffusion terms. Therefore, the velocity field is first interpolated from the dual grid to the main grid:

[TABLE]

with

[TABLE]

Next, the nonlinear convective terms can be easily discretized with a standard DG scheme on the main grid. Finally, the staggered mesh is used again in order to define the gradient of the velocity on the dual elements, which yields a very simple and sparse system for the discretization of the viscous terms. A similar procedure is employed to discretize the convective and viscous terms of the energy equation. Let us remark that for this particular equation, the temperature is already defined on the primal elements, thus interpolation between the two meshes is avoided. Let us define two auxiliary variables for the diffusion terms, namely the stress tensor $\bm{\sigma}$ and the heat flux $\mathbf{q}$ as

[TABLE]

Then, the convective and viscous subsystems of the momentum and energy equations read

[TABLE]

Defining the temperature and the velocity on the primal elements and the auxiliary variables (48) on the dual grid, we obtain a weak formulation for equations (49)-(50) and (51)-(52), that is

[TABLE]

Accounting for the discretization in time, the former systems are expressed in matrix notation as

[TABLE]

with

[TABLE]

$\bm{\Upsilon}_{i}^{\mathbf{v},c}$ and $\bm{\Upsilon}_{i}^{\theta,c}$ are standard DG discretizations of the nonlinear convective terms, and $\mathbf{v}^{-}$, $\mathbf{v}^{+}$, $\bm{\theta}^{-}$, and $\bm{\theta}^{+}$ refer to the boundary extrapolated values from within the cell and from the neighbors, respectively. Furthermore, we use the Rusanov flux [82] as approximate Riemann solver:

[TABLE]

where $s^{\mathbf{v}}_{\max}=\max\left(2|\bar{\mathbf{v}}^{+}\cdot\vec{n}_{i}|,2|\bar{\mathbf{v}}^{-}\cdot\vec{n}_{i}|\right)$ and $s_{\max}^{\theta}=\max\left(|\bar{\mathbf{v}}^{+}\cdot\vec{n}_{i}|,|\bar{\mathbf{v}}^{-}\cdot\vec{n}_{i}|\right)$ are the maximum eigenvalues of the convective operators $\mathbf{F}^{\mathbf{v}}_{c}$ and $\mathbf{F}^{\theta}_{c}$, respectively. Finally, substituting (58) into (57), and (60) into (59), we obtain

[TABLE]

In order to avoid the solution of a nonlinear system due to the presence of the nonlinear operator related to the convective terms, a fractional step scheme combined with an outer Picard iteration is used:

[TABLE]

The previous procedure constitutes a discretization of the nonlinear convective and viscous terms on the main grid, both for the momentum and the energy equations, where $\bm{\sigma}$ and $\mathbf{q}$ are computed on the face-based dual mesh. To recover the contribution of the convective and viscous terms in the dual grid, as required in (40), we perform the following projection:

[TABLE]

3.6 Eulerian-Lagrangian approach

Instead of applying the Eulerian advection scheme illustrated in Section 3.5 for the approximation of the convection and diffusion terms, $\mathbf{F}^{\mathbf{v}}_{t}\left(\mathbf{v}^{n+1,k+\frac{1}{2}}\right)$ and $\bar{\mathbf{F}}_{t}^{\theta}\left({\bm{\theta}^{n+1,k+\frac{1}{2}}},{\mathbf{v}^{n+1,k+\frac{1}{2}}}\right)$, we may take into account the Lagrangian trajectory of the flow particles and use an Eulerian-Lagrangian approach. Specifically, the departure point $x^{*}$ of the Lagrangian trajectory has to be determined in order to compute the corresponding value of the transferred quantities, namely velocity $\mathbf{v}^{*}$ and temperature $\bm{\theta}^{*}$. The Lagrangian trajectory is defined by the solution $\mathbf{x}(\tau)$ of the trajectory equation

[TABLE]

where $\mathbf{x}_{l}^{n+1}$ is the location of a generic quadrature point from which the backtracking of the trajectory is started. Furthermore, $\tau$ represents the rescaled time coordinate referred to the time step $\Delta t:=t^{n+1}-t^{n}$, and can be easily evaluated as $\tau=t-t^{n}$, while ${\mathbf{v}}(\mathbf{x}(\tau))$ is the local fluid velocity. The foot point of the characteristics, which is nothing but the sought departure point, is given by $\mathbf{x}^{*}=\mathbf{x}(\Delta t)$. In order to solve the system of ordinary differential equations (ODE) (70), we rely on the approach presented in [65, 66], hence using a high order Taylor method, which leads to the solution $\mathbf{x}_{r+1}$ at the new time $\tau_{r+1}$:

[TABLE]

Expansion (71) allows the scheme to be up to third order accurate in $\tau$. The index $r$ represents the iteration number if a sub-time stepping is chosen for the approximation of the time step interval $[0,\Delta t]$. High order time derivatives are then replaced by high order spatial derivatives using repeatedly the trajectory equation (70) via the Cauchy-Kovalevskaya procedure, which is also typical for the ADER approach of Toro and Titarev [56, 57, 58]. Thus, assuming that ${\mathbf{v}}={\mathbf{v}}(\mathbf{x})={\mathbf{v}}(\mathbf{x},t^{n})$ is frozen during one time step, one obtains

[TABLE]

In this work we follow the methodology detailed in [55], in which a staggered high order DG scheme has been used together with the Eulerian-Lagrangian technique previously illustrated. In the DG context, the Eulerian-Lagrangian scheme is used to compute the nonlinear convective terms inside a weak formulation of the governing PDE, therefore the starting points for the trajectory equation (70) are given by the Gaussian quadrature points in space within each control volume $T_{i}$. In order to guarantee sufficient accuracy, we integrate backward in time twice the minimum number of Gaussian points that ensures the formal order of accuracy of the scheme, hence performing over-integration. According to [55], the integration of the ODE (70) is carried out in a reference system with local coordinates $(\xi,\eta,\zeta)$, so that element and physical boundaries can be easily identified and treated. Finally, the high order spatial discretization of the DG scheme is employed to compute the values of the velocity and the temperature at the foot point $\mathbf{x}^{*}=\mathbf{x}(\Delta t)$ of each trajectory. Further details on the aforementioned methodology can be found in [55], whereas in [83] a study for determining the departure points of trajectories is proposed. If stiff problems are considered, in which the solution rapidly changes in time, a complete space-time method would be necessary to capture properly the flow trajectory and which will be subject of future work. However, for the applications considered in this work, that are concerned with natural convection problems, this is typically not the case, so that our Eulerian-Lagrangian algorithm can be adopted for the discretization of the convection and diffusion terms, $\mathbf{F}^{\mathbf{v}}_{t}\left(\mathbf{v}^{n+1,k+\frac{1}{2}}\right)$ and $\bar{\mathbf{F}}_{t}^{\theta}\left({\bm{\theta}^{n+1,k+\frac{1}{2}}},{\mathbf{v}^{n+1,k+\frac{1}{2}}}\right)$.

The semi-implicit scheme for the Navier-Stokes equations together with the Eulerian-Lagrangian approach becomes unconditionally stable for arbitrary high order of accuracy and thus allows large time steps. For a detailed analysis in the case of semi-implicit finite difference schemes, see [80].

3.7 Overall method

Given $\bar{\mathbf{p}}_{i}^{n}$, $\hat{\mathbf{v}}_{j}^{n}$ and $\bar{\bm{\theta}}_{i}^{n}$, the final algorithm reads as follows.

The temperature at the new iteration, $\bar{\bm{\theta}}^{n+1,k+\frac{1}{2}}_{i}$, is obtained solving (68) when the fully Eulerian method is selected. Otherwise, the Eulerian-Lagrangian approach determines the term $\bar{\mathbf{F}}^{\theta}_{t}$ in equation (41) from which the temperature is updated. 2. 2.

The nonlinear convective and viscous terms for the momentum equation, $\mathbf{F}_{t}^{\mathbf{v}}\left(\mathbf{v}^{n+1,k+\frac{1}{2}}\right)$, are computed. Let us remark that within this term we are accounting for the contribution of the pressure at the previous Picard iteration, see Eqn. (67). 3. 3.

The pressure, $\bar{\mathbf{p}}_{i}^{n+1,k+1}$, results from solving system (44) after substitution of $\mathbf{F}_{t}^{\mathbf{v}}\left(\mathbf{v}^{n+1,k+\frac{1}{2}}\right)$ and $\bar{\bm{\theta}}^{n+1,k+\frac{1}{2}}_{i}$:

[TABLE] 4. 4.

The velocity, $\hat{\mathbf{v}}_{j}^{n+1,k+1}$, is then updated from (40):

[TABLE]

Remark 1**.**

For the sake of simplicity, the former algorithm has been presented assuming the theta method is used for time discretization. As a consequence, the accuracy of the resulting scheme is of arbitrary order in space and only up to second order in time. Besides, the space-time extension for the Eulerian approach has been developed following [2, 3], hence obtaining arbitrary high order semi-implicit DG schemes at the aid of test and basis functions that depend on both space and time.

4 Numerical method for the compressible model

The staggered semi-implicit discontinuous Galerkin scheme described in the previous section for the incompressible model has been extended in [4] to solve the compressible Navier-Stokes equations at all Mach numbers. Concerning semi-implicit finite volume schemes for all and low Mach number flows, we refer the reader to [84, 85, 86, 87, 88].

To simulate natural convection problems, some modifications are needed in order to incorporate the gravitational terms. In what follows, we will provide the details associated to their inclusion. For well-balanced schemes for the compressible Euler equations with gravity source terms, see [89, 90, 91, 9, 92].

4.1 Staggered semi-implicit DG scheme

The computational domain is discretized as already explained in Section 3.1 and the basis functions are defined according to Section 3.2. Then, the discrete pressure $p_{h}$, the fluid density $\rho_{h}$ and the discrete total energy density $(\rho E)_{h}$ are computed on the main grid, while the discrete velocity vector field ${\mathbf{v}}_{h}$, the discrete momentum density $(\rho{\mathbf{v}})_{h}$, and the discrete specific enthalpy $H_{h}$ are defined on the dual grid.

The numerical solution of (7)-(9) at a given time $t$ is represented inside the control volumes of the primal and the dual grids by piecewise spatial polynomials. The discrete pressure is approximated using (12), while the total energy as well as the density on the main mesh and the momentum on the dual mesh read

[TABLE]

Similarly to what has been done in Section 3.3 for the incompressible system, a weak formulation of (7)-(9) is obtained by multiplication of the governing equations by appropriate test functions and integration over the associated control volumes:

[TABLE]

The above weak formulation of the governing PDE accounts for the discontinuities of pressure and momentum along the boundaries of the primal and dual cells, respectively. Let $\mathbf{w}_{j}:=\bm{\sigma}_{j}\mathbf{v}_{j}$ be the work of the stress tensor in the energy equation and let $\mathbf{q}_{j}:=\kappa\nabla\theta_{j}$ define the heat flux vector. Using the polynomial approximations (74)-(76) in the semi-discrete system (77)-(79) leads to

[TABLE]

Finally, discretization in time of the above system yields

[TABLE]

where the convective and diffusive terms

[TABLE]

are evaluated following the procedure introduced in Section 3.5. The density on the dual mesh, needed in (84), is recovered from its value on the primal mesh as

[TABLE]

By introducing

[TABLE]

and using the matrix definitions (28)-(35), the above system is written more compactly as

[TABLE]

4.2 Pressure system and Picard iteration

The pressure appearing in the momentum equation (91) is discretized implicitly as well as the momentum in the energy equation. Then, a pressure equation can be derived by formal substitution of the discrete momentum equation (91) into the discrete energy equation (92):

[TABLE]

Here, $(\overline{\bm{\rho e}})_{i}^{n+1}=\bm{\rho}\mathbf{e}\left(\bar{\mathbf{p}}_{i}^{n+1},\bar{\bm{\rho}}_{i}^{n+1}\right)$ is to be understood as a componentwise evaluation of the internal energy density $\rho e$ for each component of the input vectors $\bar{\mathbf{p}}_{i}^{n+1}$ and $\bar{\bm{\rho}}_{i}^{n+1}$. Due to the dependency of the enthalpy on the pressure, and due to the presence of the kinetic energy in the total energy, system (93) is highly nonlinear. In order to obtain a linear system for the pressure, a simple but very efficient Picard iteration is used, as suggested in [93, 94, 95, 47, 1, 3]. Consequently, some nonlinear terms are discretized at the previous Picard iteration, and thus become essentially explicit. Specifically, only pressure remains fully implicit. Therefore, at each Picard iteration, $m=1,\ldots,N_{Pic}$, we need to solve

[TABLE]

We typically use a total number of $N_{Pic}=2$ iterations. The above equations introduce an implicit discretization for pressure and momentum in the momentum equation and in the energy equation, respectively. Second order of accuracy in time can be easily achieved by taking $\bar{\mathbf{p}}_{i}^{n+\Theta,m+1}$ and $(\widehat{\bm{\rho\mathbf{v}}})_{j}^{n+\Theta,m+1}$ in the momentum and energy equations with $\Theta=0.5$, which yields a Crank-Nicolson discretization.

4.3 Overall method

Given $\bar{\bm{\rho}}_{i}^{n}$, $\bar{\mathbf{p}}_{i}^{n}$, $\widehat{\bm{\rho}\mathbf{v}}_{j}^{n}$, and $\overline{\bm{\rho}\bm{E}}_{i}^{n}$, the staggered semi-implicit DG algorithm for solving the compressible Navier-Stokes equations (7)-(9) can therefore be summarized as follows.

The density at the new time time step, $\bar{\bm{\rho}}_{i}^{n+1}$, is obtained using (90). This value is then projected onto the dual mesh, $\hat{\bm{\rho}}_{j}^{n+1}$, by means of (88). 2. 2.

System (94)-(95) is solved relying on a Picard iteration procedure in order to evaluate the pressure and the momentum at the new time step, $\bar{\mathbf{p}}_{i}^{n+1,m+1}$ and $\widehat{\bm{\rho}\mathbf{v}}_{j}^{n+1,m+1}$, respectively. 3. 3.

The updated total energy density, $\overline{\bm{\rho}\bm{E}}_{i}^{n+1}$, is computed from (96) once the Picard loop has finished.

5 Time step restriction

The maximum time step is restricted by a CFL-type condition based on the local flow velocity:

[TABLE]

with $\mathrm{CFL}<1/d$, $d$ the space dimension, $h_{\mathrm{min}}$ the smallest insphere diameter (in 3D) or incircle radius (in 2D) and $\mathbf{v}_{\mathrm{max}}$ is the maximum convective speed. If viscous terms are present, the eigenvalues of the viscous operator have to be considered as well (see [33, 4]).

As commented in Section 3.6, if we employ the Eulerian-Lagrangian approach, the scheme becomes unconditionally stable for inviscid fluids, so that the time restriction is no more determined by the above CFL condition and $CFL\geq 1$ may be chosen. However, the parallel version of the code requires a safety factor to be specified, in order to ensure that the Lagrangian trajectories never exit the MPI neighborhood of each region (see [55] for further details).

6 Numerical test problems

In this section, classical benchmarks for natural convection problems are used in order to verify the validity and the efficiency of the novel algorithms presented in this work. Moreover, these tests allow us to analyze the strengths and drawbacks of our numerical schemes.

6.1 Taylor-Green Vortex with gravity

To analyse the accuracy of the proposed schemes, we consider a modification of the Taylor-Green vortex benchmark by including the gravity term in the momentum equation. The exact solution of this test case reads

[TABLE]

with $g=-9.81$. In order to satisfy the governing equations (2)-(3) with the definitions (98)-(100), the following source terms need to be added to the right hand side of the momentum equations:

[TABLE]

The simulations are run on the computational domain $\Omega=[0,2\pi]^{2}$ with periodic boundaries on a sequence of successively refined unstructured grids. Two settings are considered with different viscosity coefficient, namely $\nu=0$ and $\nu=0.1$. The convergence results at $t_{\mathrm{end}}=0.1$ are shown in Table 1. We observe that the optimal convergence rates are achieved for this non-trivial test with gravity and viscosity terms in a transient regime. The space-time DG discretization of [2, 3] has been employed with $p$ and $p_{\gamma}$ denoting the polynomial degree in space and time, respectively. From the obtained results we can conclude that the scheme converges with the expected convergence rate of at least $p+\frac{1}{2}$.

6.2 Differentially heated cavity

The differentially heated cavity test (DHC) has been proposed in [96] to assess the performance of numerical methods used to solve the incompressible Navier-Stokes equations with heat transport. The problem is defined on a squared shaped domain $\Omega=[0,L]^{2}$ with two opposite differentially heated walls and characteristic length $L=1$. Due to the small temperature gap between the walls $\delta\theta$ the Boussinesq assumption can be applied, which means that changes in the density can be neglected everywhere in the incompressible Navier-Stokes equations, apart from the buoyancy forces (gravity source term) in the momentum equation.

We consider the computational domain defined in Figure 3. Five different sets of parameters have been set according to values of the Rayleigh number $Ra=\frac{\beta\delta\theta gL^{3}}{\nu\alpha}$ between $10^{3}$ and $10^{7}$ with Prandtl number $Pr=0.71$ (see Table 2 for further details on the operating conditions). Initially, we assume a constant temperature $\theta(\mathbf{x},0)=\theta_{0}=0$ and a fluid at rest, i.e. $\mathbf{v}(\mathbf{x},0)=0$. Adiabatic boundary conditions are prescribed at the bottom and top walls, whereas the exact temperature is imposed in the two heated walls, which is $\theta_{h}=0.5$ on the left and $\theta_{c}=-0.5$ on the right, hence $\delta\theta=\theta_{h}-\theta_{c}=1$.

The staggered mesh employed has $5172$ primal elements. The simulations were run considering the numerical parameters $p=2$, $p_{\gamma}=0$, i.e. we use piecewise quadratic polynomials in space and a first order scheme in time. The steady state stopping criterion is defined as

[TABLE]

In order to analyze and compare the numerical results with available data in the literature, let us introduce the Nusselt number at the heated walls:

[TABLE]

where $\Gamma$ refers to one of the heated walls, $L$ denotes the characteristic length and $\kappa$ is the thermal conductivity with $\kappa_{0}=\kappa\left(\theta_{0}\right)$.

The Nusselt numbers obtained at both heated walls using the fully Eulerian and the Eulerian-Lagrangian schemes are shown in Table 3. The average values presented in [96], [97], [98], [99] and [6] have been included for comparison purposes. To avoid overestimation of the Nusselt number by the Eulerian-Lagrangian scheme we need to bound the time step. More precisely, the time step given by (97) is halved, thus allowing the small structures embedded in the flow to be properly tracked.

Figures 4-8 depict the numerical solution of temperature and velocity. Moreover, streamlines and vorticity contours are shown in Figure 9. A good agreement is observed between the results computed with our two different advection schemes (Eulerian upwind scheme and Eulerian-Lagrangian scheme) and with the data reported in literature (see [100, 97, 101, 102, 98, 99, 103, 6]). Comparison with available data is also carried out considering the vertical and horizontal velocities in the mid planes, which have been plotted in Figures 11-10. Furthermore, the normalized maximum velocities in the mid plane are reported in Tables 4-5.

6.3 Cavity with differentially heated cylinders

We propose a modification to the previous test case on a more complex geometry where six cylinders of radius $r=0.1$ are subtracted from the square cavity, see Figure 12.

We assume an initial fluid at rest with constant temperature $\theta_{0}=293.15K$, density $\rho_{0}=1$ and $\mathbf{g}=\left(0,-9.81\right)^{T}$. Dirichlet boundary conditions are imposed on the cylinder walls, where the exact temperature is prescribed. Specifically, the left cylinders are cooled with $\theta_{c}=292.65K$, whereas the right ones are heated with $\theta_{h}=293.65K$. The remaining walls are assumed to be adiabatic. The computational domain is paved with a primal mesh of $18396$ elements and two different settings are considered, according to the diffusion coefficients: $\left(\mu_{T_{1}},\alpha_{T_{1}}\right)=\left(4.874\cdot 10^{-4},6.86\cdot 10^{-4}\right)$ and $\left(\mu_{T_{2}},\alpha_{T_{2}}\right)=\left(4.874\cdot 10^{-6},6.865\cdot 10^{-6}\right)$.

Figure 13 shows the contour plots of temperature, horizontal and vertical velocity obtained for the first setup using the fully Eulerian and the Eulerian-Lagrangian schemes with $p=2$, $\Theta=0.51$. A good agreement is observed between the two approaches. To illustrate the mesh convergence, the results corresponding to a mesh refinement factor of 2 are included as well.

The reduced viscosity coefficient defined in the second setup leads to an unsteady flow field, developing many secondary vortices and flow instabilities. The simulation was run relying on the Eulerian-Lagrangian scheme which reduces the computational cost. Nevertheless, to capture the sub-scale structures of the flow, we have fixed the time step to be $5/2$ times greater than the one computed from (97). In Figure 14 we depict the results obtained at different output times.

6.4 Warm bubble in two space dimensions

As fourth test case we propose to solve the smooth rising bubble benchmark problem introduced in [5], which has been widely used to test numerical solvers of thermal convection problems (see [104, 105, 106, 10]). This problem assumes an initial fluid at rest and a warm bubble embedded in it. During the simulation, the bubble will rise and deform gradually acquiring a mushroom-type shape and later also developing secondary Kelvin-Helmholtz-type flow instabilities. Although no exact solution is known, the analysis of the results can be qualitatively performed taking into account the symmetry of the bubble as well as other numerical results available in the literature.

The computational domain is a square of $1$km side length. The bubble is initially placed at $\mathbf{x}_{b}=(0.5,0.35)$ km and is assigned the following initial temperature fluctuation:

[TABLE]

with $r=\left\|\mathbf{x}-\mathbf{x}_{b}\right\|$ representing the distance from the center of the bubble, $r_{b}=0.25$km being its radius and $\theta_{b}=0.5$K denoting the maximal initial amplitude of the perturbation. The temperature perturbation (107) is imposed over a background temperature which is assumed to be $\theta_{0}=303.15$K.

The use of the Eulerian scheme, without any limiter, to solve Euler equations may produce high velocity gradients in the boundary of the warm bubble. To avoid the growth of instabilities in the numerical solution, we may add a small artificial viscosity. However, as already discussed in [104, 105, 107], this leads to a smoothness of the boundary of the bubble, that is affected by numerical dissipation. Figure 15 shows the results obtained at time instants $t\in\left\{600,700,800,900\right\}$s using a constant artificial viscosity with $\mu=10^{-6}$ and a very coarse primal mesh composed of only $5172$ triangles. We run the test problem using piecewise polynomials of degree $p=4$ in order to represent the discrete solution. We observe that the evolution of the main shape of the bubble agrees with the results available in literature. Moreover, for large times, long-wave oscillations arise on top of the main structure.

An alternative strategy to get a stable solution with very little numerical dissipation is the use of a Eulerian-Lagrangian advection scheme. In this case we do not need to add any artificial viscosity, even for high order DG discretization. Moreover, the use of this methodology yields to a reduction on the computational cost of the algorithm. For the test case shown in this section, the implementation based on the Eulerian-Lagrangian advection scheme was overall 4.3 times faster than the implementation based on the classical Eulerian advection scheme. This improvement is related to the increase of the time step size allowed by the Eulerian-Lagrangian approach.

The results obtained using the Eulerian-Lagrangian approach for $\mu=0$ are depicted in Figure 16 for the intermediate time instants $t\in\left\{600,700,800,900\right\}$s. We observe that the solution at the beginning of the simulation is almost symmetric, where the small discrepancies which arise are due to the unstructured and non-symmetric mesh that is used. After $800$s mushroom-shaped structures start to grow from the cap of the main shape. Moreover, for long times small vortex structures develop on the top of the long waves resulting in secondary Kelvin-Helmholtz instabilities. At this step the flow becomes completely turbulent and therefore the initial symmetry is lost. The comparison against the results depicted in Figure 15, where a non-zero artificial viscosity was employed, show that the viscosity required by the Eulerian scheme to provide a stable solution is already too high to properly capture the small vortexes generated within the fluid flow. Therefore, the use of a Eulerian-Lagrangian scheme with zero viscosity is a key point in order to avoid the damping of Kelvin-Helmholtz instabilities. Furthermore, Figure 16, which shows the results obtained considering spatial polynomial degrees $2$ and $4$, highlights the importance of employing a high order scheme in order to capture the small secondary vortex structures.

6.5 Warm bubble in three space dimensions

We now extend the warm bubble test to a three-dimensional domain, following the setup described in [107, 106]. To this end, we define the computational domain $\Omega=\left[0,1\right]\times\left[0,1\right]\times\left[0,1.5\right]$ km. The warm bubble is initially located at $\mathbf{x}_{b}=(0.5,0.5,0.35)$ and the temperature is defined by (107) with $r=\left\|\mathbf{x}-\mathbf{x}_{b}\right\|$. Furthermore, we consider $r_{b}=0.25$km, $\theta_{0}=303.15$K and zero viscosity, i.e. $\nu=0$.

The computational mesh counts a total number of $1,331,442$ primal tetrahedral elements and we use piecewise polynomials of degree $p=3$ for the representation of the discrete solution. This amounts to a total of $26,628,840$ degrees of freedom for the pressure. The computer code was parallelized using the MPI standard and the simulation was run on 960 CPU cores of the SuperMUC-NG supercomputer at the LRZ in Garching, Germany. Figures 17-18 depict the numerical results obtained at different time instants. The three-dimensional results mimic the behavior of the warm bubble already observed in the two-dimensional setting. The differences in temperature induce a velocity field in the flow that will produce an upward movement of the bubble, which yields the development of the main mushroom-shaped structure for the temperature. However, we observe that the Kelvin-Helmholtz instabilities arising in the bi-dimensional test case do not appear in 3D. This may be due to the coarse grid resolution used here. As already noticed in Figure 16, an increase in the spatial order of accuracy of the scheme could improve the final results.

6.6 Comparison of incompressible and compressible flow solvers for a Gaussian bubble in 2D

In this section we compare the results obtained with both the incompressible and the compressible flow solver by considering a rising bubble test defined in the computational domain $\Omega=\left[-0.5,1.5\right]\times\left[-0.5,1.5\right]$ with periodic boundary conditions on the lateral boundaries. The initial temperature profile is the one of a Gaussian bubble and reads

[TABLE]

where $r=\left\|\mathbf{x}-\mathbf{x}_{b}\right\|$ is the distance from the center of the bubble $\mathbf{x}_{b}=(0.5,0.35)$; $r_{b}=0.25$ is the bubble radius and $\sigma=2$ is a halfwidth; $p_{0}=10^{5}$ is the initial pressure; and $R=287$ denotes the specific gas constant. Moreover, we consider the following operating conditions:

[TABLE]

Consequently, the main dimensionless numbers are given by

[TABLE]

The simulations were run in the Eulerian framework taking $p=3$, $\Theta=0.51$. The primal mesh consists of $3602$ control volumes. Figures 19, 20 and 21 depict the numerical results at output times $t\in\left\{10,20,30\right\}$. As expected, we notice that the temperature and velocity contours obtained with the incompressible (INS) and compressible (CNS) flow solver are really close to each other. This is confirmed by a more detailed analysis of the numerical solution computed along specific cuts along the $x-$axis in the computational domain, as we can observe in the plots included at the right hand side of Figures 19-21. Therefore, for a flow regime verifying the hypothesis related to the Boussinesq approximation, we can use both sets of equations to solve natural convection problems obtaining similar results. One of the main advantages of the incompressible model is the lower computational cost compared to the compressible one, that is approximately seven times slower. Moreover, the algorithm for the incompressible model is much simpler, since we can neglect density variations, hence avoiding the non-linearity of the system presented for the compressible case (see Section 4.2). On the other hand, the compressible model is more general. Since the algorithm has been developed to perform well at all Mach numbers, see [4], we are able to simulate a broader variety of problems including both small and large temperature fluctuations.

As second part of this test case, we compare the results obtained considering two different time discretizations. Figure 22 reports the numerical solution computed with the incompressible model, when using two sets of space-time basis functions, namely $(p,p_{\gamma})=(3,1)$ and $(p,p_{\gamma})=(3,2)$. We can observe that the theta method produces the same results of the second and third order accurate space-time DG schemes.

6.7 Locally heated cavity

As last test case, we introduce a locally heated cavity test which looks similar to the classical one reported in Section 6.2. Positive and negative heat sources on the bottom and upper wall, respectively, are imposed as boundary conditions. Due to the physical instability of the chosen boundary conditions, we drive the generation of the rising/falling structures by setting sinusoidal functions. The computational domain $\Omega=[-0.8,0.8]\times[-0.45,0.45]$, is initially discretized using $5172$ triangular elements. Then, we perform some refinements in order to show mesh convergence. We use adiabatic no-slip walls for the left and right boundaries and we prescribe

[TABLE]

on the top and bottom boundaries. Finally, we set $\mathbf{g}=\left(0,-9.81\right)^{T}$, $\mu=10^{-4}$, $\alpha=10^{-3}$, and $p=3$. The resulting time evolution of the temperature for mesh refinements by a factor of $2,3,4$ and $5$ in each space dimension is reported in Figures 23-24. We can observe the generation of some mushroom-shape structures that arise from the flow that is only driven by the temperature boundary conditions. Later in time, they start to mix and to form complex substructures. Let us remark that, for refinement levels $2$ and $3$, we obtain essentially the same results up to $t=12s$. However, at $t=14s$ we encounter some small differences that become evident at the final time, see Figure 23. However, a good agreement is observed for long time solutions when increasing the mesh resolution, see Figure 24.

7 Conclusions

In this work, we have presented a new high order accurate staggered semi-implicit discontinuous Galerkin finite element scheme for the solution of natural convection problems. The algorithm is based on the work proposed in [1, 3, 4]. A unified framework for the discretization of incompressible and compressible Navier-Stokes equations with gravity terms has been introduced. The computational cost of the global algorithm has been reduced thanks to the development of a novel Eulerian-Lagrangian approach for the treatment of the nonlinear convective terms, which leads to an unconditionally stable scheme for the incompressible Navier-Stokes system. Furthermore, the new methodology is able to handle flows with large temperature and velocity gradients. The final algorithms developed within this work have been assessed using several classical benchmarks, showing good agreement with the reference data. A detailed comparison between the fully Eulerian and the Eulerian-Lagrangian approaches for advection has been performed, highlighting the main advantages and drawbacks of both methodologies. Moreover, the numerical results obtained with the incompressible Navier-Stokes equations in combination with the Boussinesq approximation have also been validated against a numerical solution of the full compressible Navier-Stokes equations. We have shown computational results for a rising bubble in three space dimensions using more than 26 million spatial degrees of freedom, which clearly shows that the proposed computational approach can also be used on modern massively parallel distributed memory supercomputers. This test case has been set to demonstrate the potential capability of the proposed high order staggered semi-implicit DG algorithm for direct numerical simulations (DNS) on unstructured grids at moderate Reynolds numbers, i.e. the full resolution of all small scale structures present in the flow without the use of any turbulence models.

Future research will involve the development of a conservative Eulerian-Lagrangian approach for the treatment of convective and viscous terms as well as its extension to the compressible regime. Another research direction will be the use of subgrid scale turbulence models for large eddy simulations (LES) of gravity-driven flows at high Reynolds numbers.

Acknowledgments

This work was financially supported by INdAM (Istituto Nazionale di Alta Matematica, Italy) under a Post-doctoral grant of the research project Progetto premiale FOE 2014-SIES; M.T. and M.D. acknowledge partial support of the European Union’s Horizon 2020 Research and Innovation Programme under the project ExaHyPE, grant no. 671698 (call FETHPC-1-2014); W.B. was partially financed by the GNCS group of INdAM and the program Young Researchers Funding 2018. The simulations were performed on the HazelHen supercomputer at the HLRS in Stuttgart, Germany and on the SuperMUC-NG supercomputer at the LRZ in Garching, Germany. The authors acknowledge funding from the Italian Ministry of Education, University and Research (MIUR) in the frame of the Departments of Excellence Initiative 2018–2022 attributed to DICAM of the University of Trento (grant L. 232/2016) and in the frame of the PRIN 2017 project Innovative numerical methods for evolutionary partial differential equations and applications. Furthermore, M.D. has also received funding from the University of Trento via the Strategic Initiative Modeling and Simulation.

The authors would like to thank the two anonymous referees for their helpful comments and suggestions, which helped to improve the quality of this manuscript.