Scaling tests for the *BoomerAMG* solver have been performed up to 4096 cores by Baker et al. [1]. Matrices arising from the finite element and finite difference discretizations of 2D and 3D scalar diffusion problems were considered. The authors used HMIS coarsening and extended+i interpolation and showed that *l*1-scaled Jacobi, *l*1-scaled Gauss-Seidel and Chebyshev smoothers are good choices for such problems.

Weak scaling up to 125*,*000 cores has been presented in Ref. [2], where *BoomerAMG* was used as a preconditioner for a conjugate gradient solver. The test case considered is that of a 3D Laplace operator. The parameters for the AMG solver were again HMIS coarsening, extended+i interpolation and symmetric hybrid Gauss–Seidel for the smoother. Aggressive coarsening with multipass interpolation was used on the finest grid, while the problem on the coarsest level was solved by Gaussian elimination. The authors show the impact of additional parameters such as the use of 64 bits for the integers or the use of an hybrid parallel strategy with OpenMP and MPI.

In the present work, we use the *BoomerAMG* from *hypre* to precondition a GMRES solver for the pressure equation arising from the spectral element discretization of the Navier–Stokes equations. We study strong scaling up to 131,072 cores on two different supercomputers: Mira, based on the IBM Blue Gene/Q architecture, and Hazel Hen, a Cray XC40 system. The first test case considered is the flow around a NACA4412 airfoil, which we use to identify a set of best parameters for the *BoomerAMG* solver. A second test case is employed for a strong scaling study: the turbulent flow in wire-wrapped pin bundles [3, 13].

The paper is organized as follows. In Sect. 2, we introduce the discretization method and describe the preconditioning strategy and the *hypre* library. In Sect. 3, we study which set of parameters gives the fastest time to solution for the problem at hand. Using those parameters, we perform a strong scaling study for the flow in wire-wrapped pin bundles in Sect. 4. We finish with conclusions and outlook in Sect. 5

#### **2 Problem Description**

Considering an operator splitting strategy to solve the Navier–Stokes equations, the consistent pressure *p* is the solution to a Laplace problem of the form *p* = *r*. In Nek5000, this equation is discretized using the spectral element method [7] and it has been shown that it is well preconditioned by an overlapping additive Schwarz method. The preconditioner combines local problems RT kA−1 k Rk and a coarse grid problem RT 0A−1 0 R0 and is expressed as

$$\mathbf{M}^{-1} = \mathbf{R}\_0^\mathsf{T} \mathbf{A}\_0^{-1} \mathbf{R}\_0 + \sum\_{k=1}^K \mathbf{R}\_k^\mathsf{T} \mathbf{A}\_k^{-1} \mathbf{R}\_k \ .$$

where R0 and Rk are restriction operators, Ak are local stiffness matrices and *K* is the total number of spectral elements. The matrix A0 corresponds to a Laplace operator defined on the element vertices only. Because of its low number of degrees of freedom and global extent, the scalability of the coarse grid problem is mostly limited by communication and latency. We note that the term "coarse grid" here refers to the fact that *A*0 is defined on the vertices of the spectral elements only. The problem is therefore "coarse" in comparison to the solution fields, which are expanded on the Gauss–Lobatto–Legendre points inside each spectral element (typically order 10 quadrature points in each direction). When talking about the different levels arising from the coarsening phase of the AMG setup, we will use the term "coarse level" to avoid confusion.

As mentioned before, the solver of choice for large problems is currently an *inhouse AMG*, developed specifically for Nek5000 [6], whose main drawback is a setup phase by an external and serial code. The default option for the setup step of the *in-house AMG* is a Matlab code, which uses Ostrowski coarsening with norm bound, a diagonal Chebyshev smoother, applied on the second branch of the V-cycle only, and an energy-minimizing interpolation, all described in Ref. [9]. Other properties of the AMG include no smoothing on the finest level, a number of smoothing steps predefined for each level during the setup phase and a coarsest level made of one variable only. The good scalability of the *in-house AMG* is due to the fact that it automatically chooses, at run time, the fastest communication strategy at each level of the coarsening process, between three options: a pairwise exchange, a crystal router method or an allreduce operation. Previous work has shown that, when far from the strong scaling limit, the total time spent in the pressure solver is typically 85–90% of the total computational time, including the time spent in the coarse grid solver, which amounts to about 5–10% of that [11].

As an intermediate step in a previous work [12], the coarsening and interpolation steps from this Matlab code have been transferred to *BoomerAMG*, while the smoother and the solver were left unchanged. This "hybrid" serial setup was shown to significantly reduce the setup time without affecting the rapidity and scalability of the *in-house AMG* solver.

In the present work, we completely replace the *in-house AMG* by *BoomerAMG*, which allows for the whole AMG problem (setup + solver) to be performed online and in parallel. The existing code requires only limited modifications. Local contributions to the coarse grid operator are built on each process and then handed over to *BoomerAMG*, which takes care of assembling the global operator and of communication. If the operator possesses a nullspace, the solution is normalized such that the mean of the solution entries is 0. Apart from that, the critical aspect of switching to *BoomerAMG* is the choice of parameters for the setup and for the solver that match the performance of the *in-house AMG*.

#### **3 Optimal Parameter Selection**

The choice of parameters for the *BoomerAMG* solver is done by testing a set of parameters on a medium-sized test case and looking for the optimal combination. We consider the turbulent flow around a NACA 4412 airfoil at *Rec* = 400*,*000 [13] on a mesh made of 253,980 elements, with polynomial order 11, and we run the simulation for 30 timesteps. The best set of parameters is defined as the one which minimizes the time to solution for the pressure equation. This is achieved by balancing two competing aspects: the accuracy of the coarse grid solution and the total number of iterations of the GMRES solver used for the pressure equation. Since the AMG is used as a preconditioner, a high level of accuracy is not paramount. Yet, it should be sufficient to ensure efficient preconditioning. Based on results obtained with the *in-house AMG*, the initial error on the coarse grid problem should be reduced by approximately one order of magnitude. While the *in-house AMG* is designed to ensure that a given reduction in the error is attained at minimal cost, the *BoomerAMG* is designed to ensure the maximum reduction in the error occurs at a given cost. Therefore, the best choice of parameters for the *BoomerAMG* is case dependent; here we optimize this choice in the case of large 3D simulations, when the use of AMG is most relevant.

All tests are run on 4096 processors on the Blue Gene Mira at the Argonne National Laboratory. We test a total of 96 combinations of the following parameters:

We assign a letter and a number to each option which will be used to identify the method when comparing the results.

**Table 1** Five best timings for the *BoomerAMG* with corresponding parameters as compared to the *in-house AMG*

3 V-cycles for the *BoomerAMG*. Total number of iterations for the pressure solver and timings are reported for 30 timesteps

**Table 2** Five best timings for the *BoomerAMG* with corresponding parameters as compared to the *in-house AMG*

1 V-cycle for the *BoomerAMG*. Total number of iterations for the pressure solver and timings are reported for 30 timesteps

In all cases, we set a relative tolerance of 0*.*1 on the solution of the coarse grid problem and a maximum of 3 V-cycles for the AMG. Moreover, the problem on the coarsest level is solved by Gaussian elimination. Total timings and number of pressure iterations for 30 timesteps are presented in Table 1 for the *in-house AMG* and the five fastest combinations of parameters for the *BoomerAMG*. We assign a letter to each run for comparison later. The time spent in the pressure solver is reported for process 0 and the time spent in the coarse grid solver (CGS) is the maximum value over all processors for a single run. Since an allocation on Mira is always made of cores that are physically contiguous and isolated from the rest of the network, the timings suffer little noise and uncertainty. We see that the number of pressure iterations is on par with the *in-house AMG*, but the time spent in the coarse grid solver is more that three times as much. As a result, the time spent in the pressure solver is between 8 and 13% higher.

To accelerate the *BoomerAMG* solver, we set the maximum number of V-cycles to 1, instead of 3, and perform another test using the optimal parameters. This should significantly accelerate the resolution of the coarse grid solver but reduce the accuracy of the solution, therefore increasing the number of pressure iterations. The corresponding results are presented in Table 2. A surprising result comes from run d, where the number of pressure iterations has actually decreased. Since this might be the sign of an unstable solver, we discard this choice of parameters. The second best set of parameters corresponds to HMIS coarsening, extended+i interpolation, *l*1-Gauss–Seidel forward solve on the down cycle + backward solve on the up cycle for the smoother and an AMG strength threshold of 0*.*5, which we use for the rest of our simulations.

We also experimented with more aggressive non-Galerkin coarsening to change the communication pattern on the largest AMG levels and reduce communication time. However, this caused a drop of accuracy for the coarse grid solver and an increase of iterations for the pressure solver, which led to slower overall timings as a result.

We note that the time to setup the coarse grid problem, which includes building the matrix *A*0 and performing the *BoomerAMG* setup, is negligible. In the present configuration, it amounts to less than a second, whereas a single timestep takes about 20 s. Since the setup is performed once at the beginning of the simulation, we do not discuss the matter further. Furthermore, it is orders of magnitude lower than the serial versions of the setup [12].

Further analysis of the results shows that the use of classical Ruge-Stueben, Falgout, CGC or CGC-E significantly slows down the solver. Moreover, another valid choice for the smoother could have been *l*1-scaled hybrid symmetric Gauss– Seidel, whose speed is on par with our choice. Other relaxation methods, Chebyshev and *l*1-scaled Jacobi, are also consistently slower. Finally, the choice of the interpolation method does not have a significant impact and both methods give very similar results.

#### **4 Scaling Results**

As is often the case with numerical simulations, we look for the fastest path to solution for a given problem. Therefore, a strong scaling study, where the total amount of work is fixed and the number of processes is increased, is a relevant measure of the efficiency of the *BoomerAMG*. This is opposed to a weak scaling analysis, where the amount of work per process is kept constant, which is not carried out here.

We consider the turbulent flow inside a reactor assembly made of 61 wirewrapped pins, a configuration appearing in a nuclear reactor core [3]. The mesh consists of 1,650,240 elements, uses polynomial order 7 and has a complex, fully three-dimensional topology, making it a relevant test case for evaluating preconditioning strategies. The initial velocity field is turbulent and we run the simulation for 10 timesteps. Two series of tests were conducted on two supercomputers: Mira and Hazel Hen. The number of compute nodes considered is 512, 1024, 2048, 4096 and 8192 on the former machine and 256, 512, 1024, 2048, 4096 on the latter one. On both computers, the number of MPI processes per node is equal to the number of available compute cores, i.e. 16 on Mira and 24 on Hazel Hen. We use our previous defined optimal parameters for the setup of the *BoomerAMG* and we also include non-Galerkin coarsening. A drop-tolerance of 0*.*05 for sparsification is set as default on all levels, with the exception of the five finest levels, which have respective droptolerances of 0*.*0, 0*.*01, 0*,* 02, 0*.*03 and 0*.*04. This choice of parameters is motivated by the fact that, unlike with the wing case, the time for the coarse grid solver is reduced by about 25%, as tests on 6,144 cores on Hazel Hen have shown, without impacting the number of pressure iterations.

First, let us mention that the setup time for the *BoomerAMG* solver is once again negligible in comparison to the time for the entire simulation, requiring less than 3 s on any number of cores on any machine. It is also significantly lower than reading the data of the *in-house AMG*, a serial process, which takes about 80–90 s. Therefore, it does not represent a bottleneck and we do not investigate timings for the setup phase in details.

Next, we present the strong scaling results, based on a single run per core count. The reported value for the time spent in the pressure solver is the timing from core 0. The time spent in the coarse grid solver is measured on each processor and we consider the maximum value among all processes. The average time per timestep for the pressure solver on Mira is shown in Fig. 1, left plot. Unlike what was observed for the wing simulation, the choice of AMG solver does not impact the number of pressure iterations. The *in-house AMG* is slightly faster than the *BoomerAMG* on all core counts and it seems to scale marginally better. Since all other timings are the same, the reason is a faster coarse grid solver, as can be seen in Fig. 1, right plot. The *in-house AMG* achieves a better performance because it optimizes the communication process independently on each level of the AMG. On the coarsest levels, it is able to take advantage of the fast allreduce operation offered by the network of the Blue Gene architecture in hardware. On the finest levels, it picks up the fastest method between a crystal-router strategy or a pairwise exchange.

On 131,072 cores, the actual speed up for the coarse grid solver is 3*.*11 and the parallel efficiency is 0*.*194 for the *in-house AMG*, when the timings on 8,192

**Fig. 1** Rod-bundle test case on Mira, with 1,650,240 elements, leading to 12 or 13 elements per core for the largest core count. AMG parameters: HMIS coarsening, extended+i interpolation, *l*1- Gauss–Seidel forward solve on the down cycle + backward solve on the up cycle for the smoother and an AMG strength threshold of 0*.*5. Left: time spent in the pressure solver; value for process 0, averaged over the total number of timesteps. Right: time spent in the coarse grid solver; maximum value over all processes, normalized by the total number of timesteps

**Fig. 2** Rod-bundle test case on Hazel Hen, with 1,650,240 elements, leading to 16 or 17 elements per core for the largest core count. AMG parameters: HMIS coarsening, extended+i interpolation, *l*1-Gauss–Seidel forward solve on the down cycle + backward solve on the up cycle for the smoother and an AMG strength threshold of 0*.*5. Left: time spent in the pressure solver; value for process 0, averaged over the total number of timesteps. Right: time spent in the coarse grid solver; maximum value over all processes, normalized by the total number of timesteps

cores are used as references. These quantities are respectively 2*.*44 and 0*.*153 for the *BoomerAMG*. As mentioned before, the network on Mira is characterized by little noise and the uncertainty on the timings is low. Yet, these numbers are based on a single run and are only indications.

The same timings on Hazel Hen are shown in Fig. 2. Unfortunately, the node allocation on this machine can be scattered and the interconnect noise, which is shared with the rest of the computer, can be high and unpredictable. Therefore, a thorough analysis of the scaling results is not possible from a single run and the data from Fig. 2 are only indicative. Nevertheless, the runs for each AMG solver on the same number of cores are obtained using the same node allocation. This makes a comparison between the two solvers on a given core count somewhat relevant. In contrast to the results on Mira, the *BoomerAMG* is slightly faster than the *inhouse AMG* on most core counts; the exception is on 6144 cores, where the timings are almost equal. Furthermore, it is quite clear that the coarse grid solver on Hazel Hen does not scale at all. Indeed, the time spent in the coarse grid solver is almost constant from the lowest amount of cores considered.

Based on the available data, we see that beyond 24*,*576 thousand cores, the time spent in the coarse grid solver accounts for roughly half of the time spent in the pressure solver. At that point, there is about 67 elements and 35*,*000 grid points per core, which is consistent with the strong scaling limit on a similar computer as identified in Ref. [11].

#### **5 Conclusions**

We used the *BoomerAMG* solver from the *hypre* library for linear algebra to solve a global coarse problem that is part of the preconditioner for the pressure equation arising when time-integrating the Navier–Stokes equations. The set of parameters for the *BoomerAMG* setup that leads to the lowest solver time for the pressure equation is HMIS coarsening, extended+i interpolation, *l*1-Gauss–Seidel forward solve on the down cycle + backward solve on the up cycle for the smoother and an AMG strength threshold of 0*.*5. We also used non-Galerkin coarsening, with more aggressive drop-tolerance on the coarser levels, to speed up the solver. This new method replaces an existing AMG solver, which is fast and scales well but requires a setup phase done externally in serial. Strong scaling was assessed for both AMG solvers on a real large-scale test case on two supercomputers: an IBM Blue Gene/Q (Mira) and a Cray XC40 (Hazel Hen). On Mira, the *in-house AMG* leads to a faster pressure solver that the *BoomerAMG* on all core counts. The maximum difference is about 10% on 131,072 cores. This is because the *in-house AMG* is able to take advantage of the fast hardware allreduce operation on this machine at the coarsest levels of the AMG solver; in that sense the present result was expected. On Hazel Hen, however, the *BoomerAMG* is consistently faster than the *in-house AMG* and we observe the strong scaling limit to be reached at about 24,576 cores. Overall, the *BoomerAMG* is a valid alternative to the *in-house AMG*; both methods are close in terms of performance and the *BoomerAMG* has the advantage to be set up online and in parallel. In particular for modern architectures, the *BoomerAMG* is even faster than the *in-house AMG* with obvious advantages in the setup phase. All the codes developed in this work are available from the https://github.com/nicooff/nek5000/ tree/amg\_hypre\_c Github repository.

Future work will extend the use of *BoomerAMG* to mesh refinement, where an online and parallel AMG setup phase is a requirement.

**Acknowledgements** We would like to thank Aleks Obabko for his advice and for sharing the wire-wrapped pin bundle case with us. Financial support by the H2020 EU Project "ExaFLOW: Enabling Exascale Fluid Dynamics Simulation" (grant reference 671571), and the Knut and Alice Wallenberg Foundation is gratefully acknowledged. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Additional computer time was provided by ExaFLOW at HLRS Stuttgart, and by resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC Stockholm.

#### **References**

Y., Saied, F. (eds.) High-Performance Scientific Computing: Algorithms and Applications, pp. 261–279. Springer, London (2012)

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Two Decades Old Entropy Stable Method for the Euler Equations Revisited**

**Björn Sjögreen and H. C. Yee**

#### **1 Introduction, Objectives and Preliminaries**

The two decades old high order central differencing via entropy splitting and summation-by-parts (SBP) difference closure of Olsson and Oliger, Gerritsen and Olsson, and Yee et al. [2, 7, 25] is revisited. The entropy splitting is a form of skewsymmetric splitting in terms of the physical entropy of the nonlinear Euler flux derivatives. Central differencing applied to the entropy splitting form of the Euler flux derivatives together with SBP difference operators will, hereafter, be referred to as **entropy split schemes**.

The objective is to prove for the first time, in the recent definition of entropy stability based on the *L*2-energy-like norm estimate, that entropy splitting for central schemes with SBP operators are entropy stable. The proof is to replace the spatial derivatives by summation-by-parts (SBP) difference operators in the entropy split form of the equations using the physical entropy of the Euler equations. The numerical boundary closure follows directly from the SBP operator. No additional numerical boundary procedure is required. In contrast, Tadmor-type entropy conserving schemes [18] using mathematical entropies do not naturally come with a numerical boundary closure. A generalized SBP operator has to be developed [8]. Standard high order spatial central differencing as well as high order central spatial DRP (dispersion relation preserving) spatial differencing is part of the entropy stable methodology. An entropy split scheme satisfies the *L*2-energy

B. Sjögreen

© The Author(s) 2020

Multid Analyses AB, Gothenburg, Sweden e-mail<EMAIL_ADDRESS>
H. C. Yee (-) NASA Ames Research Center, Mountain View, CA, USA e-mail<EMAIL_ADDRESS>
S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_21

norm estimate readily without an added numerical dissipation term for smooth flows. For flows containing discontinuities the Yee et al. nonlinear filter approach [10–12, 14, 15, 22–25] is employed at isolated computed locations. After each full time step of the entropy split method to suppress spurious oscillations while maintaining accuracy on the remaining flow field. Since the nonlinear filter step is executed as an Euler time discretization at isolated location after the completion of a full time step of the entropy stable central scheme, entropy conservation/stability is valid almost everywhere. The efficiency and performance of the entropy stable split schemes using the physical entropies are compared with Tadmor-type entropy conservative method [18] using mathematical entropies for long time integration of a 2D smooth flows and a 3D direct numerical simulation (DNS) of turbulence with shocklets. It is found that Tadmor-type entropy conservative methods required twice the CPU time than the entropy stable split schemes using the same order of the central scheme. Comparisons among the three skew-symmetric splittings (entropy splitting [19, 20, 25], Ducros et al. splitting [1] and the Kennedy and Grubber splitting [5]) on their nonlinear stability and accuracy performance without added numerical dissipations for smooth flows is included. See [16] for additional details and comparison.

*Remarks* It is noted that the Hughes et al. formulation [4] using the Harten's idea [3] but solving the flow equations in nonconservative form in terms of the entropy variables is completely different from the entropy split schemes. The entropy split scheme solve the entropy splitting form of the Euler flux derivatives consisting of a one parameter family of conservative and a non-conservative portions in terms of the entropy variables. If the parameter satisfies the energy estimate, entropy stability is immediate. The entropy split scheme has been generalized from a perfect gas to a thermally perfect gas and gas flows consisting of linear combination of perfect gases [21, 25]. In addition, these high order schemes have been formulated in time varying deforming curvilinear grids with free-stream preservation [17, 21].

#### **2 Entropy Splitting of the Euler Flux Derivatives**

We consider the 3D equations of inviscid compressible gas dynamics

$$\mathbf{q}\_l + \mathbf{f}\_\times + \mathbf{g}\_\times + \mathbf{h}\_\times = \mathbf{0}$$

with conserved variables **q** = *(ρ ρu ρv ρw e)T* and fluxes in an arbitrary direction **k** = *(k*1 *k*2 *k*3*)* with |**k**| 2 = 1, and

$$\hat{\mathbf{f}} = k\mathbf{\hat{f}} + k\mathbf{\hat{g}} + k\mathbf{\hat{y}}\\
\mathbf{h} = \left(\rho\hat{u}\ \rho u\hat{u} + k\mathbf{\hat{u}}\ p\ \rho v\hat{u} + k\mathbf{\hat{z}}\ p\ \rho w\hat{u} + k\mathbf{\hat{z}}\ p\ \hat{u}(e+p)\right)^T,\qquad(1)$$

where *u*ˆ = *k*1*u*+*k*2*v*+*k*3*w*. The total energy is related to the pressure *p* by the ideal gas law, *e* = *p γ* −1+1 2*ρ*|**u**| 2, where *γ >* 1 is a given constant, and |**u**| 2 = *u*2+*v*2+*w*2.

An entropy is a convex function, *E(***q***)*, of the conserved variables that allows an additional conservation law,

$$E\_I + F\_\chi + G\_\chi + H\_\natural = 0,\tag{2}$$

when the solution is smooth. The entropy fluxes in the *x*-, *y*-, and *z*-directions are denoted by *F*, *G*, and *H*, respectively. The entropy variables are defined by **v** = ∇**q***E* (the notation *E***q** for the gradient will sometimes be used). The convexity of *E* ensures that these are well-defined. The Entropy conservation law (2) follows if the relation **v***T ∂***f** *∂***q** = ∇**q***F* for the *x*-direction fluxes, and similarly for the *y*- and *z*directions, holds. Moreover, the entropy variables symmetrize the equations; *∂***f***/∂***v** is a symmetric matrix.

Harten [3] considered the class of entropies

$$E = -\frac{\gamma + \alpha}{\gamma - 1} \rho (p \rho^{-\gamma})^{\frac{1}{a + \gamma}},\tag{3}$$

where *α* is a parameter. To ensure that *E* is convex, i.e., that the matrix *E***q***,***q** is positive definite, *α* is required to satisfy *α >* 0 or *α <* −*γ* . The full range for *α* was given in [25], while [3] only considered *α >* 0, and [2] used only the special case *α* = 1 − 2*γ* from *α <* −*γ* . The corresponding entropy flux in the direction **k** = *(k*1 *k*2 *k*3*)T* is

$$F = \hat{u}E.$$

The entropy variables **v** = *E***q** are straightforwardly found to be

$$\mathbf{v} = \frac{\rho}{p} \mathbf{s}^{\frac{1}{\alpha+\gamma}} \left( -\frac{\alpha}{\gamma - 1} \frac{p}{\rho} - \frac{1}{2} |\mathbf{u}|^2 \ u \ v \ w \ -1 \right)^T,\tag{4}$$

where *s* denotes *pρ*−*γ* . The conserved variables are homogeneous functions of the entropy variables (4),

$$\mathbf{q}(\theta \mathbf{v}) = \theta^{\beta} \mathbf{q}(\mathbf{v}),\tag{5}$$

where *β* = *(α* + *γ )/(*1 − *γ )*. From (5) it follows that

$$\mathbf{q}\_{\mathbf{V}}\mathbf{v} = \beta\mathbf{q} \tag{6}$$

$$
\hat{\mathbf{f}}\_{\mathbf{V}} \mathbf{v} = \beta \hat{\mathbf{f}}.\tag{7}
$$

See [3, 16] for the proof. The range of *α*, where *E***q***,***q** is positive definite, translates to *β* satisfying *β <* − *γ γ* −1 or *β >* 0.

Entropy splitting of the Euler flux derivative in the x-direction with the *y*- and *z*-directions suppressed [2, 25] is written as a weighted sum of a conservative part, **f***x*, and a non-conservative part, **fvv***x*, as

$$\mathbf{f}\_{\times} = \frac{\beta}{\beta + 1} \mathbf{f}\_{\times} + \frac{1}{\beta + 1} \mathbf{f}\_{\mathbf{v}} \mathbf{v}\_{\times}.$$

Replacing **f***x* by this split flux derivative gives

$$\mathbf{q}\_{l} + \frac{\beta}{\beta + 1} \mathbf{f}\_{\times} + \frac{1}{\beta + 1} \mathbf{f}\_{\mathbf{v}} \mathbf{v}\_{\times} = \mathbf{0}. \tag{8}$$

The entropy splitting weights the non-conservative portion of the flux derivative by 1 1+*β* . This means that the range *β >* 0 corresponds to a weight that is less than 1, whereas negative *β* leads, unphysically, to a weight that is greater than 1. The global entropy conservation can be rewritten as an *L*2-like estimate. The entropy time derivative can be rewritten as

$$\frac{d}{dt}E(\mathbf{q}) = \frac{1}{\beta + 1} \frac{d}{dt} (\mathbf{v}^T (E\_{\mathbf{q}, \mathbf{q}})^{-1} \mathbf{v})$$

by using the homogeneity (5). Due to a page limit, see [16] for further discussion. Note that it is necessary to bound the eigenvalues of *E*−1 **q***,***q** in order to make the *L*2 like norm a valid estimate.

#### **3 Semi-Discrete Entropy Split Discretization of the Euler Equations**

Consider the 1D compressible gas dynamic equations discretized on a domain *a < x<b* by a uniform grid *xj* = *(j* − 1*)Δx* + *a*, *j* = 1*,...,N*, and grid spacing *Δx* = *(b* − *a)/(N* − 1*)*. Define the semi-discrete entropy split approximation

$$\frac{d}{dt}\mathbf{q}\_j + \frac{\beta}{\beta + 1}D\mathbf{f}\_j + \frac{1}{\beta + 1}(\mathbf{f}\_\mathbf{v})\_j D\mathbf{v}\_j = 0, \quad j = 1, \ldots, N,\tag{9}$$

where *D* is a SBP difference operator. With *entropy split scheme*, we will always mean the entropy split form of Eqs. (8) discretized in space by a summation-byparts finite difference operator. The flux Jacobian matrix with respect to the entropy variables, **fv**, is symmetric. The SBP scalar product is denoted by

$$(\mathbf{u}, \mathbf{v})\_h = \Delta x \sum\_{j=1}^N w\_j \mathbf{u}\_j^T \mathbf{v}\_j,$$

where *ωj >* 0 are weights that are different from 1 only at a few points near the boundaries. The operator *D* satisfies the SBP property

$$(D\mathbf{u}, \mathbf{v})\_{\hbar} = -(\mathbf{u}, D\mathbf{v})\_{\hbar} - \mathbf{u}\_{1}^{T}\mathbf{v}\_{1} + \mathbf{u}\_{N}^{T}\mathbf{v}\_{N},\tag{10}$$

but is otherwise arbitrary. In the most common case *D* is a standard SBP centered difference operator, but other operators are possible.

A zero velocity, *u*1 = 0, *uN* = 0, boundary condition is enforced, corresponding to wall boundaries. Thanks to the SBP property of the difference approximation the derivation of entropy conservation for the continuous problem can be carried over to the discretization.

**Theorem 1** *The approximation (9) together with the boundary conditions u*1 = 0 *and uN* = 0 *conserve the global entropy in the sense that d dt* <*N j*=1 *ωjEj* = 0*.*

A method is entropy dissipative, or "entropy stable", if the computed solution satisfies (2) with inequality,

*Proof* Denote

$$\mathbf{r} = -\frac{\beta}{\beta + 1} (\mathbf{v}, D\mathbf{f})\_h - \frac{1}{\beta + 1} (\mathbf{v}, (\mathbf{f}\_\mathbf{v}) D\mathbf{v})\_h.$$

The scheme (9) can be written

$$\frac{d}{dt}\mathbf{q}\_j = P\mathbf{r}\_j,\tag{11}$$

where the projection *P* sets *u*1 = 0 and *uN* = 0. Because *P*2 = *P*, applying P to both sides of (11) gives that

$$\frac{d}{dt}P\mathbf{q} = \frac{d}{dt}\mathbf{q},$$

i.e., that *P***q** = **q** if the initial data satisfy the boundary conditions. For the entropy

$$\frac{d}{dt}E = (\mathbf{v}, \mathbf{q}\_I)\_h = (\mathbf{v}, P\mathbf{r})\_h = (\mathbf{v}, \mathbf{r})\_h - (\mathbf{v}, (I - P)\mathbf{r})\_h = $$

$$(\mathbf{v}, \mathbf{r})\_h - (P\mathbf{v}, (I - P)\mathbf{r})\_h = (\mathbf{v}, \mathbf{r})\_h,\qquad(12)$$

where we use that *P***v** = **v**. This is due to the second component of **v** is zero when the *x*-velocity, *u*, is zero, and the orthogonality *(P***v***, (I* − *P )***r***)h* = 0. The entropy equation is now of the same form as for the continuous problem, but replacing with integration-by-parts by summation-by-parts gives

$$\frac{d}{dt}E(\mathbf{q}\_j) = -F\_N + F\_1.$$

Entropy conservation follows by observing that *F* = *uE*, so that the boundary conditions imply that *F*1 = *FN* = 0.

If the boundary conditions are periodic, no SBP modification of the difference operator is needed. Entropy conservation is proved with periodic boundary conditions by direct application of the same technique as above. It can be shown that the result carries over directly to the semi-discrete approximation, since only time derivatives are used in the proof. Hence, the *L*2-like estimate

$$\frac{d}{dt} \sum\_{j=1}^{N} \omega\_j \mathbf{v}\_j (E\_{\mathbf{q}, \mathbf{q}})^{-1}\_j \mathbf{v}\_j = 0$$

is obtained for the approximation (9). It can be shown that Tadmor-type entropy conservative discretization using the Harten entropy and high order central spatial differencings are also entropy conservative methods. See Sjögreen and Yee [16] for the proof.

#### **4 Numerical Experiments**

More extensive numerical experiments are reported in the extended version of this paper [16]. Previous studies using SBP boundary closures for non-periodic boundary conditions can be found in [25]. Here selected summary results are presented.

#### **Test Case 1: 2D Compressible Euler Simulation of Smooth Flow: Isentropic Vortex Convection**

The compressible Euler equations in two space dimensions are solved with initial data

$$\rho(\mathbf{x}, \mathbf{y}) = (1 - \frac{(\mathbf{y} - 1)\boldsymbol{\beta}^2}{8\chi\pi^2}e^{1 - r^2})^{\frac{1}{\chi - 1}}\tag{13}$$

$$u(\mathbf{x}, \mathbf{y}) = u\_{\infty} - \frac{\beta(\mathbf{y} - \mathbf{y}\_0)}{2\pi} e^{(1 - r^2)/2} \tag{14}$$

$$v(x,y) = v\_{\infty} + \frac{\beta(x-x\_0)}{2\pi}e^{(1-r^2)/2} \tag{15}$$

$$p(\mathbf{x}, \mathbf{y}) = \rho(\mathbf{x}, \mathbf{y})^{\mathcal{Y}},\tag{16}$$

where *r*2 = *x*2 + *y*2, *β* = 5, *γ* = 1*.*4, *u*∞ = 1, and *v*∞ = 0. The exact solution is the initial data translated, **u***(x, t)* = **u**0*(x* − *u*∞*t,y* − *v*∞*t)*.

The computational domain is 0 ≤ *x* ≤ 18, 0 ≤ *y* ≤ 18 with periodic boundary conditions. The center of the vortex is chosen to be *(x*0*, y*0*)* = *(*9*,* 9*)*. The problem is solved in time with the classical fourth-order accurate explicit Runge–Kutta method to time *t* = 72, which corresponds to four revolutions of the vortex across the domain.

Comparisons of high order classical central split schemes with high order DRP schemes with grid refinements are reported in [13]. Due a space limitation only one grid with maximum and *L*2 error norm compared with the exact solution is shown in Fig. 1. Here C08-DS represents eighth-order central differencing applied to the Ducros et al. splitting form of the Euler flux derivatives. The corresponding eighth-order entropy splitting, entropy conservative method and Kennedy Grubber splitting are indicated by "C08-ES", "C08-EC" and "C08-KGS". If the computed solutions by "C08-DS", "C08-ES", "C08-EC" and "C08-KGS" are nonlinearly filtered by a dissipative portion of WENO7 (seventh-order weighted essentially nonoscillatory spatial method) with an adaptive flow sensor, they are indicated by C08-DS+WENO7FI, C08-ES+WENO7FI, C08-EC+WENO7FI, and C08-KGS+WENO7FI [14, 15, 22–25]. For the smooth flow without any turbulent structure, *β* = 1 for the entropy split scheme. The *β* parameter studies are reported in [9, 16, 25]. In general, for compressible shock-free turbulence and turbulence with

**Fig. 1** Inviscid 2D compressible vortex convection with 1002 grid points: comparison of maximum-norm of error vs. time for C08-DS, C08-ES, C08-EC, and C08-KGS (left, top), and C08- DS+WENO7FI, C08-ES+WENO7FI, C08-EC+WENO7FI, and C08-KGS+WENO7fFI (right top). Bottom left and bottom right are the corresponding *L*2-norm of error vs. time

shocklets, *β* lies somewhere in the range 1*.*5 *<β<* 2*.*5. In general, the optimal *β* is problem dependent. A general conclusion is that *β* should not be very large or below 1.

Other high resolution dissipative shock-capturing methods are also candidates for the nonlinear filter approach as well as other optimal WENO or ENO methods. However, with a good control of the numerical dissipation away from discontinuities, there is no need to use the more complicated and more CPU intensive shock-capturing methods. The non-split C08 without any added numerical dissipation diverges shortly after time evolution. Results by WENO5 or WENO7 are very diffusive with large maximum or *L*2 errors. For this smooth long time integration flow, entropy splitting is the most accurate method.

#### **Test Case 2: 3D Isotropic Turbulence with Eddy Shocklets**

The second numerical test problem computes decaying compressible isotropic turbulence with eddy shocklets. For high enough turbulent Mach numbers weak shocks (shocklets) develop from the turbulent motion. Here the initial turbulent Mach number is 0.6. The Navier–Stokes equations are solved using *γ* = 1*.*4. The computational domain is a cube with side length 2*π* and periodic boundary conditions in all three directions. The initial datum is a random divergence free velocity field, *ui,*0, *i* = 1*,* 2*,* 3, that satisfies

$$\frac{3}{2}\mu\_{rms,0}^2 = \frac{1}{2}\langle u\_{l,0}, u\_{l,0} \rangle = \int\_0^\infty E(k) \, dk$$

with energy spectrum

$$E(k) \sim k^4 e^{-2(k/k\_0)^2}.$$

The computations were made with *urms,*0 = 1 and *k*0 = 4. The angular brackets denote averaging over the entire computational domain. The density and pressure fields are initially constant. The Taylor-scale Reynolds number, *Reλ,*0, is 100. See [6] for definitions of the quantities and more details about the set up of the problem. The simulation is run to the final time 4.

Figure 2 shows the comparison of two splitting methods (DS and KGS), ES (entropy splitting and entropy stable) and EC (entropy conservative) using the same nonlinear filter. The time evolution of the domain averaged kinetic energy (upper left), enstrophy (upper right), temperature variance (lower left), and dilatation (lower right) are compared. All four forms of the nonlinear filter method provide similar resolution.  Over all, DS splitting is slightly less CPU intensive than ES. KGS skew-symmetric splitting is more CPU intensive than DS and ES. The EC method is around two times more expensive than DS. In addition, as the order of these methods increases, the gain in efficiency (CPU) by entropy split schemes increases.

**Fig. 2** 3D Isotropic turbulence problem with 643 grid points. Comparison of two splitting method (DS and KGS), ES (entropy splitting and entropy stable) and EC (entropy conservative) using the same nonlinear filter. Evolution of kinetic energy (upper left), enstrophy (upper right), temperature variance (lower left), and dilatation (lower right) DNS computed on 2563 grid points and filtered down to 643 resolution is considered as the reference solution

Although entropy split methods are not in conservation form but entropy conservative, Sect. 4 showed that they perform well on problems with shocklets. Over all, Extension of the entropy split scheme to other equations of state (nonperfect gas) and the MHD can be found in the original 2000 Yee et al. [25] paper. The entropies (3) can be used to construct entropy conserving schemes in conservative form. See [16] for the derivation.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **A Mimetic Spectral Element Method for Free Surface Flows**

**L. Nielsen and B. Gervang**

#### **1 Introduction**

In the last decades, CFD simulations of free surface flows have become a key tool in engineering analysis in the design of marine structures. To be able to obtain valid estimates of environmental stress on ship-wave hydrodynamics, offshore wind turbines, wave energy converters, and offshore production systems the CFD tools need to be able to account for non-linear wave-wave and wave-body interaction. Traditionally free surface water simulation has been simulated using lower order methods, however recently spectral element methods have been used [2]. In contrast to earlier work, in the present article, we simulate 2D free surface waves using a mimetic spectral element method. This ensures that the invariants of the system mass, momentum, and energy are conserved throughout the simulation.

The governing equation for incompressible, Newtonian fluids is the Navier– Stokes equation. Free surface waves can be assumed to be governed by an inviscid and irrotational fluid flow. Assuming first the fluid to be inviscid we arrive at the Euler equations,

$$
\rho \left[ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right] = -\nabla p + \rho \mathbf{g},
$$

together with the continuity equation

$$\nabla \cdot \mathbf{u} = 0.$$

L. Nielsen · B. Gervang (-)

© The Author(s) 2020

Department of Engineering, Aarhus University, Aarhus, Denmark e-mail<EMAIL_ADDRESS>
S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_22

Using the vector identity 1 2∇*(***u**·**u***)* = *(***u**· ∇*)***u**+**u**×*(*∇ ×**u***)* and using that the fluid is irrotational (∇ ×**u** = 0*)*, we can rewrite the momentum and continuity equations,

$$
\rho \frac{\partial \nabla \phi}{\partial t} = -\frac{\rho}{2} \nabla |\nabla \phi|^2 - \nabla p - \rho \mathbf{g}, \tag{1}
$$

$$
\nabla^2 \phi = 0,\tag{2}
$$

where *φ* is a vector potential defined as **u** = ∇*φ, φ* = *φ(x, z, t)*. We can now rewrite the momentum equation as,

$$\nabla \left[ \rho \frac{\partial \phi}{\partial t} + \frac{\rho}{2} |\nabla \phi|^2 + p + \rho gz \right] = 0,$$

which we can integrate in space to obtain the time dependent Bernoulli's equation.

$$
\rho \frac{\partial \phi}{\partial t} + \frac{\rho}{2} |\nabla \phi|^2 + p + \rho gz = C(t),
$$

where *C(t)* is an arbitrary function of integration. We assign *C(t)* = 0 by recalling that *φ* and *φ* + ; *C(t)dt* yield exactly the same flow. Redefining *φ* and retaining the symbol *φ* := *φ* + ; *C(t)dt* we obtain the time dependent Bernoulli's equation for the problem as,

$$
\rho \frac{\partial \phi}{\partial t} + \frac{\rho}{2} |\nabla \phi|^2 + p + \rho gz = 0. \tag{3}
$$

The governing equations for inviscid and irrotational flows for an incompressible fluid are stated through (2) and (3), where the unknowns are the velocity potential, *φ*, and the pressure, *p*. Equations (2) and (3) together with proper boundary conditions constitute a well-posed problem. The velocity potential, *φ*, can be solved from the Laplace equation and then substituted into the Bernoulli's equation to obtain the pressure field.

#### *1.1 Boundary Conditions*

The physical domain is shown in Fig. 1, where the notations are also illustrated. The fluid domain *Ω* ⊂ R*, d* = 2 is a bounded, connected domain with piecewise bathymetry *Γ b* ⊂ R*d*−1. The time domain is taken as *T* : *t* ≥ 0. The unknowns for the problem become the velocity potential and the free surface elevation *η(x, t)* : *Γ F S* × *T* −→ R. The pressure can hereafter be determined through (3).

**Fig. 1** Illustration of the physical domain with notation of the relevant quantities shown

**Fig. 2** Computational domain

The unsteady kinematic and dynamic free surface boundary conditions are given by Zakharov [8],

$$
\partial\_t \eta = -\partial\_{\underline{x}} \eta \partial\_{\underline{x}} \tilde{\phi} + \tilde{\nu} (1 + \partial\_{\underline{x}} \eta \partial\_{\underline{x}} \eta) \quad \in \quad \Gamma^{FS} \times T,\tag{4}
$$

$$
\partial\_l \tilde{\phi} = -g\eta - \frac{1}{2} ( (\partial\_\mathbf{x} \tilde{\phi})^2 - \tilde{\nu}^2 (1 + \partial\_\mathbf{x} \eta \partial\_\mathbf{x} \eta) ) \quad \in \quad \Gamma^{FS} \times T,\tag{5}
$$

where ˜signify functions defined only on the free surface. The vertical component of the velocity *ν*˜ = *∂zφ*|*z*=*η* is calculated by solving the Laplace problem (2) together with the Zakharov boundary conditions (4) and (5) on the free surface. On the bottom we have the no penetration condition,

$$
\partial\_{\overline{z}}\phi + \partial\_{\overline{x}}h\partial\_{\overline{x}}\phi = 0,\text{ for }z = -h(\mathbf{x})\text{ on }\Gamma^b.\tag{6}
$$

On the inlet and outlet boundaries *(Γ* \*Γ F S* ∪ *Γ b)* the gradient of the velocity potential is specified. The computational domain is shown in Fig. 2.

#### **2 Discretization of Governing Equations**

The developed method adopts elements from differential geometry. The unknowns of our system are described by use of *differential forms*. In a three-dimensional setting we are making use of four types of sub-manifolds: points, curves, surfaces, and volumes, both as inner and outer oriented objects, see an example in Fig. 3. The mimetic spectral element method uses an approach similar to the Galerkin method of the finite element method where the numerical residual is weighted by an arbitrary weight function. In contrast to the traditional finite element method the arbitrary weight functions are taken from the dual space of the function space used by the unknowns.

**Fig. 3** Three-dimensional dual De Rahm complex showing the four types of sub-manifolds and their different orientations

#### *2.1 Basis Functions*

For the polynomial representation we use Lagrange polynomials *li(x)* and edge polynomials *ei(ξ )*, see [5]. The Lagrange polynomials are based on a Gauss-Lobatto-Legendre (GLL) point distribution for the nodal values. The Lagrange polynomials and edge polynomials satisfy the properties,

$$l\_l(\xi\_j) = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j, \end{cases} \qquad \int\_{L\_j} e\_l(\xi) = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j, \end{cases}$$

and the edge polynomials are explicitly given in terms of the nodal Lagrange basis functions *li(x*) as

$$e\_l(\xi) = -\sum\_{k=1}^{l-1} dl\_k(\xi\_l),\tag{7}$$

where *dlk(ξ )* is the exterior derivative applied to the 0-form *lk(ξ )*. This definition of the edge polynomial also implies, see [4] and [5],

$$dl\_l = e\_l - e\_{l+1} \,. \tag{8}$$

#### *2.2 Mimetic Discretization in 2D*

If we let the 0-form *φ(*0*)* ∈ *Λ*0*(M)* be expanded as

$$\phi\_h^{(0)} = \sum\_{l,j=0}^{N} \phi\_{l,j} l\_l(\xi) l\_j(\eta),\tag{9}$$

then we can write *φ(*0*) h* as a matrix-vector product

$$\boldsymbol{\Phi}\_h^{(0)} = \left[\mathbf{L} \otimes \mathbf{L}\right] \boldsymbol{\Phi} = \mathbf{M}^{(0)} \cdot \boldsymbol{\Phi},\tag{10}$$

where **L***i,j* = *li(ξj )* and *ξj* are the Gauss-Lobatto-Legendre points, GLL points.

If we let the 1-form *u(*1*)* ∈ *Λ*1*(M)* be defined as

$$
\mu^{(1)} = \mu^{\xi} \, d\xi + \mu^{\eta} \, d\eta,\tag{11}
$$

we can expand *uξ* and *uη* using edge polynomials as,

$$\mu\_h^\xi = \sum\_{l=1}^N \sum\_{j=0}^N \mu\_{l,j}^\xi e\_l(\xi) l\_j(\eta),\tag{12}$$

$$u\_h^\eta = \sum\_{l=0}^N \sum\_{j=1}^N u\_{l,j}^\eta l\_l(\xi) e\_j(\eta). \tag{13}$$

The discrete one-form *u(*1*)* can also be written as a matrix-vector product, where **u** is evaluated in the GLL points,

$$
\mu\_h^{(1)} = \begin{bmatrix} \left[\mathbf{L} \otimes \mathbf{E}\right] & \mathbf{0} \\ \mathbf{0} & \left[\mathbf{E} \otimes \mathbf{L}\right] \end{bmatrix} \cdot \begin{bmatrix} \mathbf{u}^\xi \\ \mathbf{u}^\eta \end{bmatrix} = \mathbf{M}^{(1)} \cdot \mathbf{u},\tag{14}$$

where **E***i,j* = *ei(ξj )*.

The 2-form *P(*2*)* ∈ *Λ*2*(M)* is expanded using only edge polynomials,

$$p\_h^{(2)} = \sum\_{i,j=1}^N p\_{i,j} e\_i(\xi) e\_j(\eta) \Rightarrow \left[\mathbf{E} \otimes \mathbf{E}\right] \cdot \mathbf{p} = \mathbf{M}^{(2)} \cdot \mathbf{p}.\tag{15}$$

The Laplace equation can be reformulated using a mixed formulation, see [1], where the equilibrium equation and the constitutive relationship are separated into two equations.

$$
\nabla \phi = \mu, \quad \nabla \cdot \mu = 0. \tag{16}
$$

Writing (16) using differential geometry for a 3-D geometry we obtain,

$$d\phi^{(0)} = \mu^{(1)},\tag{17}$$

$$dq^{(2)} = 0^{(3)},\tag{18}$$

$$q^{(2)} = \star u^{(\mathrm{l})},\tag{19}$$

where we have utilized the Hodge star operator. The Hodge star operator is a map, which maps *p*-forms onto *(n* − *p)*-forms, where *n* is the dimension of the domain, *Ω*. Given a *p*-form, *λ(p)*, the hodge star maps as follows:

$$
\star\lambda^{(p)}(\mathcal{Q}^n) = \tilde{\lambda}^{(n-p)}(\mathcal{Q}^n),
\tag{20}
$$

where ˜ denotes the change of orientation of the new form. The Hodge star is also the coupling between the outer oriented domain and the inner oriented dual space, as seen in Fig. 3.

In 2-D, using differential geometry, equations (16) take the form,

$$d\phi^{(0)} = u^{(1)}, \quad \star u^{(1)} = \tilde{q}^{(1)}, \quad d\tilde{q}^{(1)} = \tilde{0}^{(2)}.\tag{21}$$

When the exterior derivative is applied to the balance equation of (21) we obtain, see [5]

$$d\tilde{q}\_h^{(1)} = \sum\_{i,j=1}^N (q\_{l,j}^\xi - q\_{l-1,j}^\xi + q\_{l,j}^\eta - q\_{l,j-1}^\eta) e\_i(\xi) e\_j(\eta),\tag{22}$$

where we have utilized (8). The equilibrium equation, the first equation in (21), is equated to a zero valued 2-form. Expanding the last equation in (21) yields,

$$\sum\_{i,j=1}^{N} f\_{i,j} e\_l(\xi) e\_j(\eta) = \sum\_{i,j=1}^{N} (q\_{lj}^{\xi} - q\_{i-1j}^{\xi} + q\_{lj}^{\eta} - q\_{lj-1}^{\eta}) e\_l(\xi) e\_j(\eta), \tag{23}$$

where *fi,j* = 0. The basis can then be cancelled and we can rewrite (23) as,

$$\mathbf{f} = \mathbb{E}^{(2,1)}\mathbf{q},\tag{24}$$

where E*(*2*,*1*)* is an incidence matrix, only consisting of 0, 1 and −1. This matrix relates the fluxes of **q** to the volume integral of the balance equation, see Fig. 4.

The first step in developing the discrete system is the formulation of the weak form, where we make use of duality pairing between an arbitrary *k*-form*, α(k)*, and

an arbitrary *(n* − *k)*-form, *β(n*−*k)*. The duality pairing is defines as,

$$\left\langle \alpha^{(k)}, \tilde{\beta}^{(n-k)} \right\rangle\_{\mathcal{Q}^n} = \int\_{\mathcal{Q}^n} a^{(k)} \wedge \tilde{\beta}^{(n-k)}.\tag{25}$$

The pairing with the *(n* − *k)*-form, *β(n*−*k)*, takes the role of a weight function in traditional finite element analysis and lives in the dual space and carry the opposite orientation. The result of duality pairing can also be represented as a matrix-vector product,

$$\boldsymbol{\beta}^{T} \cdot \tilde{\mathbf{M}}^{(n-k),T} \cdot \mathbf{W} \cdot \mathbf{J} \cdot \mathbf{M}^{(k)} \cdot \boldsymbol{\alpha} = \boldsymbol{\beta}^{T} \cdot \mathbf{M}^{(k)} \cdot \boldsymbol{\alpha},\tag{26}$$

where **W** contains the Gauss weights and **J** is the Jacobian matrix. M*(k)* is a mass matrix of the corresponding discretized *k*- and *(n*− *k)*-form pairing, and ˜*( )* denotes a matrix of opposite orientation.

Using Stokes generalized theorem [3] and applying integration by parts to the balance equation (the last equation of (21)), we obtain.

$$\int\_{\mathcal{Q}} d\tilde{q}^{(1)} \wedge a^{(0)} = \int\_{\mathcal{Q}} d\left(\tilde{q}^{(1)} \wedge a^{(0)}\right) - \int\_{\mathcal{Q}} \tilde{q}^{(1)} \wedge da^{(0)} \tag{27}$$

$$=\int\_{\partial\Omega} \left(\tilde{q}^{(1)} \wedge a^{(0)}\right) - \int\_{\Omega} \tilde{q}^{(1)} \wedge da^{(0)}.\tag{28}$$

Using duality pairing, an inner product projection for the term with the Hodge star operator, the expansions in (9)–(15), and appropriate boundary conditions we can set up the matrix system for the discrete Laplace operator as shown in (29).

$$
\begin{bmatrix}
\mathbf{0} & \mathbf{0} & \mathbb{E}^{(1,0),T}\tilde{\mathbb{M}}^{(1)} \\
\mathbf{0} & \mathbb{M}^{(1)} & \tilde{\mathbb{M}}^{(1)} \\
\mathbb{M}^{(1)}\mathbb{E}^{(1,0)}\mathbb{M}^{(1)} & \mathbf{0}
\end{bmatrix} \cdot \begin{bmatrix}
\boldsymbol{\phi} \\
\mathbf{u} \\
\tilde{\mathbf{q}}
\end{bmatrix} = \begin{bmatrix}
\mathbf{0} \\
\mathbf{0} \\
\mathbf{0}
\end{bmatrix} .\tag{29}
$$

Using the forward Euler scheme for the temporal term and pairing it with an arbitrary 0-form, *α*˜*(*0*)* , we can rewrite the Bernoulli's equation as,

$$\left\langle \left( \rho^{(2)} \wedge \frac{\phi\_n^{(0)} - \phi\_{n-1}^{(0)}}{\Delta t} + \frac{1}{2} \rho^{(2)} \wedge \left( i\_{d\phi\_{n-1}^{(0)}} d\phi\_n^{(0)} \right) + p^{(2)} = -\rho^{(2)} \wedge h^{(0)} \text{ g} \right), \left. \tilde{a}^{(0)} \right\rangle\_{\Omega} . \tag{30}$$

The density is considered a 2-form, which leaves (30) Hodge invariant. The interior product *i* is defined in [7]. The discrete version of (30) takes the form,

$$\left[\underset{\Delta t}{^{\rho}}\mathbb{M}^{(2)} + \underset{\mathbf{2}}{^{\rho}}\mathbb{M}^{(2)}\mathbb{M}\_{i}^{(1,1)}\mathbb{M}^{(2)}\right] \cdot \left[\underset{\mathbf{P}}{^{\rho}}\right] = -\rho \text{ g } \mathbb{M}^{(2)}\mathbf{h} \ + \underset{\Delta t}{\underset{\Delta t}{^{\rho}}}\mathbb{M}^{(2)}\phi\_{n-1}.\tag{31}$$

M*(*1*,*1*) i* is derived from the interior product of the two 1-forms in the convective term, and contains information of *φ* from the previous time step and consequently has to be updated at each new time step.

The simulation is initialized by first solving the Laplace equation with the prescribed boundary conditions. The initial velocity potential *φ* on the free surface, is set to *φ(x, t* ˜ = 0*)* = *x*, and the free surface height is set to *η(x, t* = 0*)* = 0. At the following time steps, the Zakharov free surface equations are solved to obtain new values of *φ*˜ as well as the free surface elevation, *η*.

#### **3 Numerical Results**

The method is first applied to a non-temporal problem without a free surface. The geometry sketched in Fig. 5 contains a cylinder in the middle of a square. On the horizontal walls of the square and the cylinder wall the no penetration condition is applied. On the left vertical boundary a fully developed velocity profile is specified and on the right vertical boundary a constant velocity potential is defined. The velocity potential *φ* and streamlines are shown in the middle section of Fig. 5. In the right part of Fig. 5 the pressure field is shown. Figure 6 shows that we obtain spectral convergence for both unknowns.

Furthermore, the balance equation ∇ · **u** = 0 (conservation of mass) is satisfied both globally and point-wise independent of polynomial order as shown in Fig. 7.

Next we apply the method to a temporal and free surface problem where we have included a bump on the bottom boundary. The Zakharov free surface equations are applied on the top horizontal boundary. In Fig. 8 the pressure field and the free surface are plotted at *t* = 1*,*100*,*200.

**Fig. 5** Left: multi-element mesh of the cylinder problem with corresponding boundary conditions. Middle: solved velocity potential *φ* in black with corresponding streamlines in red. Right: Solved pressure field from the Bernoulli equation

**Fig. 6** The two unknowns of the system, the velocity potential, *φ*, and the pressure field, *P* , are shown to carry spectral convergence

**Fig. 7** The mass balance equation of (16) (∇ · **u** = 0) is satisfied both globally and locally for any order of the expanding polynomial

**Fig. 8** Time progression of the pressure fields, *P* , at time steps *t* = 1, 100 and 200 are shown to the left. To the right the height of free surface wave *η* is shown (a scaling factor of 10 is used)

#### **4 Discussion and Conclusion**

Using an isoparametric, multi-element formulation the solution of the discretized Laplace equation shows spectral convergence. In addition, we observe that mass is conserved both globally and locally.

In (31), the discretized Bernoulli equation was kept Hodge-invariant, leaving the equation metric free. This suggests that the fundamental invariant of the equation is conserved. The Bernoulli equation conserves the total energy of the system. However, in Fig. 9 it is observed that a small amount of energy is gained and lost in a periodic manner. It is also observed that the mean energy is constant. It was possible to time integrate over very long time periods without noticing any degradation of data and we conclude that energy is conserved over long time periods even though fluctuations were observed for short time periods. In the future we plan on using a mimetic time integration scheme, which was used in [6], as well as the mimetic spatial discretization that was used in the present work.

**Fig. 9** Potential and kinetic energy is summed for the entire system at every time step and plotted against time

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Spectral/hp Methodology Study for iLES-SVV on an Ahmed Body**

**Filipe F. Buscariolo, Spencer J. Sherwin, Gustavo R. S. Assi, and Julio R. Meneghini**

This work focuses on the correlation study between a computational and physical model of an Ahmed Body with slant angle of 25◦, which generates a complex flow behaviour over the slant and back, with two vortices being generated from the side combined with separation on the slant. Physical results are from a wind tunnel test, performed by Strachan et al. [12] considering moving ground and Reynolds number of 1.7M, based on the length of the body.

CFD simulations were performed using the code Nektar++, which is an open source, spectral/hp element high-order solver, which methodology combine both mesh refinement (h), with higher polynomial order (p) for higher fidelity modelling. It employs an implicit type turbulence model using a Spectral Vanish Viscosity (iLES-SVV) model, which works as a filter for high frequencies. Same physical test conditions and tunnel test section were also considered, over a total time of 4 convective lengths, with same Reynolds number of 1.7 Million from reference experiments.

Considering the drag coefficient values for fully developed cases on the 5th and 6th polynomial order, the difference observed, compared with experimental results, was a maximum difference of 16%, however the simulation does not consider the upper support used in the experimental setup. Comparing the Spectral/hp element

F. F. Buscariolo (-)

Imperial College London, NDF-USP, London, UK e-mail<EMAIL_ADDRESS>
S. J. Sherwin Imperial College London, London, UK e-mail<EMAIL_ADDRESS>
G. R. S. Assi · J. R. Meneghini NDF-USP, São Paulo, Brazil e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
© The Author(s) 2020

S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_23

297

LES-SVV case from literature, the agreement with the experimental drag coefficient has been improved, reducing the gap from 45 to 16%. For the lift coefficient the maximum difference between the simulation results compared to experimental data is only 3%. There is also a good agreement between the LDA measurements on the end of the body with the results from the simulation. It is possible to observe a more intense vortex core on the simulation results, as compared to experimental data, which might well be explained by the upper support used to fix the Ahmed body in experimental test, which weakens the vortices.

The methodology shows promising results against the open literature once an appropriate validation study has been undertaken. Despite the relatively coarse resolution adopted the results are encouraging. Having identified an appropriate resolution, we will next consider other slant angles, to see how well these correlate with the experimental studies.

#### **1 Introduction**

Among all automotive bluff bodies in literature, the most studied one is the Ahmed Body. It was first proposed by Ahmed et al. [2], based on previous work from Morel [7], which was the first to study the behavior of slanted bluff bodies. The Ahmed body was designed to have shape similar to road vehicles and generate their main flow features, such as stagnation and separation points. The main dimensions of the Ahmed body are highlighted on Fig. 1.

Based on the results found by Ahmed et al. [2] on the variation of the slant inclination angle, Huminic and Huminic [4] states that three different flow configurations are found: from 0 to 12.5◦, the airflow over the angled surface remains fully attached before separating from the model when it reaches the vertical surface of the back end. The flow from the angled section and the side walls produces a pair of counter rotating vortices, which continue downstream; from 12.5 to 30◦, the flow over the angled section becomes highly complex. Two increased counterrotating lateral vortices are shed from the sides of the angled section with increased size, which affects the flow over the whole back end, causing a three-dimensional wake. These vortices are also responsible for maintaining attached flow over angled surface up to an angle of 30◦; from 30◦ and above, the flow is fully separated. There remains though a weak tendency of the flow to turn around the side edge of the model, a result of the relative separation positions of the flow over model top and that over the backlight side edges.

Due to some limitations on the wind tunnel and resources, Ahmed performed only force measurements on the bluff body during his experiments. In order to better understand the flow phenomena on an Ahmed Body, Lienhart and Becker [6] performed a study using Laser Doppler Anemometry (LDA), Hot-Wire Anemometry (HWA) and static pressure measurements in order to investigate the flow and turbulence structure around the Ahmed Body model for two slant angle conditions: 25 and 35◦. The main scope was to supply a detailed data set acquired under

**Fig. 1** Ahmed Body schematic drawing considering its main dimensions and 3D visualization

well-defined boundary conditions, similar to Ahmed first test, which considered a Reynolds number of 4.29 Million based on the length and static floor, to be used as reference data for numerical simulations.

Aiming to reproduce the real highway conditions of a vehicle, Strachan et al. [12], performed an Ahmed Body wind tunnel test with moving road conditions and both the aerodynamic forces and flow characteristics by time-averaged LDA were recorded. The flow conditions were slightly different from the ones used on Ahmed first test, by reducing the flow velocity to 25 m/s resulting in a Reynolds number of Re = 1.7 Million based on its length and the supports on the ground were replaced by a fixing system on the top of the tunnel, due to the rolling road simulation.

The Ahmed Body stands as one of the most used validation cases for CFD codes employed for automotive applications. Simulations employing a Reynolds Averaged Navier-Stokes (RANS) methodology are able to predict with good accuracy the drag coefficient, even for cases with complex flow topology, such as the slant angle of 25◦, with correlation factor of around 95% compared to experimental results, however the flow physics does not agree, usually under-estimating the flow features. Attempts considering more refined methodologies such as Detached Eddy Simulations (DES) and Large Eddies Simulation (LES) provide better correlation with experiments when comparing the flow structures but aerodynamic quantities values lose accuracy.

A trend that rose to improve the confidence level of CFD simulations was the high-order or high-fidelity methods, such as the spectral/hp element method [5]. The spectral/hp elemental method combines, according to Xu et al. [13], the advantages of the spectral element method, in terms of the properties of accuracy and rapid convergence, with those of the classical h-version finite element method, that allows complex geometries to be effectively captured. It also provides an attractive higherprecision approximation to solve partial differential equations.

One of the software that employs the spectral/hp element methodology is Nektar++ [9]. Nektar++ is a cross-platform spectral/hp element framework which aims to make high-order finite element methods accessible to the broader community. This is achieved by providing a structured hierarchy of C++ components, encapsulating the complexities of these methods, which can be readily applied to a range of application areas, as stated by Cantwell et al. [3]. It allows the use of high complex solution such as implicit LES (iLES) using a Spectral Vanish Viscosity (SVV) technique to stabilize the solution.

The latest achievements in the high-fidelity turbulence models around an Ahmed Body with slant angle of 25◦ are summarized in the compilation work of Serre et al. [11], in which a comparative analysis of recent simulations, conducted in the framework of a French–German collaboration on LES of Complex Flows at Reynolds number of 768,000. It compares the results obtained with different eddyresolving modelling approaches, with two LES on body-fitted curvilinear grids: LES with Smagorinsky model and wall function (LES-NWM) and Wall-resolving LES with dynamic Smagorinsky model (LES-NWR), a stabilized spectral method known as iLES-SVV, similar to the one used in this present work, which is the base of the Nektar++ code and a DES-SST approach on an unstructured grid with element number ranging from 18.5 to 40 Million. Results of the flow field shows good agreement with results measured by Lienhart and Becker [6] by a gap on the drag coefficient values of 17% for the best case and 45% for the one using iLES-SVV.

#### **2 Objectives**

The main objective of this work is to evaluate the aerodynamic behaviour in terms of the drag and lift coefficients, considering an Ahmed Body with slant angle of 25◦ using a spectral/hp elements method methodology as shown on Fig. 2. To achieve this, we first present a mesh study, evaluating two different size refinements referred as h-refinement for each of those, we employ three high-order surface mesh values to improve curvature representation. As the spectral/hp element method has also the possibility to improve the solution by increasing the polynomial order and consequently the number of degrees of freedom, we also evaluated three high polynomial orders for each mesh case, in a total of eighteen load cases.

**Fig. 2** Representation of the Ahmed Body with slant angle of 25◦

All load cases employ moving ground condition and Reynolds number of 1.7 Million, based on the length of the body. Due to the same conditions considered, results are compared with experiments performed in the study of Strachan et al. [12].

#### **3 Spectral/hp iLES-SVV**

In this work, Nektar++ is used to run an implicit LES simulation using spectral/hp method. In this method, the domain is first divided into non-overlapping elements, offering geometric flexibility and allows for local refinement. Simulations were performed using the incompressible Navier–Stokes solver employing a velocity correction scheme, combined with a Continuous Galerkin (CG) projection. More details are presented by Cantwell et al. [3].

The mathematics behind Nektar++ basically considers the numerical solution of partial differential equations (PDEs) of on a domain , which may be geometrically complex, for some solution u. Practically, takes the form of a d-dimensional finite element mesh consisting of elements K*i*, embedded in a space of dimension dc, such that d ≤ dc ≤ 3, with = *uiKi* is an empty set or an interface between elements of dimension dbar *<* d. The PDE problem is solved then in the weak sense, considering that u|K*i* must be smooth with at least a 1st-order derivative. Therefore is required that u|K*i* is in the Sobolev space W1*,*2(K*i*) equivalent to H1 (K*i*), according to Adams [1]. For a continuous discretisation, we impose C0 continuity along element interfaces.

We assume the solution can be represented as u*δ*(x) = < *n u*ˆ*nn*(x), a weighted sum of N trial functions *n*(x) defined on and the problem becomes that of finding the coefficients *u*ˆ*n*. The approximation u*δ* does not directly give unique choices for the coefficients *u*ˆ*n*. To achieve this, a restriction is placed on the residual so that its L2 inner product, with respect to the test functions *n*(x), is zero. For a Galerkin projection it is chosen that the test functions are the same as the trial functions, that is *n* = *n*. As outlined previously, to construct the global basis *n* it is first considered the contributions from each element in the domain. Each K*i* is mapped from a standard reference space K is between [−1, 1] by a parametric mapping *χ*e: K becomes K*i* given by x = *χe(ξ )*, where K is one of the supported region shapes, and *ξ* are d-dimensional coordinates representing positions in a reference element, distinguishing them from x which are d-dimensional coordinates in the Cartesian coordinate space.

The next step is to construct a local polynomial basis on each reference element with which to represent solutions. For 3D regions, a tensorial basis may be used, where the polynomial space is constructed as the tensor-product of one-dimensional bases on segments, quadrilaterals or hexahedral regions.

Spectral/hp element discretisation generally lead to approximations that have low dissipation and low dispersion per degree of freedom when compared to lowerorder methods. As stated by Xu et al. [13], in solving advection-diffusion equations and nonlinear partial differential equations such as advection-dominated flows, at marginal resolutions, oscillations appear that may render the computation unstable. Artificial viscosity has been used in may discretisation methods to suppress wiggles associated with high wavenumbers has been broadly and effectively used in simulations using the Fourier method. A related concept is the so-called SVV, which was originally proposed based on a second-order diffusion operator for spectral Fourier methods. SVV has been explicitly regarded as a turbulent model of implementing iLES under the assumption that the action of subgrid scales on the resolved scales is equivalent to strictly dissipative action stated by Sagaut [10], even though SVV is not explicitly designed as a subgrid-scale model. An example of a 1-D SVV kernel is:

$$\text{Df} = \begin{cases} 0, \text{ p} \le \text{Pcut} \\ \exp\left(-\frac{(p-P^2)}{(p-Pcut)^2}\right), \text{ p} > \text{Pcut} \end{cases} \tag{1}$$

where P is the total number of modes employed and P*cut* is the cutoff polynomial order. SVV with the kernel function Df can be regarded as a low-pass filter. We see that the SVV dissipation added to the high mode numbers with respect to the spectral element discretisation does indeed yield dissipation at the global high wave number scales of the solution.

For this work, we employed a novel CG-SVV scheme with DGKernel, proposed by Moura et al. [8] where he dissipation curves of CG of order *p* are match to those of DG with order *p* − 2, eliminating non-smooth dissipation characteristics arising from CG dissipation when considering high Reynolds number.

#### **4 Simulation Methodology**

We first define the coordinate system as X the streamwise direction, Y the vertical direction and Z the spamwise direction. The Ahmed Body length of 1.044 m is defined as 1 AL. The virtual wind tunnel dimensions are 2*.*74 × 1*.*66 m for the test section and total length of the domain of 4 AL, similar to Strachan et al. [12] study. The Ahmed Body model back in placed on X = 0, inlet position at X = −2 AL and outlet position at X = 2 AL. A schematic setup is shown on Fig. 3.

In terms of boundary conditions, velocity was normalized to 1 in order to match the Reynolds number previously stated and set as the inlet boundary condition. The outlet was set as pressure high-order outlet condition and the floor was also set with the same velocity of the free stream in order to reproduce the moving floor effect. The top and outer side wall and the Ahmed Body wall are set as no slip condition and a symmetry condition. Total simulated time is 7 convective lengths AL, which means that the flow is able to cross the whole domain.

This study evaluates two mesh configuration considering different h-refinements and referred as Original and Refined meshes and for each of those, three highorder surface mesh settings: 4th, 5th, and 6th order, generating six different meshes. All mesh files were generated by NekMesh, which is Nektar++ high-order mesh generator. In both Original and Refined meshes, cases two refinement zones were generated, where the first one, defined as the Ahmed Body refinement, ranges from 0.3 AL before the beginning of the geometry and 0.3 AL after the end of the body, in a total length of 1.6 AL. The second refinement, defined as the Wake Refinement region, intercepting the first refinement in 0.3 AL before the end of the body, to 1.3 AL after the end of the body, in order to fully capture the flow phenomena in the separation region, with same total length of 1.6 AL, as illustrated on Fig. 4.

The Original Mesh has total number of elements for half model around 95,000. For the Refined mesh, the boundary layer setup was the same and the dimensions

**Fig. 3** Schematic representation of the boundary condition on the Ahmed Body simulation

**Fig. 4** Plane Z = 0 representation of mesh refinement regions. Ahmed Body Refinement region highlighted in yellow and Wake Refinement region highlighted in black

**Fig. 5** Plane Z =0 mesh refinement comparison between two h-refinement cases. (**a**) Original mesh. (**b**) Refined mesh

were kept the same in terms of sizing. Refined mesh setup, giving a total of 310,00 elements. Details of both meshes are shown in Fig. 5.

Most of the commercial CFD code employ low order methods and the highest order polynomial interpolation for the solutions usually seem is 3rd. The mesh plays the major role for complex simulations such as LES, leading to elevate number of elements to reach a reliable result. To make use of the flexibility of the spectral/hp element methods, we proposed solutions considering polynomials with order higher than 3rd within the previous mesh refinement studies as the higher order polynomials increase the degrees of freedom and resolution of the mesh. For the Nektar++ implicit LES simulations using the Incompressible Navier–Stokes solver evaluated three different polynomial expansions, 4th, 5th and 6th orders, referred here as P4, P5 and P6. In summary, 18 load cases were evaluated using HPC with 432 CPUs for each case.

#### **5 Results**

#### *5.1 Drag Coefficient Comparison Results*

The drag coefficient for the 18 cases evaluated, considering 9 from the Original mesh with 95,000 elements considering fourth, fifth and sixth polynomial order and the Refined mesh case with 310,000 elements also considering fourth, fifth and sixth polynomial order expansion, with maximum RMS and compared with experimental results are shown on Fig. 6.

From Fig. 6 it is possible to observe that for the drag coefficient, P4 polynomial expansion considering both mesh cases presented mean drag results around 35% higher than the experimental results. For the P5 cases, considering again both Original and Refined mesh cases, the error was reduced to 5% however results change the trend from over-predicted to under-predicted when the mesh is refined further. The cases considering P6 polynomial expansion presented the same trend for both mesh cases, highlighting its consistency although the mean error when compared to experiments increases to 16%.

#### *5.2 Lift Coefficient Comparison Results*

Similar to the drag coefficient graph, in Fig. 7 the lift coefficient for the all evaluated cases is shown, considering Original and Refined meshes and fourth, fifth and sixth polynomial order expansion. Maximum RMS is also plotted for all cases and compared with experimental results from Strachan et al. [12].

**Fig. 6** Drag coefficient for the 18 evaluated test cases. On the left, average values for Original mesh, considering fourth, fifth and sixth polynomial expansions (P4, P5 and P6). On the right average, values for Refined mesh, considering fourth, fifth and sixth polynomial expansions (P4, P5 and P6)

**Fig. 7** Lift coefficient for the 18 evaluated test cases. On the left, average values for Original mesh, considering fourth, fifth and sixth polynomial expansions (P4, P5 and P6). On the right average, values for Refined mesh, considering fourth, fifth and sixth polynomial expansions (P4, P5 and P6)

Analyzing Fig. 7, we observe that the h-refinement from Original mesh to Refined mesh lead to results closer to experimental values when adopting P4 as the polynomial expansion basis. For both P5 and P6 polynomial expansions, lift coefficient results present good agreement with experimental data, with maximum mean error of 5%.

#### *5.3 Flow Structure Comparison*

In terms of the polynomial order expansions for the solution, combined with the 6th order surface mesh, the results present focus on the Refined mesh case, once they improved the correlation for the P4 polynomial expansion within experimental results and kept similar trend for P5 and consistent results for P6 in terms of drag and lift coefficient prediction. An initial comparison in terms of flow structures is present by the Q-Criterion of 350 coloured by pressure, comparing the Refined mesh case, considering P4, P5 and P6 polynomial expansions in Fig. 8.

**Fig. 8** Iso-Surface of Q-Criterion = 350 colored by pressure on the Ahmed Body with slant angle of 25◦, considering Refined mesh and 6th order surface mesh for fourth, fifth and sixth polynomial expansions (P4, P5 and P6)

**Fig. 9** Contour of Lambda 2 on the plane x/L = 0 on the back of the Ahmed Body with slant angle of 25◦, considering Refined mesh and 6th order surface mesh for fourth, fifth and sixth polynomial expansions (P4, P5 and P6)

From Fig. 8 it is possible to visualize that P4 is unable to define the vortex on the side of the slant, explaining also the difference in terms of both drag and lift coefficients, compared to experimental results. Results for P5 show the side vortex clearly defined and P6 is also able to capture the lower vortex, detailed on the lower image, which is not present in the studies considering the Ahmed Body, but they are important to understand the behaviour with the moving floor. Figure 9 shows a contour of Lambda 2 to illustrate the lower vortex detail on the plane x/L = 0, on the back of the Ahmed Body.

We next focus only on the Refined mesh with 6th order surface mesh for P5 and P6, once they were able to predict both lower and top vortices. Due to nature of the wind tunnel with moving ground used by Strachan et al. [12], the model had to be fixed on the top by a steam, which can be removed in the drag coefficient calculations, however it might change the flow topology over the slant, as stated by the authors themselves.

Comparing the plane x/L = 0.076 with the measurements of the flow velocity on x direction U normalized by the free stream velocity of Lienhart and Becker [6] with static floor without the steam on the upper portion with results of Strachan et al. [12], it is possible to notice intensity changes in the U normalized velocity and this is attributed by the last due to the upper support.

As the simulations do not included the upper support, but do include the moving ground, the expected results are the top portion to be similar to Lienhart and Becker [6] measurements and the lower part, correlated with measurements of Strachan et al. [12], which both P5 and P6, proved to have good agreement in terms of normalized U velocity, shown in Fig. 10.

Similar comparison is presented on Fig. 11 for the vortex intensity on the slant, on the plane x/L = 0, on the back of the Ahmed Body for vertical velocity V normalized by the free stream velocity. In this, the simulations close correlate to Lienhart and Becker [6] study, due to the absence of the steam support but for this case, the higher polynomial order expansion P6 is able to capture more scales than the P5 for the core of the main vortex, highlighting the gain of resolution of the high-order simulations.

**Fig. 10** Contour of U velocity normalized by free stream velocity on the plane x/L = 0.076 of the Ahmed Body with slant angle of 25◦, comparing LDA measurements of Lienhart and Becker [6] with static floor without the steam (left), results of Strachan et al. [12] with moving floor and steam support (middle left), Refined mesh and 6th order surface P5 (middle right) and Refined mesh and 6th order surface P6 (right)

**Fig. 11** Contour of V velocity normalized by free stream velocity on the plane x/L = 0 of the Ahmed Body with slant angle of 25◦, comparing LDA measurements of Lienhart and Becker [6] with static floor without the steam (left), results of Strachan et al. [12] with moving floor and steam support (middle left), Refined mesh and 6th order surface P5 (middle right) and Refined mesh and 6th order surface P6 (right)

#### **6 Conclusions**

Within the advances in CFD codes, confidence level and computational power, aerodynamic simulations are applied in almost every automotive company. The reason is very simple: reduced development cost and time, which is an enormous advantage in a competitive market.

High-fidelity simulations are becoming a reality for complex industrial cases in order to improve resolution and results in a reliable response time, such as presented for the Ahmed Body on this work.

On the meshing definition study, the surface mesh order seems not to influence the results in terms of aerodynamic quantities, presenting similar trend for same polynomial order, as the Ahmed Body geometry has curved surfaces only on the front portion.

Still on the mesh definition, as the h-refinement increases from Original to Refined mesh, the drag coefficient values for P4 and P6 remains unchanged and P5 values switched from positive to negative. We conclude that consistency is shown for P4 and P6 cases but P6 presented the most reliable results, with a maximum deviation of 16%. For the lift coefficient, results for P4 improved as the h-refinement increased and kept similar values for both P5 and P6, where the best agreement was found for the case considering Refined mesh with 6th order surface mesh and P6 as the polynomial expansion.

Flow structure results focus only the Refined 6th order surface mesh, where the main expected features were captured by P5 and P6 cases. It was confirmed by those two simulation cases that the lower portion has similar behaviour of the moving ground test conducted by Strachan et al. [12], however the top portion close correlate to Lienhart and Becker [6] experiments, as the simulation cases allow the body to be fixed without the upper support used in the experiment. This fact might also explain the difference from the simulation results with the literature experiments, as the simulation allows idealized configurations.

For all simulation cases, half of the body is being simulated and a symmetry plane is set on the middle portion. From Fig. 12, which shows the normalized U velocity on a line of coordinate y/L=0.15 on the plane x/L = 0, we observe that

**Fig. 12** Normalized U velocity distribution over a line at coordinate y/L = 0.15 on the plane x/L = 0 at the back of the Ahmed Body with slant angle of 25◦, comparing LDA measurements of P5 (red), P6 (orange), Lienhart and Becker [6] (dark green) and Strachan et al. [12] (light green)

simulation has good agreement with experimental results, with a small distortion as it gets closer to the symmetry plane.

**Acknowledgements** We should acknowledge the HPC facilities at Imperial and also under the UK Turbulence Consortium.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **A High-Order Discontinuous Galerkin Solver for Multiphase Flows**

**Juan Manzanero, Carlos Redondo, Gonzalo Rubio, Esteban Ferrer, Eusebio Valero, Susana Gómez-Álvarez, and Ángel Rivero-Jiménez**

#### **1 Introduction**

Multiphase flow is not a canonical problem, therefore different models can be found in the literature. Volume Of Fluid (VOF) model [9] is amongst the simplest. It defines a single set of momentum equations shared by all phases, whilst the volume fraction (fraction of a particular infinitesimal control volume which is occupied by each phase) is tracked throughout the domain following an advection equation. Phase-field methods [11] conserve the simplicity of VOF whilst increasing the physical meaning of the evolution equation of the fluids present in the simulation. The volume fraction is substituted by a phase-field parameter, which identifies each phase. In this work, the Cahn–Hilliard equation [4] is chosen to model the evolution of the phase-field parameter.

The introduced model is discretised in space using a high-order discontinuous Galerkin method. These methods have been gaining popularity for the discretisation of conservation laws, such as the Navier–Stokes equations [5–7, 13, 16, 22, 26]. Specifically, we use a Discontinuous Galerkin Spectral Element Method (DGSEM) [2] that allows the generation of provably stable schemes [8]. These schemes provide enhanced robustness when compared to classical high-order methods [17–20]. As far as the temporal discretisation is concerned, we use an efficient implicit-

J. Manzanero (-) · C. Redondo · G. Rubio · E. Ferrer · E. Valero

ETSIAE-UPM – School of Aeronautics, Universidad Politécnica de Madrid, Madrid, Spain

Center for Computational Simulation, Universidad Politécnica de Madrid, Boadilla del Monte, Madrid, Spain

e-mail<EMAIL_ADDRESS>
S. Gómez-Álvarez · Á. Rivero-Jiménez Repsol Technology Lab, Móstoles, Madrid, Spain

S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_24

explicit approach that permits maintaining the time step restriction of a typical one phase Navier–Stokes solver. It should be noticed that similar approaches to model multiphase flows have been proposed in the past, see for example [29], where an algorithm to model N immiscible incompressible fluids with high-order methods is described. However, according to the authors knowledge, this is the first implementation using the DGSEM.

The rest of the paper is organised as follows: in Sect. 2 the governing equations of the model are described. In Sect. 3 the numerical techniques to discretise the described model are introduced. Finally, in Sect. 4 the results of two validation test cases are shown.

#### **2 Governing Equations**

In this work we model multiphase flows with a phase field approach. The flow field is modelled by means of the incompressible Navier–Stokes equations. The evolution of each of the fluids is modelled with the Cahn–Hilliard equation, which defines a phase field variable, *φ* ∈ [−1*,* 1], that identifies spatial coordinates occupied by fluid 1, *φ* = −1, fluid 2, *φ* = 1, or an interface *φ* ∈ *(*−1*,* 1*)*. The value of the thermodynamic properties of the fluids at each spatial coordinate can be computed as:

$$\rho(\phi) = \rho\_1 \left(\frac{1-\phi}{2}\right) + \rho\_2 \left(\frac{1+\phi}{2}\right), \quad \eta(\phi) = \eta\_1 \left(\frac{1-\phi}{2}\right) + \eta\_2 \left(\frac{1+\phi}{2}\right), \tag{1}$$

where *ρi* is the density of fluid *i* whilst *ηi* is the dynamic viscosity of fluid *i*. The complete system is built considering first the momentum equation,

$$\frac{\partial \left(\rho \mathbf{v}\right)}{\partial t} + \nabla \cdot \left(\rho \mathbf{v} \mathbf{v}\right) = -\nabla p + \frac{1}{Re} \nabla \cdot \left(\eta \left(\nabla \mathbf{v} + \nabla \mathbf{v}^{T}\right)\right) + \frac{3}{\sqrt{2}\kappa ReCa} \mu \nabla \phi + \frac{1}{Fr^{2}} \rho \mathbf{e}\_{\mathbf{g}},\tag{2}$$

with velocity **v**, static pressure *p*, Reynolds number *Re* = *ρ*1*u*0*L η*1 (where *uo* is a reference velocity whilst *L* is a reference length), Capillary number *Ca* = *η*1*u*0 *σ* (where *σ* represents the surface tension), Froude number *F r* = √*u*0 *gL*, (where *g* is the gravity acceleration) and **e***g* is the gravity direction. Second, an artificial compressibility method [25] is used to couple the divergence-free condition,

$$\frac{\partial p}{\partial t} + \frac{\rho\_0}{\rho\_1} \frac{1}{M\_0^2} \nabla \cdot \mathbf{v} = 0,\tag{3}$$

where *ρ*0 = max *ρ*1*, ρ*2 is a reference density, and *M*0 is the artificial compressibility Mach number. Third, the Cahn–Hilliard equation for the phase field,

$$\frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \mathbf{v}) = M \nabla^2 \mu, \quad \mu = -\phi + \phi^3 - \varepsilon^2 \nabla^2 \phi,\tag{4}$$

with *M* the mobility, and *ε* the interface width, the two free parameters of the model. In (2) and (4), *μ* represents the chemical potential. Moreover, this equation is designed to minimize the free-energy functional [4], F,

$$\mathcal{F}(\phi, \nabla \phi) = \int\_{\Omega} \left( \frac{1}{4} \left( 1 - \phi \right)^2 \left( 1 + \phi \right)^2 + \frac{1}{2} \varepsilon^2 |\nabla \phi|^2 \right) d\mathbf{x}. \tag{5}$$

Note that the set of Eqs. (2)–(4) is written in non-dimensional form, where the thermodynamic variables of fluid 1 are taken as reference values, e.g.,

$$\rho(\phi) = \left(\frac{1-\phi}{2}\right) + \frac{\rho\_2}{\rho\_1} \left(\frac{1+\phi}{2}\right), \quad \eta(\phi) = \left(\frac{1-\phi}{2}\right) + \frac{\eta\_2}{\eta\_1} \left(\frac{1+\phi}{2}\right). \tag{6}$$

The set (2)–(4) can be written as an advection-diffusion system:

$$\frac{\partial \mu}{\partial t} + \nabla \cdot \mathbf{F}(u) = \nabla \cdot \mathbf{F}\_v(u, \mathbf{g}) + S(u, \mathbf{g}), \tag{7}$$

where *u* = *(φ, ρ***v***, p)* is the state vector, **g** = *(***g***φ,* **g***v,* **g***μ)* = *(*∇*φ,* ∇**v***,* ∇*μ)* is the gradients vector, **F***(u)* and **F***v(u,* **g***)* are the inviscid and viscous fluxes respectively, and *S(u,* **g***)* is a source term,

$$\mathbf{F}(\boldsymbol{\mu}) = \begin{bmatrix} \phi \mathbf{v} \\ \rho \mathbf{v} \mathbf{v} + p \boldsymbol{I}\_3 \\ \frac{\rho\_0}{\rho\_1} \frac{1}{M\_0^2} \mathbf{v} \\\\ S(\boldsymbol{\mu}, \mathbf{g}) = \begin{bmatrix} \mathbf{g}\_{\boldsymbol{\mu}} \\ \eta \left( \mathbf{g}\_v + \mathbf{g}\_v^T \right) \\ 0 \end{bmatrix}, \end{bmatrix},$$

$$S(\boldsymbol{\mu}, \mathbf{g}) = \begin{bmatrix} 0 \\ \frac{3}{\sqrt{2} \varepsilon \boldsymbol{Re} C \boldsymbol{a}} \mu \mathbf{g}\_{\boldsymbol{\phi}} + \frac{1}{\boldsymbol{Fr}^2} \rho \mathbf{e}\_{\boldsymbol{g}} \\ 0 \end{bmatrix}. \tag{8}$$

#### **3 Numerical Methods**

The numerical implementation of (2)–(4) is performed using a high-order discontinuous Galerkin scheme for the spatial discretisation (DGSEM variant) and an implicit-explicit Euler scheme for the time discretisation.

#### *3.1 Spatial Discretisation Using a Nodal Discontinuous Galerkin Scheme (DGSEM)*

Discontinuous Galerkin (DG) schemes (see [15]) are constructed by tessellating the domain in non-overlapping elements, where the solution is approximated using polynomials of an arbitrary order, *N*. In this particular implementation, we use a nodal variant of the DG method, and we restrict ourselves to hexahedral elements.

In each element we approximate the solution using polynomials written in a set of local spatial coordinates *ξ* = *(ξ, η, ζ )* ∈ [−1*,* 1] 3, which are related to the physical space by a transfinite mapping,

$$\mathbf{x} = (\mathbf{x}, \mathbf{y}, z) = \mathbf{X}(\boldsymbol{\xi}) = \mathbf{X}\left(\boldsymbol{\xi}, \boldsymbol{\eta}, \boldsymbol{\xi}\right). \tag{9}$$

Using the local coordinates, we write the solution using tensor product Lagrange polynomials,

$$\left.u(\mathbf{x})\right|\_E \approx U(\boldsymbol{\xi}) = \sum\_{i,j,k=0}^{N} U^{ijk}(t) l\_i(\boldsymbol{\xi}) l\_j(\boldsymbol{\eta}) l\_k(\boldsymbol{\zeta}),\tag{10}$$

where the time-dependent coefficients *Uijk (t)* are the nodal values of the solution *U*, and *lj (ξ )* are the Lagrange polynomials based on a set of Gauss points {*ξj* } *N j*=0. To handle curvilinear geometries, we use a mapping **X** that transforms local and physical spaces. With this mapping, we can construct covariant **a***i* and contravariant **a***i* basis, and their associated Jacobian *J* , and metrics matrix M:

$$\mathbf{a}\_{l} = \frac{\partial \mathbf{X}(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}\_{l}}, \quad \mathbf{a}^{l} = \nabla \boldsymbol{\xi}\_{l} = \frac{1}{J} \mathbf{a}\_{j} \times \mathbf{a}\_{k}, \quad J = \mathbf{a}\_{l} \cdot \left(\mathbf{a}\_{j} \wedge \mathbf{a}\_{k}\right), \quad \mathcal{M} = [J \mathbf{a}^{\xi}, J \mathbf{a}^{\eta}, J \mathbf{a}^{\zeta}]. \tag{11}$$

Following [14], we transform the system of Eqs. (7) to local coordinates,

$$
\frac{\partial}{\partial t} \left\{ \begin{aligned} J\phi \\ J\rho \mathbf{v} \\ Jp \end{aligned} \right\} + \nabla\_{\xi} \cdot \left\{ \begin{aligned} \mathcal{M}^{T} \mathbf{v} \phi \\ \mathcal{M}^{T} \rho \mathbf{v} + \mathcal{M}^{T} p \mathbf{1}\_{3} \\ \mathcal{M}^{T} \frac{1}{M\_{0}^{2}} \mathbf{v} \end{aligned} \right\} = \nabla\_{\xi} \cdot \left\{ \begin{aligned} &M\mathcal{M}^{T} \mathbf{g}\_{\mu} \\ \frac{1}{Re} \mathcal{M}^{T} \left( \eta \left( \mathbf{g}\_{\mathbf{v}} + \mathbf{g}\_{\mathbf{v}}^{T} \right) \right) \\ &0 \end{aligned} \right\}, \\ + J \begin{Bmatrix} &0 \\ \frac{1}{Fr^{2}} \rho \mathbf{e}\_{\mathbf{g}} + \frac{3}{\sqrt{2 ReCa\varepsilon}} \mu \mathbf{g}\_{\phi} \\ 0 \end{aligned} \right\}, \tag{12}
$$

with gradients,

$$J\mathbf{g}\_{\mathbf{v}} = \mathcal{M}\nabla\_{\xi}\mathbf{v},\ \ J\mathbf{g}\_{\phi} = \mathcal{M}\nabla\_{\xi}\phi,\ \ J\mathbf{g}\_{\mu} = \mathcal{M}\nabla\_{\xi}\mu,\tag{13}$$

and the chemical potential definition,

$$J\mu = -J\phi + J\phi^3 - \varepsilon^2 \nabla\_\xi \cdot \left(\mathcal{M}^T \mathbf{g}\_\phi\right). \tag{14}$$

We obtain the DG scheme replacing the continuous solution by their polynomial counterpart (10), then multiplying (12), written in compact form (7), by a polynomial test function (with same order *N* as the solution) *ϑ*, and we integrate the result in one element *E* = [−1*,* 1] 3,

$$\int\_{E} J \vartheta \frac{\partial U}{\partial t} + \int\_{E} \vartheta \nabla\_{\xi} \cdot \mathbf{F}(U) = \int\_{E} \vartheta \nabla\_{\xi} \cdot \mathbf{F}\_{v}(U, \mathbf{G}) + \int\_{E} J \vartheta S(U, \mathbf{G}). \tag{15}$$

Next, we integrate by parts the terms containing divergences, which yields surface integrals. Since the solution is discontinuous at the inter-element faces, we replace the surface flux by a *numerical flux*, **F***0*,

$$
\int\_{E} J \vartheta \frac{\partial U}{\partial t} + \int\_{\partial E} \vartheta \mathbf{F}^{\star} \cdot \hat{n} dS - \int\_{E} \nabla\_{\xi} \vartheta \cdot \mathbf{F} = \int\_{\partial E} \vartheta \mathbf{F}\_{v}^{\star} \cdot \hat{n} dS
$$

$$
$$

where *∂E* represents the six surfaces of the element E. For the inviscid numerical flux **F***0*, we use the exact Riemann solver derived in [1], whilst for the viscous numerical flux we use the Symmetric Interior Penalty (SIP) method [27], with the penalty parameter value derived in [24] and recently discussed for the DGSEM in [21]. In (16), *n*ˆ is the surface outward normal vector in local coordinates. To obtain the evolution equations for each nodal degree of freedom *Uijk* , we let *ϑ* = *li(ξ )lj (η)lk(ζ )*, and compute the integrals using the Gauss quadrature points (and weights{*wi*}) associated to the interpolation points (which provide an accuracy of 2*N* + 1),

$$\begin{split} \left.j^{j,l}\frac{dU^{ij,k}}{dt} + \frac{F^{\star}\_{\star}}{w\_{l}}(\xi,\eta\_{j},\xi\_{k})I\_{l}(\xi)\right|\_{\xi=-1}^{\xi=1} + \frac{F^{\star}\_{\star}}{w\_{j}}(\xi\_{l},\eta\_{l},\xi\_{k})I\_{j}(\eta)\Big|\_{\eta=-1}^{\eta=1} \\ \left. + \frac{F^{\star}\_{\star}}{w\_{k}}(\xi\_{l},\eta\_{j},\xi)I\_{k}(\xi)\right|\_{\xi=-1}^{\xi=1} \end{split}$$

$$\begin{split} -\sum\_{m=0}^{N} \left(\frac{w\_{m}}{w\_{l}}D\_{ml}F^{mk}\_{x} + \frac{w\_{m}}{w\_{j}}D\_{mj}F^{imk}\_{y} + \frac{w\_{m}}{w\_{k}}D\_{mk}F^{ijm}\_{z}\right) &= \\ \frac{F^{\star}\_{\star}/k}{w\_{l}}(\xi\_{lN} - \delta\_{l0}) + \frac{F^{\star}\_{\star}/k}{w\_{j}}(\delta\_{jN} - \delta\_{j0}) + \frac{F^{\star}\_{\star}/k}{w\_{k}}(\delta\_{kN} - \delta\_{k0}) \\ - \sum\_{m=0}^{N} \left(\frac{w\_{m}}{w\_{l}}D\_{ml}F^{mk}\_{v,x} + \frac{w\_{m}}{w\_{j}}D\_{mj}F^{imk}\_{v,v} + \frac{w\_{m}}{w\_{k}}D\_{mk}F^{ijm}\_{v,z}\right) &+ J^{ijk}S^{ijk}, \end{split} \tag{17}$$

where *Fijk* = *F (Uijk )* and *Fijk v* = *Fv(Uijk ,* **G***ijk )*, being **G***ijk* the nodal values of the gradient **G**. The symbol *δik* represents the Kronecker delta. The derivation matrix *Dij* is defined as *Dij* = *l j (ξi)*. To compute the gradient **G**, we perform the weak formulation of (13),

$$\int\_{E} J \, \mathbf{r} \cdot \mathbf{G} = \int\_{\partial E} U^{\star} \mathcal{M}^{T} \cdot \mathbf{r} dS - \int\_{E} U \nabla\_{\xi} \cdot \left( \mathcal{M}^{T} \cdot \mathbf{r} \right), \tag{18}$$

where *τ* is an arbitrary vector test function (from the order *N* polynomials space). Since we use the SIP method, we use solution averages to couple inter-element fluxes, *U0* = {{*U*}}. All the integrals involved in (18) are computed discretely similar to those in (16), i.e.,

$$\begin{split} \left. J^{ljk} \mathfrak{r}\_{ljk}^{d} G\_{d}^{ljk} = \frac{U^{\star}(\boldsymbol{\xi}, \boldsymbol{\eta}\_{j}, \boldsymbol{\xi}\_{k}) J a\_{d}^{\vec{k}}(\boldsymbol{\xi}, \boldsymbol{\eta}\_{j}, \boldsymbol{\xi}\_{k})}{w\_{l}} I\_{l}(\boldsymbol{\xi}) \right|\_{\boldsymbol{\xi}=-1}^{\boldsymbol{\xi}=1} \\ + \frac{U^{\star}(\boldsymbol{\xi}\_{i}, \boldsymbol{\eta}, \boldsymbol{\xi}\_{k}) J a\_{d}^{\vec{\eta}}(\boldsymbol{\xi}\_{i}, \boldsymbol{\eta}, \boldsymbol{\xi}\_{k})}{w\_{j}} I\_{j}(\boldsymbol{\eta}) \Big|\_{\boldsymbol{\eta}=-1}^{\boldsymbol{\eta}=1} + \frac{U^{\star}(\boldsymbol{\xi}\_{i}, \boldsymbol{\eta}, \boldsymbol{\eta}, \boldsymbol{\xi}) J a\_{d}^{\vec{\xi}}(\boldsymbol{\xi}\_{i}, \boldsymbol{\eta}\_{j}, \boldsymbol{\xi})}{w\_{k}} I\_{k}(\boldsymbol{\xi}) \Big|\_{\boldsymbol{\xi}=-1}^{\boldsymbol{\xi}=1} \\ - \sum\_{m=0}^{N} \left( \frac{w\_{m}}{w\_{l}} J a\_{d}^{\vec{k}, ljk} D\_{mi} U^{mjk} + \frac{w\_{m}}{w\_{j}} J a\_{d}^{\eta, ljk} D\_{mj} U^{m\vec{k}} + \frac{w\_{m}}{w\_{k}} J a\_{d}^{\xi, ljk} D\_{mk} U^{l\boldsymbol{j}m} \right). \end{split} \tag{19}$$

The gradient nodal values *Gijk d* are introduced in the viscous fluxes *Fv(Uijk ,* **G***ijk )* of (17) hence completing the discretisation of (16). Note that one needs to compute **g***φ* before computing *μ* and its gradient **g***μ*.

#### *3.2 Time Integration Using IMplicit–EXplicit (IMEX) and Runge–Kutta Schemes*

The time integration of (17) is performed with a combination of forward and backwards Euler and explicit Runge–Kutta schemes. On the one hand, the Navier– Stokes equations are integrated by means of a third order explicit Runge–Kutta (RK3) scheme [28]. On the other hand, the Cahn–Hilliard equation is integrated with a combination of explicit RK3 for the phase field advection, forward Euler for the chemical free-energy, and backwards Euler for the interfacial energy,

$$\frac{\phi^{n+1} - \phi^n}{\Delta t} + \nabla \cdot \left(\mathbf{v}\phi\right)^{RK3} = \nabla^2 \left(-\phi^n + \left(\phi^n\right)^3 - \varepsilon^2 \nabla^2 \phi^{n+1}\right). \tag{20}$$

The reason behind this choice, is that the numerical stiffness of the bi-Laplacian (∇4*φ*) operator prevents from using an explicit method, as restricts the time-step *t*

to unpractical values. We only treat implicitly the interfacial energy since it yields a constant Jacobian matrix, represented by *J* ∇2 . In particular, the linear system to solve is,

$$\left[J^{\nabla\_2} + \frac{I}{\Delta t}\right] \phi^{n+1} = \frac{\phi^n}{\Delta t} - \nabla \cdot \left(\mathbf{v}\phi\right)^{RK3} + \nabla^2 \left(-\phi^n + \left(\phi^n\right)^3\right). \tag{21}$$

The Jacobian matrix is computed numerically (see [3]) and a LU factorisation is performed only at the first time step. In each following iteration, the RHS of (21) is computed and the linear system is solved by means of forward and backward substitutions. Both the LU factorisation and the forward and backward substitutions are performed with the library MKL-PARDISO [23].

#### **4 Validation**

The proposed methodology is tested with two test cases. First, the validity of the discontinuous Galerkin discretisation of the Cahn–Hilliard equation is tested with a benchmark spinodal decomposition problem [12]. Second, the validity of the coupled Cahn–Hilliard/Navier–Stokes system is tested with a two dimensional rising bubble test [10].

#### *4.1 Spinodal Decomposition*

This test problem considers an initial mixture of two fluids. These fluids are immiscible, therefore they tend to separate to minimise their free energy (5). As stated before, the geometry, initial condition and fluid parameters are taken from [12]. In particular, the initial condition for this benchmark problem is:

$$\phi(\mathbf{x}, \mathbf{y}) = -0.05 \left[ \cos \left( 0.105 \mathbf{x} \right) \cos \left( 0.11 \mathbf{y} \right) + \left[ \cos \left( 0.13 \mathbf{x} \right) \cos \left( 0.087 \mathbf{y} \right) \right]^2 \right. \tag{22}$$

$$+ \cos \left( 0.025 \mathbf{x} - 0.15 \mathbf{y} \right) \cos \left( 0.07 \mathbf{x} - 0.02 \mathbf{y} \right) \Big]. \tag{22}$$

The physical domain is a "T" shape with a total height of 120 units, a total width of 100 units, and horizontal and vertical section widths of 20 units (Fig. 1). No-flux boundary conditions are applied at the boundaries. Following [12] mobility is set to *M* = 10, whilst the interface width is set to *ε* = 3*.*16. The physical domain is discretised with an unstructured mesh of 326 elements and a polynomial order of *N* = 4. For the time discretization, we use a time step *t* = 10−3.

Figure 1 shows qualitatively how the different phases separate, whilst Fig. 2 shows quantitatively the evolution of the total free energy with time. In Fig. 2

**Fig. 1** "T" domain for the spinodal decomposition. Initial condition (left figure) and evolution with time (the right figure is the steady-state solution)

**Fig. 2** Evolution of total free energy (5) with time

the results of this work are compared with those obtained in [12], validating the proposed method.

#### *4.2 Rising Bubble*

This test case considers a bubble of light fluid submerged in a heavy fluid, both subjected to a gravitational field. Following [10] the initial configuration, see Fig. 3, consists of a bubble of radius *r* = 0*.*25 centred at [0*.*5*,* 0*.*5] in a [1 × 2] domain. A no-slip boundary condition is used at the top and the bottom of the domain whilst a free slip condition is enforced at the vertical walls. Following [10], the Reynolds number is set to *Re* = 35 whilst *σ* and *ε* are set to 24.5 and 0.03125 respectively (this gives a Eötvös number *Eo* = 10) whilst both density and viscosity ratios are set to *ρ*1*/ρ*2 = *μ*1*/μ*2 = 10. The gravitational acceleration is *g* = 0*.*98. The problem is discretised with 16 × 32 elements with a polynomial order of *N* = 4, and a time step *t* = 4 · 10−6.

**Fig. 3** Initial condition of the rising bubble test problem

**Fig. 4** Evolution of the center of mass of the bubble with time

This test case is quantitively compared with the results of [10] in Fig. 4 with satisfactory results. It should be mentioned that the benchmark results of [10] are obtained with a sharp-interface model which may explain the small disagreement in the evolution of the center of mass shown in Fig. 4.

#### **5 Conclusions**

A method to model incompressible two phases flows is introduced. The model solves the incompressible Navier–Stokes equations coupled with the Cahn–Hilliard equation to track the evolution of the different fluids. The model is discretised in space using a discontinuous Galerkin spectral element method (DGSEM) whilst an efficient implicit-explicit approach is used to advance in time. The validity of the model is shown with two test cases. A spinodal decomposition benchmark problem is solved to validate the Cahn–Hilliard solver whilst a rising-bubble test problem is solved to validate the coupled Cahn–Hilliard–Navier–Stokes system. Both test cases are solved showing good agreement with the literature, and proving the accuracy and robustness of the proposed method.

**Acknowledgements** This work has been partially supported by REPSOL under the research grant P180021090. This work has been partially supported by Ministerio de Economía y Competitividad under the research grant TRA2015-67679-C2-2-R and under the research grant EUIN2017-88294 (Gonzalo Rubio). The authors acknowledge the computer resources and technical assistance provided by the Centro de Supercomputación y Visualización de Madrid (CeSViMa).

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **High-Order Propagation of Jet Noise on a Tetrahedral Mesh Using Large Eddy Simulation Sources**

**M. A. Moratilla-Vega, V. Saini, H. Xia, and G. J. Page**

#### **1 Introduction**

Due to the rapid expansion of the commercial aviation industry, authorities have been tightening the legislation for aircraft noise. For instance, the European Commission has set a 65% reduction goal of overall aircraft noise from the year 2000 to 2050 [1]. The noise generated by the jet exhaust is one of the main contributors to the overall aircraft noise, especially during take-off [2]. Moreover, in new generation ultra-high by-pass ratio turbofan engines the increased interaction between the engine jet and the high-lift devices can potentially affect the noise field [3]. Thus, our overall aim is to develop and investigate an accurate and efficient method for the prediction of far-field jet noise in installed jet configurations.

Rapid growth in computing power during the last decades has enabled the use of scale resolving numerical simulations for jet noise research at a reduced cost than most experimental campaigns. Conventionally, 2nd-order numerical schemes combined with surface integral techniques, particularly the Ffowcs Williams-Hawkings (FW-H) method [4] have been widely adopted for predicting the far-field noise, due to its simplicity and low cost. However, defining the envelope surface used in the FW-H method is not always trivial in complex configurations [5], for example, installed jets on aircraft wings. Also, the results may be overly sensitive to the size, shape and location of these surfaces. Now, directly resolving the Navier-Stokes (NS) equations for sufficiently accurate far-field jet noise results is prohibitively expensive [6]. LES using finite volume 2nd-order accurate schemes has proven to be reliable and robust for solving jets' near field, but large numerical dispersion and dissipation error makes them less suitable for the propagation of the

M. A. Moratilla-Vega · V. Saini (-) · H. Xia · G. J. Page

Rolls-Royce University Technology Centre, Loughborough University, Loughborough, UK e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
© The Author(s) 2020

S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_25

sound waves to the far field. High-order methods provide more accurate propagation due to their reduced numerical error but are insufficiently robust for simulating complex jet flows. Therefore, we have used a coupled approach in which a finite volume LES solver is used to obtain the acoustic sources, which are then transferred to a high-order acoustic solver that propagates noise to the far-field.

The spectral/*hp* DG method [7] is capable of providing high-order accuracy and handling mixed mesh elements types such as tetrahedra and hexahedra, thus providing a potential solution to geometrically complex acoustic problems. The solver based on this approach is *AcousticSolver* of the Nektar++ framework [8, 9]. The LES code HYDRA and acoustic code *AcousticSolver* have been coupled and validated using hexahedral elements [10, 11]. A similar coupling strategy has been used previously for jet noise [12] and combustion noise on tetrahedral grids [9].

In this paper, our focus is on two aspects: (1) estimates of mesh design for the high-order solver using a canonical two-dimensional (2D) case and (2) comparison of three-dimensional (3D) turbulent isolated jet-noise results on a tetrahedral grid and a comparable hexahedral grid using the coupling approach. From the perspective of our near future work, the tetrahedral grid results provide motivation and parameters for the set-up of the coupled methodology for jet-flap interactions.

#### **2 Numerical Methods and Solvers**

In this section, the details of the high-order spectral/*hp* DG solver employed to solve the APE equations are provided followed by a brief description of the LES code that solves the filtered compressible NS equations. Finally, the coupling of the two is briefly mentioned.

#### *2.1 APE Solver*

*Equations for Propagation* The acoustic perturbation equations (APE) solved here are the ones proposed by Ewert and Schröder [6] in the APE-4 form. These equations describe the transport of acoustic fluctuations in a linearized form, where the source terms can be non-linear, and can be written as:

$$
\partial\_t p' + \overline{c}^2 \nabla \cdot \left(\overline{\rho} \mathbf{u}' + \overline{\mathbf{u}} \frac{p'}{\overline{c}^2}\right) = \overline{c}^2 q\_c,\tag{1}
$$

$$
\partial\_l \mathbf{u}' + \nabla \left( \overline{\mathbf{u}} \cdot \mathbf{u}' \right) + \nabla \left( \frac{p'}{\overline{\rho}} \right) = \mathbf{q}\_m,\tag{2}
$$

where *p* , **u** are the acoustic pressure and acoustic velocity vector respectively and *c* is the speed of sound. The time-averaged quantities are denoted by the over-bar and acoustic fluctuations are primed. The left-hand side of (1) and (2) represents the advection of waves in the mean flow. The right-hand side describes different sources that may be present in a generic aeroacoustic problem.

Finally, the source terms, *qc* and **q***m* are defined as:

$$q\_c = -\nabla \cdot \left(\rho' \mathbf{u}'\right)' + \frac{\overline{\rho}}{c\_p} \frac{Ds'}{Dt},\tag{3}$$

$$\mathbf{q}\_m = -\left(\boldsymbol{\omega} \times \mathbf{u}\right)' + T'\nabla\overline{s} - s'\nabla\overline{T} - \left(\frac{\nabla\left(\mathbf{u}'\right)^2}{2}\right)' + \left(\frac{\nabla \cdot \underline{\mathbf{r}}}{\rho}\right)'.\tag{4}$$

These terms are classified into four categories:

In this paper, only the Lamb vector **L** is considered as a source term because it is the dominant contributor for isothermal applications with strong vortical motions (shear layers and wakes), as demonstrated in [12, 13].

*Numerical Solver* The solver used for the above APE equations is called *AcousticSolver*, which is part of the open-source Nektar++ framework [8]. The solver employs a high-order, spectral/*hp* element method with a DG formulation [7]. In short (for details see [9]), the present DG method works as follows:

The scheme used here to solve the Riemann problem is a local Lax-Friederichs scheme as defined in [9]. The temporal discretisation is performed using a 4th-order Runge-Kutta scheme. A numerical sponge layer [14] is set up using source terms to dampen out the outgoing acoustic waves smoothly, thus minimising reflections from the boundaries of the domain.

#### *2.2 LES Solver*

The LES is performed using the in-house code of Rolls-Royce plc., HYDRA [15] that solves Favre-filtered unsteady compressible Navier-Stokes equations [10]. It is a density-based, spatially 2nd-order accurate finite volume cell-vertex code used for propulsion and turbomachinery applications. More details on the set-up of the spatial scheme used can be found in [10]. For the temporal discretisation, a 2ndorder, four-stage Runge-Kutta explicit algorithm is employed. The size of the time step is chosen to keep the Courant number less than unity. The code is capable of solving arbitrary mesh topologies which is beneficial for complex geometries. The sub-grid scale model is chosen as *σ*-model [16] with model constant *Cσ* = 1*.*35 [17].

#### *2.3 Coupling of Solvers*

The 3D data from LES mesh is transferred and interpolated onto the APE mesh in real time. The interpolation is necessary because two solvers have different meshes designed specifically to capture flow and acoustics. The transfer-interpolation process takes place in parallel. This is achieved using an MPI based coupling strategy with the open-source library CWIPI [18]. More details on the coupling mechanism are provided in Lackhove et al. [9] and Moratilla-Vega et al. [10]. Note that larger time steps can be used for *AcousticSolver* since it is not restricted to resolve the small flow structures.

#### **3 Test Cases**

Two cases are presented here. First, a canonical noise propagation case due to a well-defined vortex-pair source run on *AcousticSolver* alone. A study of numerical error by changing the mesh and polynomial expansion order (P) is performed. The second case uses the mesh parameters from the first to propagate noise generated by an isolated jet in an LES simulation. This case provides validation of the coupling for a 3D turbulent jet noise case on a fully tetrahedral mesh. The results are then compared to the ones obtained with the FW-H technique.

#### *3.1 Spinning Vortex Pair*

The case is an acoustic wave propagation problem in two dimensions where the source is mathematically well-defined, as in the original work on APE [6]. The case is run with standalone *AcousticSolver*. The source is in the form of twopoint vortices at a distance of *r*0 from the origin, rotating with a circulation . An analytical solution of the induced acoustic field was found by Müller and Obermeier [19] as:

$$
\tilde{p}' = \frac{\rho\_{\infty} \Gamma^4}{64\pi^3 r\_0^4 c\_{\infty}^2} H\_2^{(2)}(kr), \tag{5}
$$

where, *H(*2*)* 2 is the Hankel function of 2nd-order and second kind, the rotation period is defined as *T* = 8*π*2*r*2 0 */* ; the angular velocity as *ω* = */*4*π r*2 0 and the Mach number as *Mr* = */*4*π r*0*c*∞. The real part of Eq. (5) gives the pressure fluctuations. Ewert and Schröder [6] found the source-term based on the Lamb vector that represents the acoustic field for this case as:

$$\mathbf{q}\_m = -\frac{\Gamma^2 \mathbf{e}\_r(t)}{8\pi^2 \sigma^2 r\_0} \sum\_{l=1}^2 (-1)^l \exp\left(-\frac{|\mathbf{r} + (-1)^l \mathbf{r}\_0(\mathfrak{H})|^2}{2\sigma^2}\right), \quad \sigma \approx r\_0,\tag{6}$$

where, *r* = *(x, y)T* , *r*0 = *r*0*er*, *er* = *(*cos *θ ,*sin *θ )T* and *θ* = *ωt*.

The computational domain considered is circular and extends to 250*r*0. The source parameters are set as in [6] i.e. */(c*∞*r*0*)* = 1*.*6 and *Mr* = 0*.*1273. Simulations are run until the pressure fluctuations reach *r* = 200*r*0 in order to minimise the boundary effects. All the elemental meshes consist of triangles and a modified hierarchical Jacobi polynomial basis [8].

First, a reference simulation with polynomial order 4 (P4) is run on a fine uniform mesh. The resulting acoustic pressure field is plotted along a diagonal line in Fig. 1a. The result from this simulation matches its analytical counterpart well. Further comparisons are made with respect to this well-resolved P4 numerical

**Fig. 1** (**a**) Pressure field along *x* = *y* line, (**b**) solution points, and contours of the source term in the source region and the pressure field in the propagation region

**Table 1** Details of the test cases run

*Nr* is the elements in the radial direction

**Fig. 2** Acoustic pressure and relative error comparison of the test cases. (**a**) Pressure along *x*=*y* line. (**b**) Relative error (as moving average)

result, henceforth called as "P4 reference". For the test cases, elemental meshes are coarsened radially and the polynomial order is elevated in the propagation region, such that, the solution points-per-wavelength (referred as "ppw") distribution is similar in the radial direction. The radial growth rate is kept ∼1.023 with geometric distribution in all the cases. In the source region, the mesh is kept the same with P1 expansion for all cases. This allows having a smooth transition of solution points distribution when crossing from one region to the other. A sample P2 mesh and contours are shown in Fig. 1b. Table 1 summarises the test runs.

Figure 2a compares the pressure fields in different test cases with the P4 reference. As expected, the P1 simulation shows a considerable reduction in the amplitude. P2 and P4 preserve this quantity more accurately. For simplification, we unify the dissipation and dispersion error by calculating the overall relative error as a moving average (M.A.) over bins of ∼30*r*0. This is plotted in Fig. 2b. For P1 and P2 simulations, a 2% error limit is reached around 50*r*0 (ppw ∼ 9) and 85*r*0 respectively (ppw∼ 6.4). P4 simulations remain below this limit under the present conditions (note ppw ∼ 5 at 150*r*0). The values of ppw for different P agree with those suggested in [20] and provide an estimate for mesh design in different polynomial order setting.

Note that for the given ppw, mesh expansion rate and Riemann solver, we did not observe reflections of the acoustic signals on the inter-element boundaries. A caveat of the present study is that we calculate total numerical error here for brevity, however dissipation and dispersion error could be studied separately as done in [20] on a one-dimensional advection study.

#### *3.2 3D Turbulent Isolated Jet Noise*

As a step forward towards noise prediction of installed jets, an isolated jet is simulated using a tetrahedral mesh for *AcousticSolver* to verify the capabilities of the present methodology for complex 3D cases. Note that the coupling is already validated on a cylinder in cross flow and a cylinder-airfoil interaction case in [11].

*Jet Flow* The LES performed is described here briefly since the same is detailed in [10]. An isothermal turbulent jet issuing from a circular cross-section nozzle at Mach 0.9 and Reynolds number *Re* = 10*,*000 (based on jet bulk velocity *Uj* and jet diameter *Dj* ) is considered. Following Shur et al. [21], the present LES domain is cylindrical in shape and extends as *x/Dj* = [−5*,* 100] and *r/Dj* = [0*,* 50]. The mesh has 190 × 75 × 49 nodes in the axial, radial and azimuthal directions respectively. It is refined in the shear layer development area and coarsened towards the outer boundaries. Figure 3a shows a central cross-section of the LES mesh. The inlet boundary condition is a total pressure profile.

*Jet Acoustics* The acoustics domain is cubical to facilitate control on mesh growth. It extends as [−5*,* 40]*Dj* in streamwise direction and [−25*,* 25]*Dj* in transverse directions. Noise propagation on two different grids is compared: fully hexahedral ("hexa") and fully tetrahedral ("tetra"). The former mesh consists of 107 × 69 × 69 elements in the streamwise and transverse directions respectively [10]. The tetra grid is generated to give a similar distribution as the hexa mesh in the vicinity of the

**Fig. 3** Cross-section view through the centre of the jet nozzle. (**a**) LES mesh elements. (**b**) *AcousticSolver* meshes

jet, providing 300,000 elements in total. Figure 3b shows the two meshes where it is seen that the nominal element size in the tetra mesh is slightly larger away from the jet nozzle. Results on the hexa grid (P4) are available from [10] and calculations are performed on the tetra grid in this study. The expansion type utilised is a P4 modified Jacobi basis [8]. A numerical sponge layer [14] of thickness 3*Dj* is applied at the outer boundaries to avoid reflections of the outgoing waves. A factor of 3 in time step size is used as compared to the compressible LES. In line with Sect. 3.1 and [20], a value of ppw ∼ 5 is chosen for accurately resolving frequencies up to a Strouhal number *St* = 0*.*9.

It is already demonstrated in [10] that the LES flow quantities are in acceptable agreement with the high-order LES study of Shur et al. [21]. The noise propagation is calculated using the FW-H method [4] in addition to the present coupled approach. The nominal cut-off *St* for the integral surface defined is ∼0.3 based on the 22 ppw criterion [6]. Figure 4 shows a visual comparison between the acoustic pressure field computed by LES alone and coupled LES-APE (on two meshes). Figure 4a, b qualitatively show that the coupled LES-APE has retained more acoustic content (especially at higher frequencies) due to lower numerical error. This difference is more pronounced in the direction perpendicular to the jet centre-line. Qualitatively comparable results are obtained on the tetra mesh as depicted in Fig. 4c.