1 Introduction

Large eddy simulation (LES) and direct numerical simulation (DNS) of compressible flows are of great interest nowadays. There has been extensive research activity towards more realistic simulation tools, including numerical discretization, physical modeling, temporal integration, and parallel computing, among which high order method has been a lasting hot topic. Various high order methods have come out in recent years, which are generally divided into three families. The first is the finite difference method (FDM), which is widely used nowadays thanks to its excellent efficiency and accuracy. Typical high order FDMs consist of the explicit difference method, the compact difference method, and the hybrid ones of these two, to name just a few^[1-4]. Major restrictions for these methods are their sensitivity to mesh quality, and infeasibility for unstructured mesh. Another family is the finite volume method (FVM)^[5]. One of the desirable feature of the FVM is their feasibility to unstructured mesh. However, to achieve higher order accuracy, FVM would generally require a quite wide stencil which would make their application less practical. There has been effort aimed to relieve this problem^[6]. The third choice is the so-called high order unstructured schemes, including discontinuous Galerkin (DG)^[7], spectral difference (SD)^[8], spectral volume (SV)^[9], and flux reconstruction (FR)^[10]. For one thing, these methods have received great focus recently due to their desirable features in terms of accuracy, spectral properties, compactness, flexibility for unstructured mesh, and scalability. For the other, they are still undergoing rapid development in various aspects such as shock capturing, scalable implicit time integration, mesh adaptivity, and parallel computing.

For high order unstructured methods such as DG, shock capturing combined with high order approximation is still very challenging. Among the common shock capturing methods are limiting, reconstruction, and artificial viscosity models. For limiting, the solution slope of the troubled element would be limited by some limiter^[11] and the corresponding accuracy would be decreased to ensure stabilization. Although the method is quite robust, it will tend to destroy the high order accuracy. Several attempts have been made to alleviate this problem^[12-13]. The reconstruction method introduces the concept of weighted essentially non-oscillatory (WENO)^[14] or Hermite WENO (HWENO)^[15-16], and reconstructs the approximation polynomial on elements containing discontinuity. Compared with WENO, HWENO-based methods are able to effectively achieve the compact stencil, and are thus more suitable for unstructured mesh. Another large family of shock capturing for unstructured high order methods is the artificial viscosity method, which introduces explicit dissipation to smooth discontinuities and has achieved considerable success. The key for the artificial viscosity method is the right scaling, i.e., to decide the appropriate value of a shock sensor to activate the artificial viscosity. There has been plenty of work in this area^[17-19].

In this work, the performance of the dilation-based shock sensor combined with high order DG methods is investigated. The idea of using velocity dilation as shock sensor is not new. Cook and Cabot^[20] proposed to use the high order derivative of strain stress for under-resolved structures using compact difference schemes. Later Fiorina and Lele^[21], Kawai and Lele^[22], and Mani et al.^[23] constantly improved the model for LES of compressible flows. Premasuthan et al.^[24] extended this model to the spectral difference method. However, it has been found that the computation of high order derivatives can be both expensive and inaccurate especially for unstructured mesh. Moro et al.^[19] derived an artificial viscosity model based directly on the dilation. Through both the survey in the literature and our experience, the velocity dilation is quite a robust shock sensor. It is not our intention to propose a novel method in this paper. Rather, the idea of dilation is revisited to examine its accuracy and robustness for a range of benchmark problems.

The aspects which may set this work apart from those in the literature include: (ⅰ) The artificial viscosity is determined to be proportional to the velocity dilation (instead of its higher order derivatives), and is combined with the Laplacian term, which is very simple, while in the literature most dilation methods employ the Navier-Stokes viscous term. (ⅱ) A method to compute the scale of curvilinear elements is proposed. (ⅲ) A smoothing method using global C0 projection for artificial viscosity is proposed. The paper is organized as follows: governing equation and numerical discretization are briefly described in Section 2. The artificial viscosity model used in this paper is given in Section 3. Several typical benchmark test cases are illustrated in Section 4, followed by conclusions in Section 5.

2 Governing equation and discretization 2.1 Governing equation

The conservation system with diffusion terms is given as

(1)

where q is the conserved variables, and f_c and f_v are the convective and diffusive flux terms, respectively. In this paper, only the Euler system will be considered. Therefore, f_v only serves as numerical diffusion, and employs the Laplacian formulation given as

(2)

where μ_β denotes the artificial viscosity, which will be computed using the method given in Section 3.

2.2 Discretization method

The DG solver employed in this paper is developed within the framework of Nektar++^[25], to which the interested reader is referred. The underlying numerics is briefly described as follows.

For DG, the computational domain is decomposed into non-overlapping elements denoted as e, and the approximated solution on e is defined to be

(3)

where P is the order of the approximation polynomials, φ_m(x) represents the basis function chosen to be a modified hierarchical Jacobi polynomial basis^[25], and is the corresponding coefficient. The weak form of the discontinuous Galerkin formulation can be written as

(4)

where and are the inviscid and viscous flux respectively. In this work, Harten-Lax-van Leer-contact (HLLC) is used for the inviscid flux, and local discontinuous Galerkin (LDG) for the viscous one.

3 Artificial viscosity method 3.1 Dilation-based artificial viscosity

The dilation-based artificial viscosity is given as

(5)

where C_β is the user-defined parameter, u_e is the velocity vector, and h_e is the element size. C_β is fixed to be 1.5 in this paper unless specified otherwise. The superscript of μ_{β, e}^* denotes a raw value with piecewise discontinuity, and remains to be smoothed (described in the following section). Note that in (5) density is not included since the Laplacian formulation is used here, while for the physical viscous term, the addition of density would then be necessary. Also, by (5) we apply artificial viscosity on regions of expansion waves, from which test cases with strong wave interactions would benefit.

3.2 Method of dealing with curvilinear elements

From (5), it can be seen that the element size plays an important role for artificial viscosity, and an inappropriate evaluation of h_e would frequently cause μ_{β, e}^* to vary away from the intended value, especially for curvilinear elements. In this section, an estimation method of h_e is proposed for quadrilateral elements.

As shown in Fig. 1, the four edges of a quadrilateral element are denoted by E0, E1, E2, and E3 respectively. Also, the midpoint position vector and length of E0 are defined to be C0 and L0 respectively, and similar definitions apply to all the other edges as well. Then the average scale in the ξ and η directions is estimated to be

(6)

Next, the two scales h_ξ and h_η are compared. Assuming that h_ξ > h_η, we define

(7)

(8)

If h_ξ < h_η, n_η is defined first, and n_ξ is then computed through rotating n_η. With the above definition, the two symmetric scale vectors of the element can be given as

(9)

Finally, the estimated scale for a quadrilateral element can be determined by

(10)

where , and ε is a small constant to prevent division by zero.

For a triangular element, the size of the bounding box constraining the element is employed. The scale of the element in the x and y directions are given as

(11)

where "max" and "min" denote the maximum and minimum coordinates of the corresponding directions of the bounding box. The final scale of the element would be

(12)

It should be noted that the above method is not optimal for triangular elements of large aspect ratio. However, we argue that the method is reasonable since it is generally desirable that the aspect ratio of a triangular element is close to 1.

3.3 Smoothing viscosity

It has been noted that the raw values of (5) would produce viscosity which is discontinuous at the element interface. It is generally considered important to further smooth the viscosity to make the method more robust and accurate^[26], and a C0 continuation of the artificial viscosity is generally sufficient. A common method is to employ a piecewise linear reconstruction enforcing continuity at the element edge. In this section, we propose an alternative method, which achieves C0 continuity using the continuous Galerkin formulation.

Given the original field value, denoted as u^*(x), its C0 approximate solution u_h(x) is determined through solving the following projection problem:

(13)

where , and ϕ_m(x) denotes both the trial and test functions. Also, note that a C0 continuity condition at the element edges is applied to the left hand side of (13).

(13) can then be written as

(14)

where

(15)

(16)

(17)

And M_el and M_er are determined according to the element connectivity. (14) can then be solved using standard methods for linear systems. Note that the above smoothing procedure is suitable for arbitrary element type. Also, it is easy to see that (13) actually achieves the minimum error compared with u^*(x) in the sense of energy norm^[27].

4 Results

In this section, several benchmark cases are performed to evaluate the performance of the artificial viscosity model.

4.1 Convergence test of inviscid cylinder flow

This case is conducted to examine the effect of artificial viscosity on the accuracy of smooth flows. The free stream Mach number is set to be 0.1. Three sets of successively refined meshes (32×8, 64×16, and 128×32) are employed. For the baseline mesh (32×8), the distance of the first layer off the wall is 0.1d, and the exterior boundary extends to 10d, with d being the diameter of the cylinder. The other meshes are successively refined from the baseline one. Entropy production is used to indicate the numerical errors.

The convergence result of the entropy error versus the mesh resolution is plotted in Fig. 2. For comparison purpose, the exact third and fifth order slopes as well as results without artificial viscosity are given as well. As can be seen, those without artificial viscosity for both P2 and P4 agree well with the corresponding exact slopes. Once the artificial viscosity model is turned on, detrimental effects on accuracy are observed as expected. For P2, the error level increases, while the convergence rate remains unaffected; while for P4, both the error level (still lower than P2) and convergence rate are significantly affected. And for P4, a third order convergence rate can still be achieved. This can actually be obtained from (5). For smooth flow, |∇·u_e| would be a finite value, and therefore the error would scale as O(h_e²). Similar results are given by Moro et al.^[19], in which numerical error scale as O(h_e²) before the shock sensor therein is turned off as the mesh refines. Although this seems to make the higher order DG method with dilation-based artificial viscosity less accurate, we argue that the flow region near the wall contributes most of the numerical error, and in a practical case, a switch or a wall function^[28] to reduce the artificial viscosity near the wall would be able to improve the accuracy of high order methods. But this topic is not included in this paper.

The density contours for P2 with and without artificial viscosity are compared in Fig. 3 and Fig. 4. It can be seen that the results are comparable at least in a visual norm. This comparison is similar for P4, and thus omitted.

4.2 Stationary shock

In order to examine the robustness of the artificial viscosity method for shocks of various strengths, the one-dimensional stationary shock is computed. A linear distribution between the two states of a stationary shock is specified as the initial condition, which is then advanced in time till a steady solution is obtained. The final position of the shock is ambiguous due to ill-posedness of the problem. Following Moro et al.^[19], the origin of the plotted results is set to be the position where the value of the density equals the average of the two states before and after the shock. The computational domain in the flow direction is chosen to be [-1, 1]. And both one-dimensional and two-dimensional cases are performed.

4.2.1 One-dimensional case

For this case, the element size h is fixed to be 0.08. Three Mach numbers (i.e., 1.4, 5, and 10) are computed, and the artificial viscosity model combined with P2, P4, and P7 is employed.

The computational results including both density and artificial viscosity profiles are shown in Figs. 5-7, where x_c on the horizontal axis denotes the coordinate where density achieves the average of the two states ahead of and behind the shock. All the results with different shock strengths are effectively stabilized. It is worth noting that the same set of empirical parameters is used for all the computation, indicating good robustness of the artificial viscosity model. It can be seen that the shock width steepens with the increasing of approximation order, and for P7 the shock is captured within one element. Also, the artificial viscosity is well restricted, and decreases significantly away from the shock. In order to gain a better intuition, the original value of the artificial viscosity is plotted throughout this paper, and shows large variance among different cases due to various field values as will be seen in the following. We note that this is reasonable since the comparison on the value of the artificial viscosity between different cases makes little sense due to their great distinction. On the contrary, we are more concerned with the correspondence of the artificial viscosity with typical flow structures, such as shocks, contact discontinuities, and small fluctuations.

Furthermore, in order to examine the performance of the method, the density profiles of different shock strengths are compared in Fig. 8. All the results are smoothed with similar width, independent of the shock strength.

4.2.2 Two-dimensional unstructured mesh

The same problems as the above section are performed using P4 on a two-dimensional unstructured mesh as shown in Fig. 9, and the horizontal and vertical sizes of the mesh element are both equal to 0.08. Results for shock strengths of Ma=1.4, Ma=5, and Ma=10 are shown in Figs. 10-12. It can be seen that the shocks are smeared to a width similar to the one-dimensional results. Also, the contours of artificial viscosity are smoothed well indicating the effectiveness of the C0 projection methods on unstructured mesh.

4.3 Lax shock-tube problem

The Lax shock-tube problem is conducted in this section. The computational domain is taken to be [0, 1], and the element size is chosen to be 0.01. The initial condition is

(18)

The computation is stopped at t=0.13. Density and velocity profiles are given in Fig. 13, in which the reference denotes the exact solution to this problem. It can be seen that results by all the approximation orders agree well with the exact solution, and the spurious oscillations are well suppressed. The slight overshooting at the tail of the expansion wave is attributed to insufficient viscosity of the model for contact discontinuity at the start-up phase. We further note that this model may encounter stability issues for problems with strong contact discontinuities. The superior accuracy of higher order methods over low ones can also be observed even for this shock-dominated case. Note that although this dilation-based shock sensor cannot trigger any artificial viscosity for contact discontinuity, as can be seen in Fig. 14, the results show little spurious oscillations around that region.

4.4 Shu-Osher problem

The Shu-Osher problem is the one-dimensional model of shock/turbulence interaction, which is quite challenging since the entropy waves could be easily over-dissipated to ensure the strong shocks stabilized. This case is conducted here to examine the potential of the artificial viscosity model for shock/turbulence interaction problems. The computational domain is [-5, 5], and the element size is uniform and chosen to be 1/15. The initial condition is given as follows:

(19)

The computation is stopped at t=1.8. The density profiles are given in Fig. 15. The reference solution is obtained on a mesh of 4 000 grid cells with a fifth order WENO scheme. It can be seen that approximation order shows a major effect on this complex problem. The entropy waves are preserved significantly better with higher order approximation. The glitch around x=-2.8 is attributed to the start-up error. Artificial viscosity of different approximation orders is plotted in Fig. 16. Similar to the Lax problem, artificial viscosity corresponds well to the flow region where stabilization is needed.

It is worth noting that many artificial viscosity models in the literature result in piecewise-constant artificial viscosity which are then smoothed into piecewise linear distribution. This could possibly cause over-dissipation for delicate structures like the entropy waves in this case. The embedded sub-cell variation of the dilation-based artificial viscosity is able to adapt the flow in a sub-cell level better, although this remains to be further validated.

4.5 Double Mach problem

This case is computed to test the performance of the artificial viscosity model for two-dimensional strong shocks. The computational domain is [0, 4]×[0, 1]. The initial condition is an oblique shock with a Mach number of 10, making an angle of 60° with the horizontal line and intersecting with the x-axis at x=1/6. The computation is stopped at t=0.2.

Density results are given for approximation orders of P2 and P4 in Fig. 17, Fig. 18, and Fig. 19. Three sets of meshes are employed here. One is a uniform structured mesh with the element size being h=1/30, while the other two are uniform unstructured meshes obtained through dividing each element of a uniform structured mesh (h=1/102, h=1/204) into two triangles. It can be seen that the spurious oscillations are effectively suppressed, and the resolution is significantly improved with higher accuracy orders. For finer meshes, delicate structures behind the reflected shock are resolved. Also, the corresponding artificial viscosity to Fig. 17 is plotted in Fig. 20. For this two-dimensional case, the artificial viscosity is shown to be generally smooth, and is well restricted around shocks. Note that the different ranges in which the value of the artificial viscosity vary, indicating a much smaller numerical dissipation required by higher order methods.

4.6 Supersonic cylinder flow

In this section, the supersonic flow over cylinder is carried out to investigate the capability of the method for curvilinear meshes. A configuration of half cylinder is employed, and the mesh used is shown in Fig. 21, which contains 32 elements along the wall, and the minimum distance of the first layer of the grid points off the wall surface is 0.05d (d being the diameter of the cylinder). The free-stream Mach number is taken to be 3 and 8, and the computations have been conducted using P4.

Contours for density, Mach number, and artificial viscosity are shown in Fig. 22 (Ma=3) and Fig. 23 (Ma=8), respectively. No significant spurious oscillation has been observed. It is worth noting that the artificial viscosity varies smoothly on this curvilinear grid, which is critical for robustness.

5 Conclusions

The dilation-based artificial viscosity model is investigated combined with high order DG methods. Aspects to improve accuracy and robustness of the model are addressed including the estimation of element size and smoothing. Several benchmark test cases are conducted ranging from subsonic flows to hypersonic flows involving strong shocks. For smooth flows, the artificial viscosity is observed to cause detrimental effects for accuracy, which is expected to be relieved with various techniques such as a switch or a wall function. For cases containing shocks, the artificial viscosity model is able to suppress spurious oscillations well, and the shock is observed to be resolved within one cell for P7. The proposed method for estimating the element size is shown to be able to provide appropriate scale for a quadrilateral element of large aspect ratio. And despite being more expensive compared with a piecewise linear method, the proposed C0 projection smoothing method is able to achieve a smooth distribution of artificial viscosity, which is more helpful for unstructured mesh.

A most desirable feature of the dilation-based sensor is its robustness as well as high accuracy. All the cases in this paper are computed with a single and fixed parameter. We do not adjust the parameter to achieve optimum results for each case. Rather, we are interested in using one single constant to handle various problems, which would be more practical from a user's point of view. Also, the high order feature of the dilation-based model is able to adapt the flow better in a sub-cell level. This may be more important when it comes to turbulence flows, which are rich in delicate structures. According to the above analysis, the dilation-based model serves as a good candidate for compressible flows. However, we note that this model might encounter stability issues for strong contact discontinuities, and produce significant dissipation within regions of fluctuations. Therefore, modifications are apparently required to make the model more practical with the probability of introducing more empirical parameters, which constitute our future work.