2014, Vol. 35
Shanghai University
Article Information
- Zhen-hua JIANG, Chao YAN, Jian YU. 2014.
- High-order discontinuous Galerkin solver on hybrid anisotropic meshes for laminar and turbulent simulations
- Appl. Math. Mech. -Engl. Ed., 35(7): 799-812
- http://dx.doi.org/10.1007/s10483-014-1834-9
Article History
- Received 2013-1-1;
- in final form 2013-8-1
During the last decade,interest in the use of the discontinous Galerkin (DG)[1, 2] methods for compressible flow computations has become more widespread in aerodynamic applications. In the DG methods,the uniform high-order accuracy is obtained using the high degree polynomial approximation within elements as in the classical finite element methods,and the physics of wave propagation is accounted by solving the Riemann problems at the element interfaces as in the upwind finite volume methods. The DG methods have many desirable features,such as the capacity to handle complicated geometries,the compact stencils,the flexibility for the h/p adaptation,and the nice mathematical properties with respect to conservation,stability,and convergence.
However,the implementation of DG can be still subject to many limitations. The develop- ment of optimal or near optimal solution strategies for higher-order operators remains one of the key determining factors in devising DG methods which are not just competitive but superiorto other lower-order methods in terms of overall accuracy and efficiency[3]. In particular,for high Reynolds number flows anisotropic meshes are highly clustered to resolve the thin bound- ary layer,inducing additional stiffness and degrading the convergence performance of common point-implicit solvers dramatically. Fidkowski et al.[4] represented an elemental-line smoother to alleviate the stiffness associated with stretched grids,which was also carried out in [5, 6, 7, 8] for high-order spatial discretizations. The line-implicit solvers are attractive because they result in block-tridiagonal matrices which can be solved very efficiently. However,in the line-implicit methodology,off-diagonal Jacobian matrices need to be stored,and special techniques need to be developed to construct lines for the unstructured (mixed) meshes[4, 9].
Moreover,it is well known that the development of high-order methods for the solution to the Reynolds-averaged Navier-Stokes (RANS) equations is a challenging task because of the extreme stiffness related to the turbulence closure equations. While most efforts have been de- voted to robust implementation of one-equation or two-equations turbulence models in the high order discretizations[10, 11, 12, 13],less attention has been paid to the extension of algebraic turbulence model for the DG methods[13],especially on the unstructured meshes. In practice,algebraic turbulence models are relatively inexpensive to compute and have demonstrated generally su- perior accuracy and reliability for high-Reynolds number attached flows[14].
These two observations have motivated the group to develop a novel high-order DG solver with the capability to perform two-dimensional (2D) laminar and turbulent flow simulations. The distinguishable features of the present method include the application of the line-implicit algorithms based on structured/unstructured hybrid meshes to reduce the memory gain and meanwhile to guarantee the efficiency for highly anisotropic meshes; the extension of the Baldwin-Lomax (BL) model for the DG methods on the hybrid meshes with a simplified and ef- fective approach proposed to obtain the eddy viscosity distribution of the unstructured meshes; and the employment of an improved compact Hermite weighted essentially non-oscillatory (HWENO) limiter[15] to treat flow solution discontinuities.
2 Governing equationsSince the BL turbulence model is used in the present work,for both laminar and turbulent computations,the governing equations are
where ρ is the density,u is the velocity vector,p is the static pressure,e is the internal energy,T is the temperature,κ is the thermal conductivity coefficient,τ is the viscous stress tensor,and δ is the Kronecker tensor. This set of equations is completed by the addition of the equation of state which is valid for perfect gas,where is the ratio of the specific heats. 3 DG spatial and temporal discretization on hybrid meshes 3.1 DG formulation
To formulate the DG method,Eqs. (1)-(3) are rewritten as
where q is the conservative state vector,s is the auxiliary variable,fc(q) is the inviscid flux tensor,and fv(q,s) is the viscous flux tensor. For numerical discretization,the approximate solution qh,sh and the piecewise polynomial test function vh are introduced on the element K. The weak formulation of Eqs. (5)-(6) can then be given as where n is the outward unit normal vector to the boundary. The flux fc · n is replaced by a monotone numerical flux function,and the flux fv ·n is approximated by the second Bassi-Rebay scheme[16]. We express qh as
Quadrilateral elements are more suitable in tackling viscous boundary layers and more flex- ible for use with curved geometries in anisotropic regions of the mesh. Therefore,the present work utilizes high aspect-ratio quadrilaterals in the boundary layer zones around the solid sur- faces and then fills the major portion of the domain with isotropic triangular cells to handle the complex geometry. Furthermore,we treat quadrilateral cells and triangular cells separately by connecting their interfaces with the same distribution of flux points. The advantages of this method are quite clear.
High-order elements with curved boundaries can only be used to quadrilaterals,and those triangles attached to the surfaces of quadrilaterals increase the computational efficiency without mapping the entire mesh to high-order and meantime avoid inconsistent mesh cross-overs for the anisotropic cells.
One does not require line creation process for the line-implicit method in the current hybrid meshes because coordinate lines of the structured meshes are readily available for the aniso- tropic regions and point-implicit method can be applied for the isotropic regions filled with the unstructured meshes.
It facilitates the implementation of algebraic turbulence model and ensures the accuracy of the boundary layer where the structured quadrilateral elements are clustered.
3.3 Turbulence modeling peculiaritiesIn the BL model,the distribution of eddy viscosity is calculated along a constant index line. Hence,it is straightforward to extend the BL model for the DG methods in the context of structured meshes. High order eddy viscosity solutions can be used to resolve the boundary layer and wake regions. However,the implementation of the BL model on the unstructured meshes has not been entirely straightforward and usually involved interpolation of variables between the unstructured meshes and the background structured meshes[14]. Here,we propose a different approach to obtain the eddy viscosity values on the unstructured meshes for the DG methods. The recipe[17],which has been designed to improve the robustness of the BL model, is also adopted in our approach with small modifications using
instead of the entropy s used in Ref. [17],the purpose of which is to make the recipe less dependent on the free flow conditions,in particular,on the Mach number Ma. However,here, we note that accurate and reliable identification for the boundary layer edge is still an open issue especially for those complex flows. The final method can then proceed as follows.
(i) Construct background structured meshes and further equally divide each of them into N new quadrilateral cells,where N is the order of discretization (for example,see Fig. 1 for 2nd-order DG).
 |
| Fig. 1 Illustration of building index lines in hybrid meshes for 2nd-order DG scheme. |
(ii) Build a set of lines for identifying the boundary layer edge with Eq. (10) based on the location of cell centroid of unstructured meshes in the whole (undivided) background structured meshes and the location of the cell centroid of structured meshes. As shown in Fig. 1,the one resulting line is · · · 0495· · · .
(iii) Build another set of lines for executing the BL model based on the location of cell centroid of unstructured meshes in each new quadrilateral and the location of all quadraturepoints of structured meshes. As shown in Fig. 1,the two resulting lines are · · · 12345· · · and · · · 6789· · · .
(iv) Obtain the eddy viscosity values for cell centroids of unstructured meshes and all quadra- ture points of structured meshes by performing the modified BL model along the index lines.
The procedure only utilizes the background mesh to locate the unstructured mesh and still performs the unstructured computation in order to less weaken the merit of the unstructured meshes. Moreover,the formed index lines in the unstructured meshes are consistent with the volume quadrature point locations in the structured meshes. Assuming that the eddy viscosity is constant over the triangular cells,which has proven to be robust and accurate in Ref. [12], the proposed approach is actually a simple and effective way to determine the distribution of the eddy viscosity for the unstructured meshes.
Despite the order of discretization,only one line is built in each cell for identifying the boundary layer edge. Furthermore,only volume quadrature points are used in searching Ymax for those index lines inside the structured meshes. We remark that these are two important components in utilizing the modified BL model robustly and accurately for the DG methods, especially for relatively complex flow problems.
3.4 Temporal discretization and HWENO limiters on hybrid meshesUsing the backward Euler difference and linearizing the residual,the fully linearized formu- lation of Eq. (9) can be given as
where We use a symmetric Gauss-Seidel strategy with k iterations and express Eq. (11) as follows.
The forward Gauss-Seidel iteration is
and the backward Gauss-Seidel iteration is where D,L,and U represent diagonal,lower,and upper matrices,respectively. One can observe that off-diagonal matrices L and U are required only in the form of matrix-vector product,and such an operation can be approximated by computing increments of the flux vector. Thus,we are able to avoid forming off-diagonal matrices by the following finite difference approximation: where " is a small parameter[18]. The resultant scheme now involves a dual iteration strategy with high-order residual evaluation. To further reduce the computational cost,we choose an unfixed inner sweep strategy that we increase 1 inner sweep each time when the nonlinear residuals reach a drop of 2 orders of magnitude. The technique is quite simple and easy to implement in practical DG codes.
This modified elemental iterative method is referred to as the linearized Gauss-Seidel (LGS) method,which needs to sweep the grid elements in an orderly fashion. Meanwhile,to alleviatethe stiffness with stretched meshes,the line-implicit methodology is required to treat implicitly clusters of elements with strongest coupling,which also needs to perform a line creation process. The beauty of the line-implicit LGS method built upon the present hybrid meshes is that the orderings and coordinate lines are readily available for all anisotropic structured meshes. For isotropic unstructured meshes,the point-implicit method can then be applied. Therefore,the line creation process of the line-implicit methodology is avoided,and the structured and unstruc- tured meshes can be treated completely separately. The solver reduces the memory gain by only storing the block diagonal and off-diagonal matrices in the line and meantime guarantees the efficiency on the highly anisotropic meshes,which will be shown in the following section.
This work has also investigated the implementation of high-order HWENO limiters on the hybrid meshes to compute flows containing shock waves. We make a small modification of linear weight λ1 on the structured meshes by taking into account the element aspect-ratio (λAR) as follows:
which is observed to yield more robust and accurate results for the anisotropic meshes. 4 Numerical results 4.1 Laminar flow over flat plate
Firstly,we consider a laminar flow on an adiabatic flat plate characterized by a free-stream Mach number Ma∞ = 0.2 and a Reynolds number Re∞ = 105. For detailed validation of the developed method,we calculate this flow on three different meshes shown in Fig. 2. As can be seen,the hybrid mesh is composed of regular quadrilateral cells near the flat plate and irregular triangular cells for the remaining zone,and the unstructured mesh is built from the hybrid one with each quadrilateral divided into two triangles. The computed v-velocity profiles compared with the Blasius solution are plotted in Fig. 3. Clearly,for all meshes,the solutions get more accurate with increasing order of the discretization. Moreover,results by the 2nd-order DG scheme on the structured and hybrid mesh show improved agreement with Blasius solution than those on the unstructured mesh. It is noted that the hybrid case represents a 20% reduction of the total number of cells than the structured case,indicating that the hybrid mesh is capable of delivering equivalent accuracy compared with the structured mesh at a lower computational cost.
 |
| Fig. 2 Meshes used for computing laminar flow over flat plate. |
 |
| Fig. 3 vy profiles for laminar boundary layer problem obtained by 2nd- and 4th-order DG,where vy denotes normalized velocity component in y-direction. |
The laminar flow past a single NACA0012 airfoil is conducted with three different initial conditions.
Case 1 Ma∞ = 0.5,an angle of attack α = 0,and Re∞ = 5 000.
Case 2 Ma∞ = 0.8,α = 10,and Re∞ = 73.
Case 3 Ma∞ = 2.0,α = 10,and Re∞ = 106.
The effectiveness of the line-implicit solver is studied using two hybrid anisotropic meshes which consist of 3 098 elements with different maximum aspect-ratio of 30 and 300,respectively. Figure 4(a) shows the computational mesh with higher aspect-ratio cells,and Figs. 4(b)-4(d)
 |
| Fig. 4 Computational mesh and Mach number contours for laminar flow past NACA0012 airfoil problem. |
 |
| Fig. 5 Heat flux coefficient distributions along airfoil. |
 |
| Fig. 6 Efficiency and effectiveness comparison of line-implicit scheme and point-implicit scheme. |
A two-element airfoil test case taken from Ref. [8] is simulated to examine the robustness of the DG solver dealing with the complex configuration. The mesh used in the computation consists of 6 051 triangles and 3 132 quadrilaterals with a maximum λAR of 320. Figure 7 displays the computed Mach number contours obtained by the 4th-order DG scheme without any limiting or local order-reduction techniques. Note that the wake is very well preserved due to the anisotropic mesh and high-order accurate solutions. The density residuals depicted in Fig. 8 exhibit no rises through the convergence history and again a weak P-dependence on the highly anisotropic mesh,demonstrating the robustness and fast convergence performance of the current solver.
 |
| Fig. 7 Mach number contours for laminar flow past two-element airfoil problem. |
 |
| Fig. 8 Density residual convergence history forlaminar flow past two-element airfoilproblem. |
The case is selected to validate the implementation of the BL model and HWENO limiters on the hybrid mesh for transonic turbulent simulations. The flow conditions are Ma∞ = 0.734, α = 2.79◦,and Re∞ = 6.5×106. Computational meshes composed of 14 136 elements are shown in Fig. 9. The O-type background mesh is built from the structured mesh and clustered around the wake region to locate fine unstructured grids. Figure 10 presents the obtained Mach number contours using the 2nd,3rd,and 4th DG schemes. High-order schemes have better resolution for the captured shock wave,and HWENO limiters are indeed able to eliminate the spurious oscillations for all test cases. The surface pressure and skin friction distributions computed by 4th-order DG are displayed in Fig. 11,where we can see the accurate shock position capture and the structure of separation and reattachment flows induced by the transonic shock. Overall numerical solutions are smooth and compared favorably with the experiment[20],although there is a slightly low prediction for the leading edge suction peak[12].
 |
| Fig. 9 RAE2822 airfoil mesh zoom-in. |
 |
| Fig. 10 Mach number contours for transonic flow past RAE2822 airfoil problem. |
 |
| Fig. 11 Pressure coefficient distribution (left) and skin friction coefficient distribution (right) fortransonic flow past RAE2822 problem. |
High-lift multi-element airfoil flows involve strong but smooth flow field gradients which require a very robust implementation for both the turbulence model and solution strategies. Such challenging flows therefore provide a suitable test case to demonstrate the capability of our high-order DG solver to simulate complex turbulent flow problems. The numerical solutions are presented for Ma∞ = 0.2,α = 16.2◦,and Re∞ = 9 × 106,and the computational meshes in Fig. 12 contain 45 896 elements.
 |
| Fig. 12 Multi-element airfoil mesh zoom-in. |
Figure 13(a) presents the surface pressure distributions for the 3rd-order DG scheme,which demonstrates that the 3rd-order DG solutions provide a good agreement on the surface pressure compared to the experimental results[21]. The computed skin friction coefficients with the results of the NSU2D in Ref. [12] are plotted in Fig. 13(b). Overall,a relative smooth skin friction profile is obtained with the exception of small glitches due to the geometrical slope discontinuities,which is also slightly revealed with the NSU2D calculations. The Mach number contours are depicted in Fig. 14(a),where the slat wake persists all the way over the flap. Streamlines around the slat are illustrated in Fig. 14(b). With the high flow incidence angle and high streamline curvature over the gap,the flow is accelerated around the leading edge of the slat,and smooth reverse flow appears near the slat lower surface. Furthermore,Fig. 14(c) shows a strong recirculation zone in the flap-cove and the flow acceleration through the vortex and the flap. The DG solver coupled with the modified BL turbulence model displays no adverse robustness implications as a result of these smooth high gradient phenomena.
 |
| Fig. 13 Surface pressure coefficients and skin friction coefficients obtained by 3rd-order DG scheme. |
 |
| Fig. 14 Mach number contours and streamlines for turbulent flow past 30P30N airfoil problem. |
In this article,primary efforts are devoted to the development of efficient and robust DG solver based on structured/unstructured hybrid anisotropic meshes. Techniques are adopted to alleviate the memory gain and improve the computational efficiency for the line-implicit algo- rithms. We propose simplified and effective approaches to obtain the eddy viscosity distribution and further extend the BL model on the current hybrid meshes for the DG methods. Both lam- inar and turbulent RANS calculations are presented using the high-order DG discretization. Numerical results demonstrate convergence rates which are weakly dependent on the aspect- ratio of the cell and the order of the discretization. For viscous flows with shocks,HWENO limiters are shown to be effective in maintaining the accuracy and achieving the non-oscillatory property. Furthermore,using the designed turbulence model,we are able to robustly capture various flow structures and give reliable qualitative numerical predictions for turbulent flow past a complex geometry.
| [1] | Cockburn, B. and Shu, C. W. The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems. J. Comput. Phys., 141, 199-224 (1998) |
| [2] | Cockburn, B., Karniadakis, G., and Shu, C. W. The development of discontinuous Galerkin method. Discontinuous Galerkin Methods, Lecture Notes in Computational Science and Engineering, Springer, New York, 3-50 (2000) |
| [3] | Nastase, C. R. and Mavriplis, D. J. High-order discontinuous Galerkin methods using an hp- multigrid approach. J. Comput. Phys., 213, 330-357 (2006) |
| [4] | Fidkowski, K. J., Oliver, T. A., Lu, J., and Darmofal, D. L. P-multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations. J. Comput. Phys., 207, 92-113 (2005) |
| [5] | Shahbazi, K., Mavriplis, D. J., and Burgess, N. K. Multigrid algorithms for high-order discontinu- ous Galerkin discretizations of the compressible Navier-Stokes equations. J. Comput. Phys., 228, 7917-7940 (2009) |
| [6] | Diosady, L. T. and Darmofal, D. L. Preconditioning methods for discontinuous Galerkin solutions of the Navier-Stokes equations. J. Comput. Phys., 228, 3917-3935 (2009) |
| [7] | Haga, T., Gao, H., and Wang, Z. J. Efficient solution techniques for high-order methods on 3D anisotropic hybrid meshes. The 49th AIAA Aerospace Sciences Meeting Induding the New Horizons Forum and Aerospace Exposition, Curran Associates, Florida (2011) |
| [8] | Burgess, N. K., Nastase, C. R., and Mavriplis, D. J. Efficient solution techniques for discontinuous Galerkin discretizations of the Navier-Stokes equations on hybrid anisotropic meshes. The 48th AIAA Aerospace Sciences Meeting Induding the New Horizons Forum and Aerospace Exposition, American Institute of Aeronautics and Astronautics, Florida (2010) |
| [9] | Mavriplis, D. J. Multigrid strategies for viscous flow solvers on anisotropic unstructured meshes. J. Comput. Phys., 145, 141-165 (1998) |
| [10] | Bassi, F., Crivellini, A., Rebay, S., and Savini, M. Discontinuous Galerkin solution of the Reynolds- averaged Navier-Stokes and k-! turbulence model equations. Computers and Fluids, 34, 507-540 (2005) |
| [11] | Landmann, B., Kessler, M., Wagner, S., and Kramer, E. A parallel, high-order discontinuous Galerkin code for laminar and turbulent flows. Computers and Fluids, 37(4), 427-438 (2008) |
| [12] | Burgess, N. K. and Mavriplis, D. J. Robust computation of turbulent flows using a discontinuous Galerkin method. The 50th AIAA Aerospace Sciences Meeting Induding the New Horizons Forum and Aerospace Exposition, American Institute of Aeronautics and Astronautics, Tennessee (2012) |
| [13] | Li, X., Yang, Y., Hao, H., and Jiao, J. Exploring high-order accurate discontinuous Galerkin method for numerical solution of compressible Reynolds-averaged Navier-Stokes (RANS) equa- tions (in Chinese). Journal of Northwestern Polytechnical University, 30(3), 407-411 (2012) |
| [14] | Mavriplis, D. J. Algebraic turbulence modeling for unstructured and adaptive meshes. AIAA Journal, 29, 2086-2093 (1991) |
| [15] | Jiang, Z., Yan, C., Yu, J., and Yuan, W. Hermite WENO-based limiters for high order discontinuous Galerkin method on unstructured grids. Acta Mechanica Sinica, 28(2), 1-12 (2012) |
| [16] | Bassi, F. and Rebay, S. GMRES discontinuous Galerkin solution of the compressible Navier-Stokes equations. Lecture Note in Computational Science and Engineering, Springer-Verlag, New York, 197-208 (2000) |
| [17] | Zhao, R., Yan, C., and Yu, J. New kind Baldwin-Lomax turbulence model under the limit of entropy (in Chinese). Journal of Beijing University of Aeronautics and Astronautics, 38(2), 175- 180 (2012) |
| [18] | Nielsen, E. J., Anderson, W. K., and Walters, R. W. Application of Newton-Krylov methodol- ogy to a three-dimensional unstructured Euler codes. The 12th Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, California (1995) |
| [19] | Bassi, F. and Rebay, S. A high-order accurate discontinuous finite element method for the nu- merical solution of the compressible Navier-Stokes equations. J. Comput. Phys., 131, 267-279 (1997) |
| [20] | Cook, P. H., McDonald, M. A., and Firmin, M. C. P. Aerofoil RAE2822 Pressure Distributions, Boundary Layer and Wake Measurements, AGARD AR-138, Research and Technology Organi- sation, Neuilly-sur-Seine (1979) |
| [21] | Klausmeyer, S. M. and Lin, J. C. Comparative results from a CFD challenge over a 2D three- element high-lift airfoil. NASA Technical Memorandum, Microfoche, Washington, D. C. (1997) |