Next Article in Journal
Layout Design of Strapdown Array Seeker and Extraction Method of Guidance Information
Previous Article in Journal
Analysis of Bilateral Air Services Agreement Liberalization in Australia
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

High-Order CFD Solvers on Three-Dimensional Unstructured Meshes: Parallel Implementation of RKDG Method with WENO Limiter and Momentum Sources

1
School of Aeronautics Science and Engineering, Beihang University, 37 Xueyuan Road, Haidian District, Beijing 100191, China
2
Ningbo Institute of Technology, Beihang University, Ningbo 315800, China
*
Author to whom correspondence should be addressed.
Aerospace 2022, 9(7), 372; https://doi.org/10.3390/aerospace9070372
Submission received: 9 June 2022 / Revised: 4 July 2022 / Accepted: 7 July 2022 / Published: 11 July 2022
(This article belongs to the Section Aeronautics)

Abstract

:
In aerospace engineering, high-order computational fluid dynamics (CFD) solvers suitable for three-dimensional unstructured meshes are less developed than expected. The Runge–Kutta discontinuous Galerkin (RKDG) finite element method with compact weighted essentially non-oscillatory (WENO) limiters is one of the candidates, which might give high-order solutions on unstructured meshes. In this article, we provide an efficient parallel implementation of this method for simulating inviscid compressible flows. The implemented solvers are tested on lower-dimensional model problems and real three-dimensional engineering problems. Results of lower-dimensional problems validate the correctness and accuracy of these solvers. The capability of capturing complex flow structures even on coarse meshes is shown in the results of three-dimensional applications. For solving problems containing rotary wings, an unsteady momentum source model is incorporated into the solvers. Such a finite element/momentum source hybrid method eliminates the reliance on advanced mesh techniques, which might provide an efficient tool for studying rotorcraft aerodynamics.

1. Introduction

Computational fluid dynamics (CFD) solvers, widely used in aerospace engineering, are based either on finite difference (FD) schemes or on finite volume (FV) schemes. However, FD schemes are restricted on structured meshes, while most FV schemes have only first- or second-order accuracy. One way to develop higher-order schemes suitable for unstructured meshes is to use the discontinuous Galerkin (DG) finite element (FE) methods. As with other FE methods, DG methods assume continuous approximation in each element (i.e., cell). However, they allow discontinuities to exist on intercell boundaries. Such discontinuities are then resolved by numerical fluxes, e.g., exact or approximate Riemann solvers, as with Godunov-type schemes [1]. The original DG scheme was not proposed for CFD but for solving the linear neutron transport equation [2]. For nonlinear conservation laws, Chavent and Salzano [3] designed an explicit Euler DG cheme, which is second-order accurate in space but only first-order accurate in time. To improve the accuracy of time discretization, Cockburn and Shu [4] replaced the first-order forward Euler method with high-order Runge–Kutta (RK) methods, which led to the family of Runge–Kutta discontinuous Galerkin (RKDG) schemes for scalar conservation laws [5,6], for one-dimensional conservation law systems [7], for multi-dimensional scalar conservation laws [8] and for multi-dimensional conservation law systems [9]. The last case includes the Euler system that depicts inviscid compressible flows, which is the first application of RKDG in CFD. We will give a concise derivation of this scheme in Section 2.1 and Section 2.2.
The Euler system is hyperbolic, and might have discontinuities in its solutions. Approximating such solutions with finite terms of continuous functions would generate numerical oscillations near discontinuities (known as the Gibbs phenomenon [10]). Therefore, any useful high-order scheme must incorporate some mechanism for suppressing such spurious oscillations. For second-order FD and FV schemes, there are slope and flux limiters [11]. These limiters work well on shock capturing but might reduce the order of accuracy when used in smooth regions. Furthermore, they are not applicable to higher-order schemes. To overcome these drawbacks, the essentially non-oscillatory (ENO) [12] and weighted ENO (WENO) [13] limiters for high-order FD schemes, which can essentially suppress spurious oscillations and maintain high-order accuracy in smooth regions, were developed. Similar ideas were later incorporated into FV [14] and DG [15,16] schemes, which are suitable for unstructured meshes. However, such early versions of unstructured WENO limiters have to query not only a cell’s immediate neighbors, but also the neighbors’ neighbors, and so on. To preserve the compactness of RKDG schemes, Zhong [17] and Zhu [18] proposed a compact WENO limter for two-dimensional structured and unstructrued meshes. In our early work [19], we extended it to three-dimensional and presented an efficient parallel implementation. We will summarize the key steps of this compact WENO limiter in Section 2.3.
The parallel implementation of the RKDG scheme with a compact WENO limiter can be used for solving many practical problems in aerospace engineering. However, for problems with moving boundaries (such as rotary wings on a helicopter or blades on a wind turbine), advanced meshing techniques (such as dynamic meshing, overset meshes or adaptive mesh refinement) have to be applied. On the other hand, typical lifting bodies (such as rotary wings) have been successfully modeled as various momentum sources in FD/FV schemes. Rajagopalan and Mathur [20,21] proposed the first FD/FV-based momentum source model, which models the time-averaged effect of a rotor as an actuator disk. Such a steady model was used by Kang and Sun [22,23] for solving the flow fields of rotors in ground effect. Shi [24] compared this model with an unsteady overset mesh method, and found that the steady momentum source model can predict complex aerodynamic interactions between a rotor and a flight deck at a much cheaper cost. Shen [25] and Kim [26] refined the model by replacing the actuator disk for the entire rotor with an actuator surface for each blade of the rotor, so that unsteady effects of rotary wings can be captured. Such early works are either oversimplified by ignoring unsteady effects (such as [20,21,22,23]), or restricted to structured meshes (such as [24,25]). To overcome these drawbacks, our parallel RKDG solvers are augmented with an unsteady actuator line model, which successfully captures unsteady effects of a rotor on unstructured meshes. The formulation of this finite element/momentum source hybrid method will be given in Section 2.4. Some demonstrative results are presented in Section 3.3.

2. The RKDG Method with WENO Limiters and Momentum Sources

Mathematically, the problems that we are going to solve in Section 3 are defined by the three-dimensional Euler system
t ρ ρ u x ρ u y ρ u z ρ e 0 + x ρ u x ρ u x u x + p ρ u y u x ρ u z u x ρ h 0 u x + y ρ u y ρ u x u y ρ u y u y + p ρ u z u y ρ h 0 u y + z ρ u z ρ u x u z ρ u y u z ρ u z u z + p ρ h 0 u z = 0 b x b y b z 0
with certain boundary and initial conditions, which describes inviscid compressible flows. In Equation (1),
  • Variables in the first column are called conservative variables;
  • The density ρ , the velocity components u x , u y , u z and the pressure p are called primitive variables;
  • The thermodynamic quantities e 0 and h 0 are the specific total energy and the specific total enthalpy, respectively; for ideal gases, they can be written as functions of the five primitive variables:
    e 0 = p / ρ γ 1 + u x 2 + u y 2 + u z 2 2 , h 0 = e 0 + p ρ ,
    in which γ is the heat capacity ratio of the gas, and we use γ = 1.4 in this article;
  • The body force components b x , b y , b z are assumed to be 0s, except in Section 3.3.
The method that we use for solving the Euler system (1) is the RKDG method with a WENO limiter. A detailed formulation of this method for sourceless problems has been given in our previous work [19]. In this section, we give a more complete formulation, which takes source terms into consideration. If all source terms were canceled, it would reduce to the form in [19].

2.1. The DG Space Discretization

Equation (1) can be written as
t U ̲ + · F ̲ = Q ̲ ,
in which F ̲ = F x ̲ e x + F y ̲ e y + F z ̲ e z is the flux vector whose components F x ̲ , F y ̲ , F z ̲ , and U ̲ , Q ̲ are all 5 × 1 matrices, each row of which is a scalar function depending on position x and time t. By multiplying both sides with an arbitrary test function ψ ( x ) and integrating the products on the i-th element (i.e., cell) E i , and applying integral by parts and Gauss’s divergence theorem, Equation (2) can be turned into a weak form:
E i ψ t U ̲ F ̲ · ψ + E i ν · F ̲ ψ = E i Q ̲ ψ ,
where ν is the outer normal unit vector of E i (which is the boundary of E).
By choosing an orthonormal basis of the linear space spanned by polynomials less than the p-th degree over E, denoted as ϕ ̲ ( x ) : = ϕ 1 ( x ) ϕ L ( x ) , the unknowns U ̲ and the arbitrary test function ψ can be approximated as
U ̲ ( x , t ) U ̲ h ( x , t ) = l = 1 L U ̲ ^ l ( t ) ϕ l ( x ) , ψ ( x ) ψ h ( x ) = l = 1 L ψ ^ l ϕ l ( x ) ,
in which each U ̲ ^ l ( t ) is a 5 × 1 matrix of temporal functions, and each ψ ^ l is an arbitrary constant number. Substituting them into the weak form (Equation 2) gives
l ψ ^ l k E i ϕ l ϕ k d U ^ ̲ k d t + E i F ̲ · ϕ l Q ̲ ϕ l + E i ϕ l F ν ̲ = O ̲ ,
where ν · F ̲ = : F ν ̲ is the normal flux at some point on E i , whose value could be given by an exact or approximate Riemann solver of Equation (2); see [27] for details.
Recall the arbitrariness of ψ ^ l l = 1 L and adopt the inner product notation
f | g : = E i f ( x ) g ( x ) ,
and we could turn Equation (4) into a system of ordinary differential equations:
d U ^ ̲ d t = R ̲ ( U ^ ̲ ) ,
in which
U ^ ̲ ( t ) = U 1 | ϕ 1 U 1 | ϕ L U 5 | ϕ 1 U 5 | ϕ L
is a 5 × L matrix of temporal functions (which will be solved in Section 2.2), and 
R ̲ = E i F x ̲ F y ̲ F z ̲ x ϕ ̲ y ϕ ̲ z ϕ ̲ Q ̲ ϕ ̲ E i F ν ̲ ϕ ̲
is a variable matrix depending on U ^ ̲ . The integrals in Equation (6) are evaluated by Gaussian quadrature rules, such as those given by [28].

2.2. The RK Time Discretization

Equation (5) is a typical system of nonlinear ordinary differential equations. In this article, we use the following explicit Runge–Kutta methods:
first-order: 
U ^ ̲ n + 1 = U ^ ̲ n + R ̲ n Δ t ;
second-order: 
U ^ ̲ n + 1 / 2 = U ^ ̲ n + R ̲ n Δ t , U ^ ̲ n + 1 U ^ ̲ n + 2 / 2 = 1 2 U ^ ̲ n + 1 2 U ^ ̲ n + 1 / 2 + R ̲ n + 1 / 2 Δ t ;
third-order: 
U ^ ̲ n + 1 / 3 = U ^ ̲ n + R ̲ n Δ t , U ^ ̲ n + 2 / 3 = 3 4 U ^ ̲ n + 1 4 U ^ ̲ n + 1 / 3 + R ̲ n + 1 / 3 Δ t , U ^ ̲ n + 1 U ^ ̲ n + 3 / 3 = 1 3 U ^ ̲ n + 2 3 U ^ ̲ n + 2 / 3 + R ̲ n + 2 / 3 Δ t ;
in which the R ̲ μ : = R ̲ ( U ^ ̲ μ ) s are the values mapped by the nonlinear operator R ̲ from the corresponding U ^ ̲ μ s.
It has been proven by Gottlieb and Shu that they are all total variation diminishing (TVD) [29] or strong stability preserving (SSP) [30] (which is a desired feature for hyperbolic problems), and fourth-order Runge–Kutta methods cannot be TVD or SSP without introducing the adjoint operator of R ̲ , denoted as R ̲ , which satisfies
U ^ ̲ n + 1 = U ^ ̲ n R ̲ ( U ^ ̲ n ) Δ t .
It is non-trivial to implement such an adjoint operator and thus we do not implement Runge–Kutta methods whose accuracy order is higher than three.

2.3. The Compact WENO Limiter

Following our previous work [19], we use the three-dimensional extension of the two-dimensional WENO limiters by Zhong [17] and Zhu [18]. For simplicity, we denote the index set of E i ’s neighbors as K i . For each k K i , we first transform Equation (2) to the ν -split form on the interface shared by E i and its neighbor E k :
t U ̲ + ν F ν ̲ t + A ν ̲ ν U ̲ = Q ̲ ,
in which the flux Jacobian can be approximated by the average value of U ̲ on E i :
A ν ̲ F ν ̲ U ̲ F ν ̲ U ̲ U ̲ i .
From the hyperbolicity of Equation (2), there always exists the eigenvalue decomposition
A ν ̲ = R ν ̲ u ν a u ν u ν u ν u ν + a ( R ν ̲ ) 1 ,
in which
R ν ̲ = 1 1 0 0 1 u x a ν x u x σ x π x u x + a ν x u y a ν y u y σ y π y u y + a ν y u z a ν z u z σ z π z u z + a ν z h 0 u ν a u x 2 + u y 2 + u z 2 2 u σ u π h 0 + u ν a , u ν u σ u π = ν x σ x π x ν y σ y π y ν z σ z π z 1 u x u y u z ,
( R ν ̲ ) 1 = 1 2 B 2 + u ν a 1 2 B 1 u x + ν x a 1 2 B 1 u y + ν y a 1 2 B 1 u z + ν z a 1 2 B 1 1 B 2 B 1 u x B 1 u y B 1 u z B 1 u σ σ x σ y σ z 0 u π π x π y π z 0 1 2 B 2 u ν a 1 2 B 1 u x ν x a 1 2 B 1 u y ν y a 1 2 B 1 u z ν z a 1 2 B 1 ,
where B 1 ( γ 1 ) / a 2 and B 2 : = B 1 ( u ν 2 + u σ 2 + u π 2 ) . The original conservative variable U ̲ is then projected into the space spanned by the columns of R ̲ , which defines the characteristic variable:
V ̲ : = ( R ν ̲ ) 1 U ̲ .
Each scalar component of V ̲ can now be treated as independent functions, to which any suitable scalar WENO limiter (such as the one in [17,18]) can be applied. Once obtaining the reconstructed characteristic variable V ̲ new , it is turned back to the original conservative variable:
U ̲ i | k new : = R ν ̲ V ̲ new ,
The subscript i | k emphasizes the fact that it is a function defined on E i with the help of E k . There is one such reconstructed U ̲ i | k new for each k K i , so the final step is to weight them by the volume of the corresponding adjacent element:
U ̲ i new : = k K i U ̲ i | k new | E k | k K i | E k | .

2.4. The Momentum Source Model

To model a fixed or rotary wing with a high aspect ratio, we use the actuator line model, which is a special case of the more general momentum source model. In this model, a wing is sliced into a series of thin pieces, each of which is perpendicular to the wing’s axis. When moving relative to the air, a piece feels a pair of aerodynamic forces:
L d s D d s = 1 2 ρ ( u 2 + w 2 ) C L ( α ) C D ( α ) c d s ,
in which
  • c is the chord length of the piece, and s is the arc length parameter of the axis;
  • α is the angle of attack, which is related to u and w (velocity components resolved in the airfoil’s local frame) by tan α = w / u ;
  • C L and C D are the lift and drag coefficients of the airfoil, respectively;
  • L and D are the lift and drag per unit length, respectively.
Using Newton’s third Law, the force per unit length felt by the air surrounding the piece is
b L ( r ( s ) ) = e L ( s ) e D ( s ) L ( s ) D ( s ) ,
in which e L and e D are unit vectors along and perpendicular to the direction of airflow relative to the airfoil. To avoid ambiguity, we denote the force per unit length as b L and use b V to denote the force per unit volume. If the intersection of a wing’s axis P Q ¯ and a DG element E is the line segment R T ¯ , then the body integral of the source term in Equation (6) is actually a line integral:
E i Q ̲ ϕ ̲ = E i 0 b V ( r ) 0 ϕ ̲ ( r ) = 0 ̲ R T ¯ b L ( r ) ϕ ̲ ( r ) 0 ̲ ,
By parameterizing the line segment R T ¯ , the line integral in Equation (13) can be evaluated as a definite integral:
R T b L ( r ) ϕ ̲ ( r ) = ξ R ξ T b L ( r ( ξ ) ) ϕ ̲ ( r ( ξ ) ) r ξ ξ ,
to which the standard Gaussian quadrature can be applied.

3. Results of Model Problems and Engineering Problems

In this section, we give the results of various problems to show the accuracy and performance of the methods described in Section 2.

3.1. Lower-Dimensional Model Problems

3.1.1. Shock Tube Problems

These problems are usually defined as one-dimensional problems, but we treat them as three-dimensional ones. All these problems are considered in a [ 0.0 , 5.0 ] × [ 0.0 , 1.0 ] × [ 0.0 , 0.5 ] box with all boundaries closed but the left and right ends open. Using the method described in [27], their exact solutions can be obtained, which can be used for assessing the accuracy of our solvers. In our earlier work [19], we have given the results obtained by running our solvers on an unstructured hexahedral mesh. This time, we use an unstructured tetrahedral mesh instead.
Problem 1
(Sod). Solve the Euler system Equation (1) for t [ 0.0 , 1.0 ] with the initial condition
ρ u v w p t = 0 = 1 0 0 0 1 , x < 2.5 ; 1 8 0 0 0 1 10 , x > 2.5 .
This is a classical problem of inviscid compressible flows. It contains all three types of discontinuities: a shock wave and a contact discontinuity running towards the right and an expansion wave running towards the left. In Figure 1, we plot the density contour given by our third-order solver with an EigenWeno limiter. In Figure 2, we compare the density distributions given by various solvers along the longitudinal axis (on which y = 0.5 and z = 0.25 ) of the box. The accuracy of our solvers and the effect of p-refinement are clear in these figures.
It is predictable that both mesh refinement (decreasing h) and order increment (increasing p) can help to improve accuracy. To compare the performance of solvers with different orders more fairly, it is better to use finer meshes for running lower-order solvers. After a few trials, we find that the solutions given by the h-p pairs listed in Table 1 roughly have the same level of accuracy, as shown in Figure 3. It is clear, at least for Problem 1, that high-order solvers are better than low-order ones in the sense of obtaining the same level of accuracy with less time and space costs. Similar conclusions were drawn in our previous work [19], in which the p = 3 solution of a linear advection problem on an h 1 / 4 mesh defeated the p = 1 solution of the same problem on an h 1 / 32 mesh in accuracy but saved quite a lot of time. These encouraging results justify our efforts to implement higher-order solvers.
Problem 2
(Lax). Solve the Euler system Equation (1) for t [ 0.0 , 0.6 ] with the initial condition
ρ u v w p t = 0 = 0.445 0.698 0.000 0.000 3.528 , x < 2.5 ; 0.500 0.000 0.000 0.000 0.571 , x > 2.5 .
This is another problem that involves all three types of discontinuities. It is more difficult than the previous one in the sense that its solution contains values beyond the range of initial values and the discontinuities are much steeper. As before, we plot the density contour of our third-order solution in Figure 4 and compare the density distributions given by various solvers in Figure 5. Both figures demonstrate the effect of p-refinement and the ability of our high-order solvers to suppress numerical oscillations.
Problem 3
(Vacuum). Solve the Euler system Equation (1) for t [ 0.0 , 0.4 ] with the initial condition
ρ u v w p t = 0 = 1 4 0 0 0.4 , x < 2.5 ; 1 + 4 0 0 0.4 , x > 2.5 .
This problem is easier than the previous two in the sense that no strong discontinuities (shock or contact discontinuity) occur in the solution. However, there is a region of vacuum, which does not occur in the previous problems, generated between the left- and right-running expansion waves. The velocity u inside the vacuumed region (denoted by a ∗ in subscript) is undefined, but its boundary values should be determined from the left (denoted by an L in subscript) and right (denoted by an R in subscript) initial values by the two Riemann invariants:
u * L = u L + 2 a L γ 1 , u * R = u R 2 a R γ 1 .
See [27] for more detailed discussions. Unfortunately, approximate Riemann solvers usually (have to) neglect this condition. For this reason, we use an exact Riemann solver to obtain fluxes on cell boundaries. As before, we plot the results in Figure 6 and Figure 7, which again validate the correctness and robustness of our solvers.

3.1.2. The Forward Step Problem

This is a classical two-dimensional problem [31], but we treat it as a three-dimensional one:
Problem 4.
Solve the Euler system (Equation 1) in a [ 0 , 3 ] × [ 0 , 1 ] × [ 0 , Z ] box (representing a wind tunnel), where the thickness Z could be any positive value, with a [ 0.6 , 3 ] × [ 0 , 0.2 ] × [ 0 , Z ] box removed (representing a forward-facing step). The x = 0 surface is open as an inlet and the x = 3 surface is open as an outlet. All other boundary surfaces are closed as solid walls. The initial condition is given as a uniform state:
ρ u v w p t = 0 = 1.4 3.0 0.0 0.0 1.0 .
To show the applicability of our solvers on structured meshes, we solve this problem on meshes such as the one in Figure 8. As a common practice, we plot the density contours at various time steps in Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19.
The flow starts from a uniform supersonic ( M = 3.0 ) state. In Figure 9 and Figure 10, a curved bow shock wave generates in front of the forward-facing step and a fan-shaped expansion wave generates at the corner ( x = 0.6 , y = 0.2 ). The bow shock then hits the top of the tunnel (Figure 11) and reflects from it (Figure 12 and Figure 13). At approximately t = 1.2 (Figure 14), the reflected shock hits the bottom of the tunnel (the top of the step) and reflects from it. Both of the two reflections are regular at these moments. Shortly before t = 2.0 (Figure 14), a Mach stem that characterizes a Mach reflection is generated from the first reflection point. A wavy contact discontinuity generating from the triple point is captured by Figure 16, Figure 17, Figure 18 and Figure 19. This flow structure would be smeared off if we solved Problem 4 using lower-order solvers or coarser meshes.

3.2. Three-Dimensional Engineering Problems

In this section, we solve the Euler system Equation (1) on a three-dimensional unstructured mesh, which discretizes the space around a YF-17 aircraft by tetrahedral cells, as shown in Figure 20.

3.2.1. YF-17 in Subsonic Flight

Problem 5.
Solve the Euler system Equation (1) in the surroundings of a YF-17 aircraft. The initial condition is given as a subsonic flow ( M = 0.3 ) with a 20-deg angle of attack:
ρ u v w p t = 0 = 1.4 0.2819 0.0 0.1026 1.0 .
All boundaries of the aircraft are defined as solid walls, except the intake of the engine, which is defined as a subsonic outlet of the condition equal to the initial condition, and the exhaust of the engine, which is defined as a supersonic ( M = 3.0 ) inlet of the condition:
ρ u v w p ex . = 1.4 2.4 0.0 0.0 1.44 .
We solve this problem via two of our finite element solvers and plot the streamlines released from the strake and the wing at t = 10 in Figure 21 and Figure 22. It can be seen that the streamlines obtained from the high-order ( p = 3 ) solution are smoother than those from the low-order ( p = 1 ) solution, which shows the benefit of p-refinement. The difference in accuracy is more obvious in Figure 23, which clearly shows a more detailed flow structure in the third-order solution (left half) and the piecewise constantness of the first-order solution (right half). Both solutions are able to capture the vortex trailing from the wing tip, which is generated from the the pressure difference between the upper and lower surfaces of the wing. When the wing generates positive lift, the pressure on the lower wing surface is higher than that on the upper wing surface. Under the action of this pressure difference, the air under the wing rolls up around the tip and flows backward, and the tip vortex is thus formed. The vortex trailing from the strake is generated from a similar mechanism.

3.2.2. YF-17 in Supersonic Flight

Problem 6.
Solve the Euler system Equation (1) in the surroundings of a YF-17 aircraft. The initial condition is given as a supersonic flow ( M = 2.0 ) with a 0-deg angle of attack:
ρ u v w p t = 0 = 1.4 2.0 0.0 0.0 1.0 .
All boundaries of the aircraft are defined as solid walls, except the intake of the engine, which is defined as a supersonic outlet of the condition equal to the initial condition, and the exhaust of the engine, which is defined as a supersonic ( M = 3.0 ) inlet of the condition:
ρ u v w p ex . = 1.4 2.4 0.0 0.0 1.44 .
We solve this problem using the same solvers as for Problem 5. Density contours on the y = 0 surface obtained from the first- and the third-order solutions are given in Figure 24 and Figure 25, respectively. It is obvious that shock waves are captured well (without spurious oscillations) by both of them. To show the difference more clearly, we compare the two solutions in Figure 26 and Figure 27. As in the subsonic case, the third-order solution (which is piecewise quadratic) outperforms the first-order one (which is piecewise constant). Shock waves (red) are generated on surfaces facing the wind, since the relative speed of air is larger than the speed of sound. Expansion waves (blue) are generated on leeward surfaces, since the solid body contracts there, which leaves more room for the supersonic flow.

3.2.3. Parallel Efficiency

Since the mesh in Figure 20 is highly unstructured and non-uniform, simple geometric partitioning cannot achieve a relatively balanced distribution of computational effort. In our code, we use the Metis library [32] to partition the dual graph of the mesh, which gives the results in Figure 28 and Figure 29. The fluctuation of cell numbers is under 2 % , which is quite good since the optimal partitioning of an unstructured mesh is an NP-hard problem. With such an approximately optimal partitioning, the parallel efficiency, which is defined as
E = T serial P T parallel × 100 % ,
in which P is the number of processes, could surpass 99 % in theory. In practice, however, the E values given by our numerical experiments are only around 80 % , as shown in Figure 30, which is derived from the measured time costs given in Table 2. The main reason for the gap between theory and practice is that in the derivation of the ideal E value, inter-process communications are assumed to overlap in time with inner-cell computations, which is an over-optimistic assumption. Nevertheless, the acceleration is still significant, which reduces the wall clock time to an acceptable level.

3.3. Problems with Momentum Sources

In this section, we solve two problems of rotorcraft aerodynamics using the momentum source model described in Section 2.4. The problems are set to simulate wind tunnel tests of a rotor (a pair of rotary wings), using the mesh shown in Figure 31. All of the boundaries of the outside box are solid walls, except he left and right ends, which are set as the inlet and outlet, respectively. The spherical region is the circumscribed sphere of the rotor, whose rotating axis could point in any direction. The smaller box enclosing the sphere is used for refining the mesh in the surrounding and downstream regions of the rotor.

3.3.1. A Climbing Rotor

Problem 7.
Solve Equation (1) in the surroundings of a rotor, whose rotating axis points in the direction of ( 1 , 0 , 0 ) . The initial condition is given as a uniform flow:
ρ u v w p t = 0 = 1.29 1.0 0.0 0.0 101325.0 ,
which is also the background state at the two ends of the wind tunnel.
Third- and first-order solutions at the same moment ( t = 1 ) are shown in Figure 32 and Figure 33, respectively, in which density contours and velocity directions are plotted on selected slices. Large-scale flow structures, such as the contraction of airflow near the rotor disk and the rolled-up airwake, are captured in both figures. However, finer details, such as the ripples generated from the rotor disk and the strong vortex ring in the downstream, are only visible in the third-order solution. The nearly axisymmetric flow structure is caused by the periodic movement of the rotary wings. The density of air immediately below the rotor (red) is greater than that above the rotor (blue), which means that the rotor is compressing the air flowing across it. According to Newton’s third law, the rotor must be experiencing a force in the opposite direction exerted by the air. This is where the rotor thrust comes from.

3.3.2. A Rotor in Forward Flight

Problem 8.
Solve Equation (1) in the surrounding of a rotor, whose rotating axis points in the direction of ( 0 , 0 , 1 ) . The initial condition is given as a uniform flow:
ρ u v w p t = 0 = 1.29 10.0 0.0 0.0 101325.0 ,
which is also the background state at the two ends of the wind tunnel.
As with the climbing problem, we plot the third- and first-order solution at the same moment ( t = 1 ) in Figure 34 and Figure 35, respectively. As before, both of them can capture large-scaled structures, such as the skewed and rolled-up airwakes. The rolled-up structure is more clear in Figure 36, in which streamlines are plotted and colored by the magnitude of velocity. The piecewise constant first-order solution is much vaguer and loses many details, which shows again the value of developing high-order solvers.

4. Discussion

In this article, we give a concise but complete formulation of an RKDG scheme with a compact WENO limiter for solving inviscid compressible flows, possibly with momentum sources, on three-dimensional unstructured meshes. The algorithms are implemented on top of public available libraries, which support unstructured mesh partitioning, inter-process communication and distributed memory parallelization. The correctness and accuracy of the solvers are validated by lower-dimensional reference problems. Results of real three-dimensional applications are also presented, in which complex flow structures are captured by the third-order solver even on very coarse unstructured meshes. Future works may include implementing higher-order solvers under this framework and testing them on larger meshes and on larger high-performance computing platforms.

Author Contributions

Conceptualization, W.P. and S.L.; methodology, W.P.; software, W.P. and Y.J.; validation, W.P. and Y.J.; formal analysis, W.P. and Y.J.; investigation, W.P.; resources, W.P.; data curation, W.P.; writing—original draft preparation, W.P.; writing—review and editing, W.P., Y.J. and S.L.; visualization, W.P.; supervision, S.L.; project administration, S.L.; funding acquisition, S.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National High-Tech R&D Program of China (863 Program), grant number 2012AA112201.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this article are all generated from the source code publicly available in our GitHub repository at https://github.com/pvc1989/miniCFD accessed on 6 July 2022 (or the mirror site https://gitee.com/pvc1989/miniCFD accessed on 6 July 2022), except the mesh used in Section 3.2, which was downloaded from https://cgns.github.io/CGNSFiles.html accessed on 13 May 2022.

Acknowledgments

The authors would like to express their deepest appreciation to Shen Xin for his professional guidance on computer science, and to Zhou Yukai for his patient assistance with mathtyping.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations and Nomenclature

The following abbreviations are used in this manuscript:
CFDComputational Fluid Dynamics
FDFinite Difference
FVFinite Volume
FEFinite Element
DGDiscontinuous Galerkin
RKRunge–Kutta
TVDTotal Variation Diminishing
SSPStrong Stability Preserving
ENOEssentially Non-Oscillatory
WENOWeighted ENO
t , x , y , z partial derivative with respect to t , x , y , z
α angle of attack
γ heat capacity ratio of the gas
ν outer normal vector at some point on E i
π , σ tangential vectors at some point on E i
ρ density of the gas
ϕ k k-th function in an orthonormal basis
ψ an arbitrary test function
ψ h approximation of ψ
ψ ^ k projection of ψ on ϕ k
A ν ̲ Jacobian of F ν ̲ with respect to U ̲
b L body force per unit length
b V body force per unit volume
b x , b y , b z x , y , z -component of body force per unit volume
cchord length of a blade
E i and E i i-th element (cell) and its boundary
e 0 specific total energy of the gas
F ̲ flux vector whose components are matrices
F ν ̲ projection of F ̲ on ν
F x ̲ , F y ̲ , F z ̲ x , y , z -component of F ̲
h 0 specific total enthalpy of the gas
ppressure of the gas
Q ̲ column matrix of source terms
R ̲ changing rate of U ^ ̲
R ν ̲ matrix whose j-th column is the j-th eigenvector of A ν ̲
U ̲ column matrix of conservative variables
U ̲ h and U ̲ new approximation of U ̲ and its WENO reconstruction
U ̲ i | k new WENO reconstruction of U ̲ h on E i with the help of E k
U ^ ̲ coefficient matrix whose k-th column is U ^ ̲ k
U ^ ̲ k projection of U ̲ on ϕ k
U ^ ̲ n and R ̲ n U ^ ̲ and R ̲ at the n-th time step
u x , u y , u z x , y , z -component of velocity
V ̲ column matrix of characteristic variables
V ̲ h and V ̲ new approximation of V ̲ and its WENO reconstruction

References

  1. Harten, A.; Lax, P.D.; van Leer, B. On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws. SIAM Rev. 1983, 25, 35–61. [Google Scholar] [CrossRef]
  2. Reed, W.H.; Hill, T.R. Triangular Mesh Methods for the Neutron Transport Equation; Technical Report LA-UR-73-479; Los Alamos Scientific Lab.: Los Alamos, MN, USA, 1973.
  3. Chavent, G.; Salzano, G. A finite-element method for the 1-D water flooding problem with gravity. J. Comput. Phys. 1982, 45, 307–344. [Google Scholar] [CrossRef]
  4. Cockburn, B.; Shu, C.W. The Runge-Kutta local projection P1-discontinuous-Galerkin finite element method for scalar conservation laws. In Proceedings of the 1st National Fluid Dynamics Conference, Cincinnati, OH, USA, 25–28 July 1988; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 1988. [Google Scholar] [CrossRef] [Green Version]
  5. Shu, C.W.; Osher, S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J. Comput. Phys. 1988, 77, 439–471. [Google Scholar] [CrossRef] [Green Version]
  6. Shu, C.W.; Osher, S. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II. J. Comput. Phys. 1989, 83, 32–78. [Google Scholar] [CrossRef]
  7. Cockburn, B.; Lin, S.Y.; Shu, C.W. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems. J. Comput. Phys. 1989, 84, 90–113. [Google Scholar] [CrossRef] [Green Version]
  8. Cockburn, B.; Hou, S.; Shu, C.W. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV. The multidimensional case. Math. Comput. 1990, 54, 545–581. [Google Scholar] [CrossRef] [Green Version]
  9. Cockburn, B.; Shu, C.W. The Runge–Kutta Discontinuous Galerkin Method for Conservation Laws V: Multidimensional Systems. J. Comput. Phys. 1998, 141, 199–224. [Google Scholar] [CrossRef]
  10. Gottlieb, D.; Shu, C.W. On the Gibbs Phenomenon and Its Resolution. SIAM Rev. 1997, 39, 644–668. [Google Scholar] [CrossRef] [Green Version]
  11. Sweby, P.K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws. SIAM J. Numer. Anal. 1984, 21, 995–1011. [Google Scholar] [CrossRef]
  12. Harten, A.; Engquist, B.; Osher, S.; Chakravarthy, S.R. Uniformly High Order Accurate Essentially Non-oscillatory Schemes, III. J. Comput. Phys. 1997, 131, 3–47. [Google Scholar] [CrossRef]
  13. Jiang, G.S.; Shu, C.W. Efficient Implementation of Weighted ENO Schemes. J. Comput. Phys. 1996, 126, 202–228. [Google Scholar] [CrossRef] [Green Version]
  14. Hu, C.; Shu, C.W. Weighted Essentially Non-oscillatory Schemes on Triangular Meshes. J. Comput. Phys. 1999, 150, 97–127. [Google Scholar] [CrossRef] [Green Version]
  15. Qiu, J.; Shu, C.W. Runge–Kutta Discontinuous Galerkin Method Using WENO Limiters. SIAM J. Sci. Comput. 2005, 26, 907–929. [Google Scholar] [CrossRef]
  16. Zhu, J.; Qiu, J.; Shu, C.W.; Dumbser, M. Runge–Kutta discontinuous Galerkin method using WENO limiters II: Unstructured meshes. J. Comput. Phys. 2008, 227, 4330–4353. [Google Scholar] [CrossRef]
  17. Zhong, X.; Shu, C.W. A simple weighted essentially nonoscillatory limiter for Runge–Kutta discontinuous Galerkin methods. J. Comput. Phys. 2013, 232, 397–415. [Google Scholar] [CrossRef]
  18. Zhu, J.; Zhong, X.; Shu, C.W.; Qiu, J. Runge–Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes. J. Comput. Phys. 2013, 248, 200–220. [Google Scholar] [CrossRef]
  19. Pei, W.; Jiang, Y.; Li, S. An Efficient Parallel Implementation of the Runge–Kutta Discontinuous Galerkin Method with Weighted Essentially Non-Oscillatory Limiters on Three-Dimensional Unstructured Meshes. Appl. Sci. 2022, 12, 4228. [Google Scholar] [CrossRef]
  20. Rajagopalan, R.G.; Mathur, S.R. Three dimensional analysis of a rotor in forward flight. In Proceedings of the 20th Fluid Dynamics, Plasma Dynamics and Lasers Conference, Buffalo, NY, USA, 12–14 July 1989; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 1989. [Google Scholar] [CrossRef]
  21. Rajagopalan, R.G.; Mathur, S.R. Three Dimensional Analysis of a Rotor in Forward Flight. J. Am. Helicopter Soc. 1993, 38, 14–25. [Google Scholar] [CrossRef] [Green Version]
  22. Kang, N.; Sun, M. Technical Note: Prediction of the Flow Field of a Rotor in Ground Effect. J. Am. Helicopter Soc. 1997, 42, 195–198. [Google Scholar] [CrossRef]
  23. Kang, N.; Sun, M. Simulated Flowfields in Near-Ground Operation of Single- and Twin-Rotor Configurations. J. Aircr. 2000, 37, 214–220. [Google Scholar] [CrossRef]
  24. Shi, Y.; Xu, Y.; Zong, K.; Xu, G. An investigation of coupling ship/rotor flowfield using steady and unsteady rotor methods. Eng. Appl. Comput. Fluid Mech. 2017, 11, 417–434. [Google Scholar] [CrossRef] [Green Version]
  25. Shen, W.Z.; Zhang, J.H.; Sørensen, J.N. The Actuator Surface Model: A New Navier–Stokes Based Model for Rotor Computations. J. Sol. Energy Eng. 2009, 131, 011002. [Google Scholar] [CrossRef]
  26. Kim, T.; Oh, S.; Yee, K. Improved actuator surface method for wind turbine application. Renew. Energy 2015, 76, 16–26. [Google Scholar] [CrossRef]
  27. Toro, E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar] [CrossRef]
  28. Zhang, L.; Cui, T.; Liu, H. A Set of Symmetric Quadrature Rules on Triangles and Tetrahedra. J. Comput. Math. 2009, 27, 89–96. [Google Scholar]
  29. Gottlieb, S.; Shu, C.W. Total variation diminishing Runge-Kutta schemes. Math. Comput. 1998, 67, 73–85. [Google Scholar] [CrossRef] [Green Version]
  30. Gottlieb, S.; Shu, C.W.; Tadmor, E. Strong Stability-Preserving High-Order Time Discretization Methods. SIAM Rev. 2001, 43, 89–112. [Google Scholar] [CrossRef]
  31. Woodward, P.; Colella, P. The numerical simulation of two-dimensional fluid flow with strong shocks. J. Comput. Phys. 1984, 54, 115–173. [Google Scholar] [CrossRef]
  32. Karypis, G.; Kumar, V. A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs. SIAM J. Sci. Comput. 1998, 20, 359–392. [Google Scholar] [CrossRef]
Figure 1. Third-order solution of ρ ( t = 1.0 ) in Problem 1.
Figure 1. Third-order solution of ρ ( t = 1.0 ) in Problem 1.
Aerospace 09 00372 g001
Figure 2. Comparison between solutions of ρ ( t = 1.0 ) in Problem 1.
Figure 2. Comparison between solutions of ρ ( t = 1.0 ) in Problem 1.
Aerospace 09 00372 g002
Figure 3. Comparison between solutions given by various hp pairs.
Figure 3. Comparison between solutions given by various hp pairs.
Aerospace 09 00372 g003
Figure 4. Third-order solution of ρ ( t = 0.5 ) in Problem 2.
Figure 4. Third-order solution of ρ ( t = 0.5 ) in Problem 2.
Aerospace 09 00372 g004
Figure 5. Comparison between solutions of ρ ( t = 0.5 ) in Problem 2.
Figure 5. Comparison between solutions of ρ ( t = 0.5 ) in Problem 2.
Aerospace 09 00372 g005
Figure 6. Third-order solution of ρ ( t = 0.4 ) in Problem 3.
Figure 6. Third-order solution of ρ ( t = 0.4 ) in Problem 3.
Aerospace 09 00372 g006
Figure 7. Comparison between solutions of ρ ( t = 0.4 ) in Problem 3.
Figure 7. Comparison between solutions of ρ ( t = 0.4 ) in Problem 3.
Aerospace 09 00372 g007
Figure 8. A coarse ( h = 1 / 20 ) mesh for solving Problem 4. This mesh is too coarse to capture details of the flow field but clear enough to demonstrate the distribution of nodes and cells. Actually, we use a much finer ( h = 1 / 200 ) mesh to produce Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19.
Figure 8. A coarse ( h = 1 / 20 ) mesh for solving Problem 4. This mesh is too coarse to capture details of the flow field but clear enough to demonstrate the distribution of nodes and cells. Actually, we use a much finer ( h = 1 / 200 ) mesh to produce Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19.
Aerospace 09 00372 g008
Figure 9. Third-order solution of Problem 4 at t = 0.2 on an h = 1 / 200 mesh.
Figure 9. Third-order solution of Problem 4 at t = 0.2 on an h = 1 / 200 mesh.
Aerospace 09 00372 g009
Figure 10. Third-order solution of Problem 4 at t = 0.4 on an h = 1 / 200 mesh.
Figure 10. Third-order solution of Problem 4 at t = 0.4 on an h = 1 / 200 mesh.
Aerospace 09 00372 g010
Figure 11. Third-order solution of Problem 4 at t = 0.6 on an h = 1 / 200 mesh.
Figure 11. Third-order solution of Problem 4 at t = 0.6 on an h = 1 / 200 mesh.
Aerospace 09 00372 g011
Figure 12. Third-order solution of Problem 4 at t = 0.8 on an h = 1 / 200 mesh.
Figure 12. Third-order solution of Problem 4 at t = 0.8 on an h = 1 / 200 mesh.
Aerospace 09 00372 g012
Figure 13. Third-order solution of Problem 4 at t = 1.0 on an h = 1 / 200 mesh.
Figure 13. Third-order solution of Problem 4 at t = 1.0 on an h = 1 / 200 mesh.
Aerospace 09 00372 g013
Figure 14. Third-order solution of Problem 4 at t = 1.2 on an h = 1 / 200 mesh.
Figure 14. Third-order solution of Problem 4 at t = 1.2 on an h = 1 / 200 mesh.
Aerospace 09 00372 g014
Figure 15. Third-order solution of Problem 4 at t = 2.0 on an h = 1 / 200 mesh.
Figure 15. Third-order solution of Problem 4 at t = 2.0 on an h = 1 / 200 mesh.
Aerospace 09 00372 g015
Figure 16. Third-order solution of Problem 4 at t = 2.4 on an h = 1 / 200 mesh.
Figure 16. Third-order solution of Problem 4 at t = 2.4 on an h = 1 / 200 mesh.
Aerospace 09 00372 g016
Figure 17. Third-order solution of Problem 4 at t = 2.8 on an h = 1 / 200 mesh.
Figure 17. Third-order solution of Problem 4 at t = 2.8 on an h = 1 / 200 mesh.
Aerospace 09 00372 g017
Figure 18. Third-order solution of Problem 4 at t = 3.2 on an h = 1 / 200 mesh.
Figure 18. Third-order solution of Problem 4 at t = 3.2 on an h = 1 / 200 mesh.
Aerospace 09 00372 g018
Figure 19. Third-order solution of Problem 4 at t = 4.0 on an h = 1 / 200 mesh.
Figure 19. Third-order solution of Problem 4 at t = 4.0 on an h = 1 / 200 mesh.
Aerospace 09 00372 g019
Figure 20. A zoomed-in view of the most important boundaries of the mesh for solving problems in Section 3.2.
Figure 20. A zoomed-in view of the most important boundaries of the mesh for solving problems in Section 3.2.
Aerospace 09 00372 g020
Figure 21. Streamlines given by the first-order solution of Problem 5.
Figure 21. Streamlines given by the first-order solution of Problem 5.
Aerospace 09 00372 g021
Figure 22. Streamlines given by the third-order solution of Problem 5.
Figure 22. Streamlines given by the third-order solution of Problem 5.
Aerospace 09 00372 g022
Figure 23. Comparison between the third-order solution (left half) and the first-order solution (right half) of Problem 5.
Figure 23. Comparison between the third-order solution (left half) and the first-order solution (right half) of Problem 5.
Aerospace 09 00372 g023
Figure 24. Density contour obtained from the first-order solution of Problem 6.
Figure 24. Density contour obtained from the first-order solution of Problem 6.
Aerospace 09 00372 g024
Figure 25. Density contour obtained from the third-order solution of Problem 6.
Figure 25. Density contour obtained from the third-order solution of Problem 6.
Aerospace 09 00372 g025
Figure 26. Comparison between the third-order solution (upper half) and the first-order solution (lower half) of Problem 6 on the z = 0 surface (top view).
Figure 26. Comparison between the third-order solution (upper half) and the first-order solution (lower half) of Problem 6 on the z = 0 surface (top view).
Aerospace 09 00372 g026
Figure 27. Comparison between the third-order solution (upper half) and the first-order solution (lower half) of Problem 6 on the z = 0 surface (bottom view).
Figure 27. Comparison between the third-order solution (upper half) and the first-order solution (lower half) of Problem 6 on the z = 0 surface (bottom view).
Aerospace 09 00372 g027
Figure 28. A 100-part partitioning of the mesh used for solving problems in Section 3.2.
Figure 28. A 100-part partitioning of the mesh used for solving problems in Section 3.2.
Aerospace 09 00372 g028
Figure 29. Distribution of cells in the mesh partitioning given in Figure 28.
Figure 29. Distribution of cells in the mesh partitioning given in Figure 28.
Aerospace 09 00372 g029
Figure 30. Parallel efficiency of the third-order solver for solving Problem 6.
Figure 30. Parallel efficiency of the third-order solver for solving Problem 6.
Aerospace 09 00372 g030
Figure 31. A schematic mesh for solving problems in Section 3.3.
Figure 31. A schematic mesh for solving problems in Section 3.3.
Aerospace 09 00372 g031
Figure 32. Third-order solution of Problem 7 at t = 1 .
Figure 32. Third-order solution of Problem 7 at t = 1 .
Aerospace 09 00372 g032
Figure 33. First-order solution of Problem 7 at t = 1 .
Figure 33. First-order solution of Problem 7 at t = 1 .
Aerospace 09 00372 g033
Figure 34. Third-order solution of Problem 8 at t = 1 .
Figure 34. Third-order solution of Problem 8 at t = 1 .
Aerospace 09 00372 g034
Figure 35. First-order solution of Problem 8 at t = 1 .
Figure 35. First-order solution of Problem 8 at t = 1 .
Aerospace 09 00372 g035
Figure 36. Streamlines of the solution in Figure 34.
Figure 36. Streamlines of the solution in Figure 34.
Aerospace 09 00372 g036
Table 1. Comparison between time costs and file sizes of various hp pairs.
Table 1. Comparison between time costs and file sizes of various hp pairs.
hp#Cells#StepsTime CostFile Size
1/5012,469,67515001492.16 s1.23 GB
1/182117045500279.119 s213 MB
1/10319,815250201.034 s89.1 MB
Table 2. Time costs of the third-order solver running on different numbers of cores.
Table 2. Time costs of the third-order solver running on different numbers of cores.
P T 100 T 399 T 400 P T 399 T 100 399 100 P T 400 T 100 400 100
110,511.342,466.542,593.0106.873106.939
52389.259740.139769.96122.925123.012
101246.375086.625101.95128.436128.519
15839.5543420.173430.75129.462129.560
20645.7372632.322640.57132.882132.989
40325.2691310.721319.73131.833132.595
60222.890923.043931.754140.499141.773
80168.862672.439681.445134.736136.689
100137.787543.081552.682135.550138.298
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Pei, W.; Jiang, Y.; Li, S. High-Order CFD Solvers on Three-Dimensional Unstructured Meshes: Parallel Implementation of RKDG Method with WENO Limiter and Momentum Sources. Aerospace 2022, 9, 372. https://doi.org/10.3390/aerospace9070372

AMA Style

Pei W, Jiang Y, Li S. High-Order CFD Solvers on Three-Dimensional Unstructured Meshes: Parallel Implementation of RKDG Method with WENO Limiter and Momentum Sources. Aerospace. 2022; 9(7):372. https://doi.org/10.3390/aerospace9070372

Chicago/Turabian Style

Pei, Weicheng, Yuyan Jiang, and Shu Li. 2022. "High-Order CFD Solvers on Three-Dimensional Unstructured Meshes: Parallel Implementation of RKDG Method with WENO Limiter and Momentum Sources" Aerospace 9, no. 7: 372. https://doi.org/10.3390/aerospace9070372

APA Style

Pei, W., Jiang, Y., & Li, S. (2022). High-Order CFD Solvers on Three-Dimensional Unstructured Meshes: Parallel Implementation of RKDG Method with WENO Limiter and Momentum Sources. Aerospace, 9(7), 372. https://doi.org/10.3390/aerospace9070372

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop