Aerodynamics simulations of three-dimensional inviscid flow using curvilinear discontinuous Galerkin method on unstructured meshes

Abstract

Over the last decades, the discontinuous Galerkin (DG) method has demonstrated its excellence in accurate, higher-order numerical simulations for a wide range of applications in aerodynamics simulations. However, the development of practical, computationally accurate flow solvers for industrial applications is still in the focus of active research, and applicable boundary conditions and fluxes are also very important parts. Based on curvilinear DG method, we have developed a flow solver that can be used for solving the three-dimensional subsonic, transonic and hypersonic inviscid flows on unstructured meshes. The development covers the geometrical transformation from the real curved element to the rectilinear reference element with the hierarchical basis functions and their gradient operation in reference coordinates up to full third order. The implementation of solid wall boundary conditions is derived by the contravariant velocities, and an enhanced algorithms of Harten-Lax-van Leer with contact (HLLC) flux based on curved element is suggested. These new techniques do not require a complex geometric boundary information and are easy to implement. The simulation of subsonic, transonic and hypersonic flows shows that the linear treatment can limit the accuracy at high order and demonstrates how the boundary treatment involving curved element overcomes this restriction. In addition, such a flow solver is stable on a reasonably coarse meshes and finer ones, and has good robustness for three-dimensional flows with various geometries and velocities. For engineering practice, a reasonable accuracy can be obtained at reasonably coarse unstructured meshes.

1 Introduction

Computational Fluid Dynamics (CFD) has become a very important and useful design tool in aerodynamics [1, 2]. The theoretical basis of CFD includes fluid mechanics, numerical analysis, computer technology and other fields. Based on the various calculation method [3–8], the approximate solution of the fluid governing equation is obtained by using the powerful computing capability of computer. In the past decades, due to the increasing demand for complex fluid flow simulations, great effort has been done by the CFD community in order to increase the accuracy and reduce the calculation costs [9–15]. It is now well understood that numerical errors play a key role in the accuracy of CFD results [10, 12]. In recent years, there has been an increasing interest in the application of high-order discretization methods with curved meshes, which can reduce the numerical errors to aerodynamic flows [16]. Curved meshes are commonly employed to provide a satisfactory representation of the domain boundary with only a moderate number of mesh elements, so that the polynomial degree can be increased while keeping the global number of degrees of freedom (DOFs). The role that curved meshes play in obtaining accurate solutions when combined with high-order numerical discretization methods has been demonstrated, such as [16–20]. Therefore, numerical discretization methods with curved meshes supporting arbitrary approximation orders have received an increasing amount of attention.

Initially developed for advection problems, the discontinuous Galerkin (DG) method has recently been applied in the field of CFD with great success [21–24]. The DG method combine two advantageous features commonly associated with finite element and finite volume methods [25]. As in classical finite element methods, accuracy is obtained by means of high-order polynomial approximation within an element rather than by wide stencils as in the case of finite volume methods. Conservation is, similarly to finite volume methods, enforced through the use of a flux at element boundaries, where a Riemann problem is solved [17]. Furthermore, the discretization lends itself to local mesh adaptation and efficient parallelization on modern distributed memory computer architectures. The implementation of efficient DG method on curved meshes is an open field of research and crucial for developing accurate numerical schemes. In the framework of the steady-state compressible inviscid flow, the necessity of a higher-order treatment of curved wall boundaries was put in evidence by Bassy and Rebay, and is now generally accepted [26, 27]. The curved meshes method via an isoparametric parametric approximation to curve boundary is rather convenient to use for DG schemes [22, 26]. A low-storage curvilinear DG method was proposed and analyzed in [28], where the geometric factors were included in both solution and test function spaces with a provable convergence under a mild condition on the mesh. Recent works by Chan et al. [29, 30] rely on reference frame polynomial spaces introducing weight-adjusted L²-inner products in order to recover high-order accuracy. In [20], a simple boundary treatment, which can be implemented as a modified DG scheme defined on triangles adjacent to boundary, was proposed. Even though integration along the curve is still necessary, integrals over any curved element are avoided. Blended isogeometric DG method formulated on elements that exactly preserve the CAD geometry have also been recently proposed in [31].

However, three-dimensional flows with curved meshes are rarely studied and the way to impose solid wall boundary conditions in curved meshes is not mentioned. Proper boundary treatment is critical in CFD. Its importance is obvious because it is the boundary conditions that determine the flow characteristics once the governing equations are given, at least for steady flow problems [32]. Even though the boundary conditions via a correct representation of the normal to the geometry is rather convenient to use for DG schemes [17, 27], the acquisition of normal is not very easy and needs to provide geometric information, which is not provided by the general mesh generator, especially when the boundary geometry is three-dimensional. Therefore, it is difficult to extend this curved meshes method to engineering applications. A highly adaptable numerical flux format is needed in aerodynamic simulations to accommodate a wide variety of flow velocities, including subsonic flow, transonic flow, supersonic flow and hypersonic flow, and so on. However, the implementation of numerical fluxes generally depends on Lax-Friedrichs scheme and Roe flux in the curved meshes. Although there are references to other (Godunov flux) format numerical fluxes, but display expression is not given. The Harten-Lax-van Leer with contact (HLLC) has gathered interest because of their accuracy, mathematical simplicity, inherent entropy satisfying property, positivity, lack of demand for knowledge of complete eigenstructure of flux Jacobians, ability to satisfactorily handle shocks and expansion fans and their easy extensibility to various hyperbolic systems of governing equations [33–35]. The HLLC scheme is one of the simplest known Riemann solver that can resolve both linearly degenerate and genuinely nonlinear wavefields accurately. Therefore, it is very significant for researchers to find a display expression for HLLC flux on curved meshes.

In recent years, open source solvers have also become very popular. Such as famous open source libraries libMesh [36] and deal.II [37]are both excellent, which have brought inspiration and convenience to many researchers. However, Local Lax-Friedrichs (LLF) flux was used in deal.II [37, 38]. LLF has a larger numerical viscosity, which is not as applicable as HLLC in practical engineering such as supersonic and transonic. In [39], Qiu also point out that the HLLC fluxes might be good choices as fluxes for the RKDG method when all factors such as the cost of CPU time, numerical errors and resolution of discontinuities in the solution are considered. The libMesh is more of a software framework than a complete solver. In practical problems, libMesh has almost no complete functional modules that can be directly applied, and researchers need to do secondary development according to the specific problems. Therefore, the development of three-dimensional solver of aerodynamics simulation is very meaningful work.

In the present work, for aerodynamics simulation, a simulation tool based on curvilinear DG method will be provided to include both three-dimensional solid walls boundary conditions on curved meshes and HLLC flux on curved meshes. First, we present the relevant closed form expressions for a widely used set of hierarchical basis functions up to full high order for a reference element, as well as the necessary expressions for the Jacobian for a high order polynomial geometrical transformation. Given these results, implementation is then straightforward. The derived expressions are used to implement the curved meshes, which are then applied to integration along curved element and curved face for complex geometries. Second, we present the implementation process of HLLC flux in the curved element. Third, the implementation of solid wall boundary conditions depends on contravariant velocities which include the Jacobian of the transformation function of the map from each curved element to a rectilinear reference element. Based on the above three points, we have developed a simulation tool that can be used for solving the three-dimensional subsonic, transonic and hypersonic inviscid flows. Academic problem of subsonic flow past a sphere is used to show how the linear treatment of wall boundaries limits the accuracy, and demonstrate how the boundary treatment involving curved element overcomes this restriction. The results suggest that such a scheme is stable on a reasonably coarse meshes and finer ones. Engineering problems of transonic flow past an ONERA M6 wing, hypersonic flow past a blunt cone, and hypersonic flow past a HB-2 ballistic model are used to demonstrate the ability of this method for engineering practice. Numerical tests suggest that the curvilinear DG method has great robustness for three-dimensional inviscid flows with various geometries and velocities on unstructured meshes.

The paper is organized as follows: we first discuss the main idea of curvilinear DG method in Section 2 for three-dimensional inviscid flows. We present an implementation process of HLLC flux with curved element in Section 3. As a demonstration, we discuss the boundary treatment in Section 4. Numerical tests are shown in Section 5. Section 6 consists of concluding remarks.

2 Discontinuous Galerkin discretization

2.1 Three-dimensional compressible Euler equations

The governing equations for the three-dimensional inviscid flow can be written as follows:which is defined on a boundary domain with appropriate well-posed boundary data prescribed on . The conservation variables and the Cartesian components , , and of the flux function are given by

As is customary, is the density of the fluid, , , and are components of the velocity vector , is the pressure, and is the total internal energy per unit mass. The total enthalpy per unit mass is defined as . We assume that the fluid is perfect gas satisfying the equation of statewhere indicates the ratio between the specific heats of the fluid, which is set to 1.4 for the numerical experiments.

2.2 Discretization

The governing equation Eq. 1 is discretized using a DG finite element formulation which originally proposed in reference [40]. Here we provide a brief synopsis of the numerical scheme. For numerical discretization, we divide the problem domain into a collection of non-overlapping elements

Then, approximate solution is introduced on element and we construct a Galerkin problem by multiplying Eq. 1 with a test function , integrating over the domain , and using an integration by parts to obtainwhere is the outward unit normal vector to the boundary . The solution is approximated by the combination of shape functions,where is the shape function of the polynomial of degree and . Then the Eq. 5 becomes the following system of equations

Due to the discontinuous nature of the numerical solution, the normal flux is not defined on . The usual strategy is to define it in terms of a numerical flux that depends on the solution on and on the neighboring element sharing the portion of the boundary common to both elements. In our experiments, we use the HLLC flux which is necessary to consider the property of curved element discussed in Section 3. In Eq. 7, the method of calculating the integral of the volume and surface will be introduced in the following sections.

By assembling all the elemental contributions together, the system of ordinary differential equations which govern the evolution in time of the discrete solution can be written aswhere denotes the mass matrix and denotes the residual vector. The above system of ordinary differential equations is discretized in time by a three-stage TVD Runge-Kutta method [40].

2.3 Geometric mapping of curved element

In order to improve the representation of curved wall boundaries, elements with second-order shape can be used in the vicinity of the geometry. The geometric transformation of the real curved element in space to a rectilinear reference element in space, shown in Figure 1, is defined (via the Jacobian) in terms of second order Lagrangian interpolation functions. Figure 1B is a standard reference element and the lengths of the tetrahedron edges along the , and axes are unitary. Hence, it is apparent to define three of simplex coordinates , , on this tetrahedron [41, 42], and the last is then . The second order Lagrangian interpolation functions are constructed by , , , and identical to those used as nodal basis functions in finite elements. For the ten node tetrahedron, these transformation functions are defined by

FIGURE 1

The transformation from reference to real coordinates is then accomplished bywhere are the nodal coordinates in real space.

In order to define and manipulate vector quantities within a curved element, the unitary base vectors (co-variant) are introducedand the reciprocal base vectors are

Then, a vector can be represented by both its covariant components and the reciprocal base vectors

Differentiation in the reference coordinate system is connected to that in the real coordinate system as follows [22, 41]:where is the Jacobian of the transformation function and the inverse of the Jacobian can be written as

Then, Eq. 14 can be written aswhere , , and are the covariant components in the reference cell. Finally, the function Eqs 10, 11 are differentiated to obtain the Jacobian, which is shown in Table 1.

Although the expression is given, it is worth pointing out that the Jacobian is position-dependent within the element, and must be computed at each quadrature point when computing the integral of the volume and surface in the next subsection. It should be noted that the second order polynomial geometric transformation used here does not exactly model curved. More accurate geometrical modeling can be obtained by either increasing the order of the polynomial mapping or using tetrahedral rational Bezier volumes as shown in [43]. Although there are still some geometrical approximation errors, the fit is far better than that achieved using rectilinear elements.

2.4 Gradient operation of basis functions in reference coordinates

The basis functions and their gradient operation are derived for the parent element in terms of simplex coordinates. The actual set of basis functions used is the hierarchical type [44]. The hierarchical basis functions mean that the low order basis functions are maintained as a subset of the next higher-order basis functions. In order to derive the gradient on the reference element, the following associations are again made:

In the reference coordinates, , the reciprocal base vector, is simply , and similarly and . The hierarchical basis functions and their gradient operation complete to full third order are presented in Table 2.

2.5 Computation of integral of the volume and surface

The above subsections discussed geometric mapping of curved element. The real curved element (in the global coordinate system) is mapped to a rectilinear parent element (in the local ) coordinate system. Therefore, the integral of the volume and surface is changed, which should take into account the factors of the curved element.

In Eq. 8, the mass matrix elements are computed asand the residual vector elements are separated into the following two partswhere indicates the element, and indicates the neighboring element. In the reference coordinate system, considering Eqs 17, 21 becomes

In Eq. 22, is a HLLC schemes numerical flux which will be discussed in the next section. Considering Eq. 17, the calculation of has the following four situations.

In the case where , the Eq. 22 becomesin the case where , the Eq. 22 becomesin the case where , the Eq. 22 becomesand in the case where , the Eq. 22 becomeswhere is the Jacobian matrix of face, the superscript indicates the parameter in the reference element, and expressions will be discussed in the next section.

The integral of the volume and surface cannot be evaluated analytically and must be computed numerically using quadrature techniques [45]. The use of high order quadrature rules results in increased computational cost when using curvilinear as opposed to rectilinear elements. Therefore, it is advisable that curved elements only be used along the curved boundaries of the domain.

3 The implementation process of Harten-Lax-van Leer with contact flux in curved element

Amongst various approximate Riemann solvers that exist in the literature, the HLLC schemes have gained popularity because of their simplicity and accuracy. In this work, we intend to closely study the numerical schemes of HLLC flux in curved element. The original Harten-Lax-van Leer (HLL) scheme was devised by Harten, Lax and van Leer [46]. It assumes a wave structure consisting of two waves that separates three constant states. Although quite accurate in resolution of nonlinear waves, a major drawback of the HLL scheme is its inability to exactly resolve the contact and shear waves. The loss of accuracy on these waves occurs because of the assumption of constant average state between the two wave structure [34]. The inability of HLL solver to resolve contact and shear waves was mitigated through the development of the HLLC Riemann solver (C for Contact) by Toro et al. [47]. This improvement was achieved by considering two averaged intermediate states, separated by adding a contact wave with speed to the pre-existing two wave HLL structure.

The two-state approximate Riemann solution is defined as [46]

The corresponding interface flux, denoted , is defined aswhere , ; and are numerical approximations to the speeds of the left most and right most running characteristics that emerge as the solution of the Riemann problem at an interface.

In the rectilinear reference element, through Eq. 17, the unit normal vector is represented aswhere , and are the component of the unit normal vector of the interface in rectilinear reference element. Therefore, the directed velocity, in real curved element, can be rewritten as followwhere , and are contravariant velocities, which can be defined as

The choice of wavespeeds are given as [33]where , , ; and are the normal velocities across an interface; and are the respective sonic speeds; superscript “∼” indicates the standard Roe averaged quantities at the interface [48]. Applying the Rankine-Hugoniot conditions across the wave gives [33]

Similarly, the jump relations across the wave gives

For three-dimensional Euler equations with states and separated by an interface with unit normal vector , Eq. 34 becomes [33]

Similarly, Eq. 35 becomes

To determine , (and hence , ), Batten et al. [33] made the assumption thatwhere is the average directed velocity between the two acoustic waves. Multiplying the second, third, and fourth equations of Eq. 36 by , and respectively, and summing minus multiplying the second, third, and fourth equations of Eq. 37 by , and respectively, and summing, a closed form expression for the contact wavespeed, , can be written aswhere

Once is obtained, Eq. 36 can be manipulated to find all remaining components of . Using Eq. 38, the first equation in Eq. 36 gives

Multiplying the second, third, and fourth equations of Eq. 36 by , and respectively, and summing gives

Using Eq. 41 this simplifies to giveWith and specified, , , and may be obtained from Eq. 36 as

From Eqs 41–47, the flux vector may be obtained directly from Eq. 36. In the case where , the equations for follow by simply changing the or subscripts to and , respectively.

4 Boundary conditions with curved element

4.1 Mesh generation

When exploiting the ability of higher-order discretization methods to generate accurate approximations on coarse meshes, a crucial point is to provide a proper representation of curved wall boundaries [23]. Using a finer mesh to fit the boundary is a common method. However, inserting additional elements close to a boundary with the sole purpose of resolving the geometry impedes the aims of higher-order methods.

A more adequate approach we pursue within the proposed method is to represent the boundary by Simmetrix’ Modeling Suite [49], which is able to cope with high-order meshes. We gain detailed information on the geometry via additional point data included in the mesh.

The curved mesh generation process is briefly illustrated by an example in Figure 2. Firstly, solver read the model and mesh parameters, then call the mesh generation module, and use a certain method (such as Delaunay triangulation) to generate the first-order meshes, as shown in Figures 2A,D. According to geometry information, second-order meshes are then produced by projecting the mid-point of boundary edges on the geometry, resulting in curved elements on the boundary, as shown in Figures 2B and Figure 2E.

FIGURE 2

4.2 Boundary conditions

When the boundary face of an element belongs to , the normal flux function must be consistent with the boundary condition to be imposed on and will be denoted by , where is the internal boundary state at the current time level and is chosen according to the conditions that must be satisfied on the boundary. At inflow and outflow boundaries, the state is computed by means of a local one-dimensional characteristic analysis in a direction normal to the boundary by imposing the available data and the Riemann invariants associated to outgoing characteristics [26]. The most popular way to impose boundary conditions at solid wall is the reflection technique, where an extra row (rows) of ghost cells is added behind the wall. All interior solution components are reflected symmetrically to ghost states except for the normal velocity which is negated [27].

The reflecting boundary conditions state that no flow penetrates a solid wall, i.e., the normal velocity at the wall is zero. Depending on the numerical scheme, a ghost state or cell is created on the part of the numerical boundary corresponding to the solid wall. With the DG methods, a ghost state is created at every integration point on , where all components of the ghost solution are set equal to the corresponding interior values at the same point except for the normal velocity, which is negated. Then, the interior and ghost states are passed to a Riemann solver. Due to the symmetry of the reflection, the solution of the Riemann problem at integration point satisfies [27]where is velocity vector, is the inward unit normal vector to the boundary face of an element belongs to which is shown in Figure 3A. The geometric description of Eq. 48 is shown in Figure 3C.

FIGURE 3

A general process for implementing boundary conditions is shown below. A ghost state is created at each boundary integration point . ViaFigure 3C, by the addition of vectors, the relative to the element at the ghost state is given bywhere is the interior state, and is the normal velocity. The is defined aswhile the density and pressure are copied exactly from the interior values at the same point. We obtain the ghost state vector as

Finally, the Riemann problem is solved. From Eq. 51 we found that the ghost state vector at each integration point of the boundary face of an element belongs to is the same.

This approach works well for straight-sided bodies. However, results are inferior when a physical geometry is more complex [27]. Even more, the DG method is highly sensitive to the error due to approximation of a curved geometry by a straight-sided element meshes. This error may dominate the discretization error of the scheme and lead to a wrong solution. We seek to impose boundary conditions which would take this into consideration and model flow more accurately. When considering the factors of the curved element, we proceed as follow.

In Section 4.1, we know that the solid wall boundary region is discretized by curved elements. When the boundary face of a curved element belongs to , the normalized inward unit normal vector at each quadrature point is different which is shown in Figure 3B. Therefore, the solid wall boundary conditions, Eq. 51, cannot be directly applied. Even though the boundary conditions via a correct representation of the normal to the geometry is rather convenient to use for Eq. 51 [17, 27], the acquisition of normal is not very easy and needs to provide geometric information. Here, we show a direct method to implement solid wall boundary conditions which only requires some variables of rectilinear reference element.

Using the variables of rectilinear reference element, the normalized inward unit normal vector can be written aswhere

Using Eqs 52, 50 becomeswhere

From Eq. 49 we obtain the ghost state vector as

Since has a different value at each integration point, the ghost state vector at each quadrature point of the boundary face of an element belongs to is different.

5 Numerical examples

In order to demonstrate the performance of the proposed flow solver, we present several examples. All simulations were carried out on a Windows 10 64-bit Intel Core E5-2697 2.3-GHz and 256-GB RAM small workstation. For transonic flow and supersonic flow, the HWENO limiter in reference [50] is applied to the proposed flow solver. All computations were performed until solution coefficients reached a steady state, defined as the residual in reference [27].

5.1 Subsonic flow

We consider a subsonic flow at March number on four meshes having 10,670, 13,944, 271,327, 901,175 elements (the average spherical surface mesh sizes are 0.4, 0.2, 0.1, 0.05, defined as the ratio to the radius of the sphere). In order to check the accuracy and the convergence properties of the flow solver we have performed various computations using different combinations of interpolation functions for the unknowns and for the geometric mapping. Different elements are denoted by , where indicates the order of the complete polynomials used to approximate the unknowns and indicates the order of the complete polynomials used for the geometric mapping [26].

First, we solve the problem on the sequence of meshes with element and plot Mach distribution in Figure 4A. These solutions are very inaccurate, as put in evidence by the nonphysical boundary which develops along the solid wall and by the associated wake in the downstream region of the sphere. Notice that the solution obtained on the finest mesh is asymmetric and has a visible wake.

FIGURE 4

The picture changes completely when using , and elements, as reported in Figures 4B–D, which shows the Mach distribution computed on the four successively refined meshes. The quality of the solution clearly improves as and increases, which is asymmetric in the upstream-downstream direction. The solution corresponding to element on the coarsest mesh is similar to one obtained on the finest mesh with .

To quantify our findings, we measure errors in entropy defined aswhere and are pressure and density of the free stream, respectively. Table 3 reports the entropy error and the order of accuracy for the , , and computations. The convergence rate obtained by comparing the solutions of mesh and of mesh is computed aswhere is the mesh size. These computations allow us to appreciate the potentialities offered by second-order accurate elements in the numerical solutions of the three-dimensional Euler equations on unstructured mesh.

Further, we present two aerodynamic quantities: the pressure coefficient and the total pressure loss coefficient defined as a ratio of the total pressure at a point to the pressure of the free stream. The accuracy of the method is also evidenced by two aerodynamic quantities. The distribution of the pressure coefficient around the surface for four different meshes are shown in Figures 5A–D. The distribution of the total pressure loss coefficient around the surface for various meshes are shown in Figures 5E–H. We specify that mesh one to four refer to the coarse mesh to the fine mesh.

FIGURE 5

5.2 Transonic flow past ONERA M6 wing

A transonic flow over the ONERA M6 wing at a Mach number of and an attack angle of is considered in this example [51, 52]. This test case is chosen to demonstrate that the proposed flow solver with curvilinear DG method is able to enhance the accuracy of the underlying (rectilinear) methods for solving transonic flow problems and engineering problems. The mesh used in this computation consists of 682,726 elements, 122,679 points. As with the first example, indicates that the order of polynomials used to approximate the unknowns is one, and the order of the complete polynomials used for the geometric mapping is one too. The definition of is similar. In other words, is rectilinear element, and is curved element. The computed pressure distribution obtained by the DG solution on the wing surface are shown in Figure 6. In Figure 7 compares the pressure coefficient distributions at six span-wise locations (Figures 7A–F) on the wing surface between the numerical results and the experimental data. Under the same number of meshes, the pressure coefficient calculated by the curvilinear DG method is closer to the experimental value than the rectilinear DG method.

FIGURE 6

FIGURE 7

5.3 Hypersonic flow past blunt cone

This test case is selected to demonstrate that the proposed flow solver is able to solve hypersonic flow problems. The computation are performed for a high Mach number flow past the blunt cone model [53] using curvilinear and rectilinear DG method. The free stream Mach number, attack, pressure and temperature [53] for the blunt cone model correspond to , , and .

A three-dimensional axi-symmetric simulation is considered with the left and right boundaries defined as supersonic inflow and outflow boundaries respectively. The wall is modeled as a slip surface. The mesh used in this computation consists of 2,039,128 elements, 357,006 points. As with the previous definition, is rectilinear element, and is curved element. The computed density and pressure distribution near the blunted cone for the blunt cone model with is shown in Figure 8. The strong bow shock upstream of the body that is characteristic for this type of flow is well captured from our computations. The predicted wall pressure distribution non-dimensionalized with respect to the computed free stream pressure is compared with experimental data and shown in Figure 9. The computed pressure profile is slightly below predicted for a small region behind blunt cone, but in general the predicted pressure distribution is in good agreement with the experimental results. In Figure 9, we found that the results obtained by the curvilinear DG method are smoother, although the trends of wall pressure of curvilinear and rectilinear DG method are consistent.

FIGURE 8

FIGURE 9

5.4 Hypersonic flow past HB-2 ballistic model

Hypersonic flow past a HB-2 ballistic model is studied in this subsection. We carry out simulations at the flow conditions for which experimental setup details are available in the literature [54]. The free stream Mach number, attack, pressure and temperature [55, 56] for the HB-2 ballistic mode are , , and .

A three-dimensional axi-symmetric simulation is considered with the left and right boundaries defined as supersonic inflow and outflow boundaries respectively. The wall is modeled as a slip surface. The mesh used in this computation consists of 1,489,139 elements, 259,493 points. As with the previous definition, is rectilinear element, and is curved element. The computed density on the full computational domain as well as the pressure distribution near the blunted cone for the HB-2 ballistic mode with is shown in Figure 10. The strong bow shock upstream of the body as well as a weaker shock over the cylinder–cone juncture that is characteristic for this type of flow is well captured from our computations. The predicted wall pressure distribution non-dimensionalized with respect to the computed stagnation point pressure is compared with experimental data and shown in Figure 11. The computed pressure profile is slightly below predicted for a small region over the cylinder, but in general the predicted pressure distribution is in good agreement with the experimental results.

FIGURE 10

FIGURE 11

6 Conclusion

A flow solver based on curvilinear DG method has been presented for solving the three-dimensional subsonic, transonic and hypersonic inviscid flows on curved meshes. The method aims to present the relevant closed form expressions for a widely used set of hierarchical basis functions up to full high order for a reference element, as well as the necessary expressions for the Jacobian for a high order polynomial geometrical transformation. Given these results, the implementation of solid walls boundary conditions is very simple because it is based on the contravariant velocities and does not require a complex geometric boundary information. Based on covariant vectors, an improved algorithm of HLLC flux for solving the subsonic, transonic and hypersonic flows is suggested. A number of numerical experiments for a variety of flow conditions have been conducted to demonstrate the superior performance of this flow solver (curvilinear) over the underlying method (rectilinear). The solution of subsonic flow past a sphere indicates that such a flow solver is stable on a reasonably coarse mesh and finer ones, and can get accurate results at reasonably coarse tetrahedral mesh. The solution of transonic flow past an ONERA M6 wing, hypersonic flow past a blunt cone and hypersonic flow past a HB-2 ballistic model verify that our flow solver has engineering practice value. The simulation of subsonic, transonic and hypersonic flows show that the flow solver has great robustness for three-dimensional fluids with various geometries and velocities. In future work, the developed method will be applied to more complicated geometries, and extended to viscosity flows with curve boundary layer.

References

• 1.

  OktayEAlemdaroğluNTarhanEChampignyPEspineyP. Unstructured euler solutions for missile aerodynamics. Aerosp Sci Technol (2000) 4(7):453–61. 10.1016/s1270-9638(00)01072-5

  • CrossRef
  • Google Scholar

• 2.

  PanJLiuF. Wing flutter prediction by a small-disturbance euler method on body-fitted curvilinear grids. AIAA J (2019) 57(11):4873–84. 10.2514/1.j058120

  • CrossRef
  • Google Scholar

• 3.

  BeghamiWMaayahBBushnaqSArqubOA. The laplace optimized decomposition method for solving systems of partial differential equations of fractional order. Int J Appl Comput Math (2022) 8(2):52–18. 10.1007/s40819-022-01256-x

  • CrossRef
  • Google Scholar

• 4.

  ArqubOAHayatTAlhodalyM. Analysis of lie symmetry, explicit series solutions, and conservation laws for the nonlinear time-fractional phi-four equation in two-dimensional space. Int J Appl Comput Math (2022) 8:145. 10.1007/s40819-022-01334-0

  • CrossRef
  • Google Scholar

• 5.

  SrivastavaHMIqbalJArifMKhanAGasimovYSChinramR. A new application of gauss quadrature method for solving systems of nonlinear equations. Symmetry (2021) 13(3):432. 10.3390/sym13030432

  • CrossRef
  • Google Scholar

• 6.

  PourghanbarSManafianJRanjbarMAliyevaAGasimovYS. An efficient alternating direction explicit method for solving a nonlinear partial differential equation. Math Probl Eng (2020) 2020:1–12. 10.1155/2020/9647416

  • CrossRef
  • Google Scholar

• 7.

  ArqubOA. Numerical simulation of time-fractional partial differential equations arising in fluid flows via reproducing kernel method. Int J Numer Methods Heat Fluid Flow (2019) 30:4711–33. 10.1108/hff-10-2017-0394

  • CrossRef
  • Google Scholar

• 8.

  ArqubOA. Reproducing kernel algorithm for the analytical-numerical solutions of nonlinear systems of singular periodic boundary value problems. Math Probl Eng (2015) 2015:518406–13. 10.1155/2015/518406

  • CrossRef
  • Google Scholar

• 9.

  MascarenhasBSHelenbrookBTAtkinsHL. Application of P-multigrid to discontinuous Galerkin formulations of the euler equations. AIAA J (2009) 47(5):1200–8. 10.2514/1.39765

  • CrossRef
  • Google Scholar

• 10.

  AbbasAVicenteJDValeroE. Aerodynamic technologies to improve aircraft performance. Aerosp Sci Technol (2013) 28(1):100–32. 10.1016/j.ast.2012.10.008

  • CrossRef
  • Google Scholar

• 11.

  KouhiMOnateEMavriplisD. Adjoint-based adaptive finite element method for the compressible euler equations using finite calculus. Aerosp Sci Technol (2015) 46:422–35. 10.1016/j.ast.2015.08.008

  • CrossRef
  • Google Scholar

• 12.

  PonsinJFraysseFGomezMCordero-GraciaM. An adjoint-truncation error based approach for goal-oriented mesh adaptation. Aerosp Sci Technol (2015) 41:229–40. 10.1016/j.ast.2014.10.021

  • CrossRef
  • Google Scholar

• 13.

  DengXXieBXiaoF. Multimoment finite volume solver for euler equations on unstructured grids. AIAA J (2017) 55(8):2617–29. 10.2514/1.j055581

  • CrossRef
  • Google Scholar

• 14.

  PagliucaGTimmeS. Model reduction for flight Dynamics simulations using computational fluid Dynamics. Aerosp Sci Technol (2017) 69:15–26. 10.1016/j.ast.2017.06.013

  • CrossRef
  • Google Scholar

• 15.

  LiuYZhangWKouJ. Mode multigrid-a novel convergence acceleration method. Aerosp Sci Technol (2019) 92:605–19. 10.1016/j.ast.2019.06.001

  • CrossRef
  • Google Scholar

• 16.

  BottiLPietroDAD. Assessment of hybrid high-order methods on curved meshes and Comparison with discontinuous Galerkin methods. J Comput Phys (2018) 370:58–84. 10.1016/j.jcp.2018.05.017

  • CrossRef
  • Google Scholar

• 17.

  ToulorgeTDesmetW. Curved boundary treatments for the discontinuous Galerkin method applied to aeroacoustic propagation. AIAA J (2010) 48(2):479–89. 10.2514/1.45353

  • CrossRef
  • Google Scholar

• 18.

  BassiFBottiLColomboARebayS. Agglomeration based discontinuous Galerkin discretization of the euler and Navier-Stokes equations. Comput Fluids (2012) 61:77–85. 10.1016/j.compfluid.2011.11.002

  • CrossRef
  • Google Scholar

• 19.

  ZhangXTanS. A simple and accurate discontinuous Galerkin scheme for modeling scalar-wave propagation in media with curved interfaces. Geophysics (2015) 80(2):T83–T9. 10.1190/geo2014-0164.1

  • CrossRef
  • Google Scholar

• 20.

  ZhangX. A curved boundary treatment for discontinuous Galerkin schemes solving time dependent problems. J Comput Phys (2016) 308:153–70. 10.1016/j.jcp.2015.12.036

  • CrossRef
  • Google Scholar

• 21.

  CockburnBKarniadakisGEShuCW. Discontinuous Galerkin methods: Theory, computation and applications. New York: Springer Science & Business Media (2012).

  • Google Scholar

• 22.

  HesthavenJSWarburtonT. Nodal discontinuous Galerkin methods: Algorithms, analysis, and applications. New York: Springer Science & Business Media (2007).

  • Google Scholar

• 23.

  HartmannRHeldJLeichtTPrillF. Discontinuous Galerkin methods for computational aerodynamics-3d adaptive flow simulation with the dlr padge code. Aerosp Sci Technol (2010) 14(7):512–9. 10.1016/j.ast.2010.04.002

  • CrossRef
  • Google Scholar

• 24.

  AntoniettiPFCangianiACollisJDongZGeorgoulisEHGianiSet alReview of discontinuous Galerkin finite element methods for partial differential equations on complicated domains. Cham: Springe (2016).

  • Google Scholar

• 25.

  LuoHXiaYSpiegelSNourgalievRJiangZ. A reconstructed discontinuous Galerkin method based on a hierarchical weno reconstruction for compressible flows on tetrahedral grids. J Comput Phys (2013) 236:477–92. 10.1016/j.jcp.2012.11.026

  • CrossRef
  • Google Scholar

• 26.

  BassiFRebayS. High-order accurate discontinuous finite element solution of the 2d euler equations. J Comput Phys (1997) 138(2):251–85. 10.1006/jcph.1997.5454

  • CrossRef
  • Google Scholar

• 27.

  KrivodonovaLBergerM. High-order accurate implementation of solid wall boundary conditions in curved geometries. J Comput Phys (2006) 211(2):492–512. 10.1016/j.jcp.2005.05.029

  • CrossRef
  • Google Scholar

• 28.

  WarburtonT. A low-storage curvilinear discontinuous Galerkin method for wave problems. SIAM J Sci Comput (2013) 35(4):A1987–A2012. 10.1137/120899662

  • CrossRef
  • Google Scholar

• 29.

  ChanJHewettRJWarburtonT. Weight-adjusted discontinuous Galerkin methods: Curvilinear meshes. SIAM J Sci Comput (2017) 39(6):A2395–A421. 10.1137/16m1089198

  • CrossRef
  • Google Scholar

• 30.

  ChanJWilcoxLC. On discretely entropy stable weight-adjusted discontinuous Galerkin methods: Curvilinear meshes. J Comput Phys (2019) 378:366–93. 10.1016/j.jcp.2018.11.010

  • CrossRef
  • Google Scholar

• 31.

  MichoskiCChanJEngvallLEvansJA. Foundations of the blended isogeometric discontinuous Galerkin (bidg) method. Comput Methods Appl Mech Eng (2016) 305:658–81. 10.1016/j.cma.2016.02.015

  • CrossRef
  • Google Scholar

• 32.

  WangZJSunY. Curvature-based wall boundary condition for the euler equations on unstructured grids. AIAA J (2003) 41(1):27–33. 10.2514/2.1931

  • CrossRef
  • Google Scholar

• 33.

  BattenPClarkeNLambertCCausonDM. On the choice of wavespeeds for the hllc Riemann solver. SIAM J Sci Comput (1997) 18(6):1553–70. 10.1137/s1064827593260140

  • CrossRef
  • Google Scholar

• 34.

  SimonSMandalJC. A cure for numerical shock instability in hllc Riemann solver using antidiffusion control. Comput Fluids (2018) 174:144–66. 10.1016/j.compfluid.2018.07.001

  • CrossRef
  • Google Scholar

• 35.

  SimonSMandalJC. A simple cure for numerical shock instability in the hllc Riemann solver. J Comput Phys (2019) 378:477–96. 10.1016/j.jcp.2018.11.022

  • CrossRef
  • Google Scholar

• 36.

  KirkBSPetersonJWStognerRHCareyGF. Libmesh : A C++ library for parallel adaptive mesh refinement/coarsening simulations. Eng Comput (2006) 22(3-4):237–54. 10.1007/s00366-006-0049-3

  • CrossRef
  • Google Scholar

• 37.

  BangerthW. Deal.Ii (2021). Available from https://www.dealii.org/ (Accessed August 16, 2022).

  • Google Scholar

• 38.

  BangerthW. Deal.Ii the step-67 tutorial program (2021). Available from https://www.dealii.org/current/doxygen/deal.II/step_67.html (Accessed August 17, 2022).

  • Google Scholar

• 39.

  QiuJKhooBShuCW. A numerical study for the performance of the Runge-Kutta discontinuous Galerkin method based on different numerical fluxes. J Comput Phys (2006) 212(2):540–565. 10.1016/j.jcp.2005.07.011

  • CrossRef
  • Google Scholar

• 40.

  CockburnBShuCW. The Runge–Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems. J Comput Phys (1998) 141(2):199–224. 10.1006/jcph.1998.5892

  • CrossRef
  • Google Scholar

• 41.

  SwartzJPDavidsonDB. Curvilinear vector finite elements using a set of hierarchical basis functions. IEEE Trans Antennas Propag (2007) 55(2):440–6. 10.1109/tap.2006.888448

  • CrossRef
  • Google Scholar

• 42.

  FahsH. Improving accuracy of high-order discontinuous Galerkin method for time-domain electromagnetics on curvilinear domains. Int J Comput Math (2011) 88(10):2124–53. 10.1080/00207160.2010.527960

  • CrossRef
  • Google Scholar

• 43.

  MartiniEPelosiGSelleriS. A hybrid finite-element-modal-expansion method with a new type of curvilinear mapping for the analysis of microwave passive devices. IEEE Trans Microw Theor Tech (2003) 51(6):1712–7. 10.1109/tmtt.2003.812571

  • CrossRef
  • Google Scholar

• 44.

  YinJHXuLWangHXiePHuangSCLiuHXet alAccurate and fast three-dimensional free vibration analysis of large complex structures using the finite element method. Comput Struct (2019) 221:142–56. 10.1016/j.compstruc.2019.06.002

  • CrossRef
  • Google Scholar

• 45.

  ZhangLCuiTLiuH. A set of symmetric quadrature rules on triangles and tetrahedra. J Comput Math (2009) 27(1):89–96.

  • Google Scholar

• 46.

  HartenALaxPDLeerBV. On upstream differencing and godunov-type schemes for hyperbolic conservation laws. SIAM Rev (1983) 25(1):35–61. 10.1137/1025002

  • CrossRef
  • Google Scholar

• 47.

  ToroEFSpruceMSpearesW. Restoration of the contact surface in the hll-riemann solver. Shock Waves (1994) 4(1):25–34. 10.1007/bf01414629

  • CrossRef
  • Google Scholar

• 48.

  RoePL. Approximate Riemann solvers, parameter vectors, and difference schemes. J Comput Phys (1981) 43(2):357–72. 10.1016/0021-9991(81)90128-5

  • CrossRef
  • Google Scholar

• 49.

  SimmetrixI. Simmetrix’ modeling suite (2013). Available from http://www.simmetrix.com (Accessed August 6, 2021).

  • Google Scholar

• 50.

  LuoHBaumJDLöhnerR. A hermite weno-based limiter for discontinuous Galerkin method on unstructured grids. J Comput Phys (2007) 225(1):686–713. 10.1016/j.jcp.2006.12.017

  • CrossRef
  • Google Scholar

• 51.

  SchmittVCharpinF. Pressure distributions on the Onera-M6-Wing at transonic mach numbers. Exp Data base Comput Program Assess. Report of the Fluid Dynamics Panel Working Group 04, AGARD AR 138 (1979). Available from http://www.ttctech.com/samples/onera_m6_wing_multiblock/onera_mb.htm.

  • Google Scholar

• 52.

  BatinaJT. Accuracy of an unstructured-grid upwind-euler algorithm for the onera M6 wing. J Aircr (1991) 28(6):397–402. 10.2514/3.46040

  • CrossRef
  • Google Scholar

• 53.

  ZhongXMaY. Boundary-layer receptivity of Mach 7.99 flow over a blunt cone to free-stream acoustic waves. J Fluid Mech (2006) 556:55–103. 10.1017/s0022112006009293

  • CrossRef
  • Google Scholar

• 54.

  DamljanovićDRašuoBĐVMandićSIsakovićJ. Hypervelocity ballistic reference models as experimental supersonic test cases. Aerosp Sci Technol (2016) 52:189–97. 10.1016/j.ast.2016.02.035

  • CrossRef
  • Google Scholar

• 55.

  TisseraS. Assessment of high-resolution methods in hypersonic real-gas flows. Cranfield: Cranfield University (2010).

  • Google Scholar

• 56.

  JohnBEmersonDRGuXJ. Parallel Navier-Stokes simulations for high speed compressible flow past arbitrary geometries using flash. Comput Fluids (2015) 110:27–35. 10.1016/j.compfluid.2014.12.008

  • CrossRef
  • Google Scholar