A Third Order Accurate Cellwise Relaxation Implicit Discontinuous Galerkin Scheme for Unstructured Hybrid Meshes

Abstract

A third order accurate cellwise relaxation implicit Discontinuous Galerkin (DG) scheme for RANS simulations using unstructured hybrid meshes is presented. A scalar parabolic equation is first examined to clarify what is really important in construction of implicit matrix to keep its diagonal dominance for the third and fourth order cellwise relaxation implicit DG schemes. In addition, discussions are given to approximated construction of implicit matrix for reducing computational cost. Then, the third order accurate cellwise relaxation implicit DG scheme for RANS simulations is successfully developed utilizing the expertise learned in the study of solving the parabolic equation. Superior spatial accuracy of the third order accurate cellwise relaxation implicit DG scheme for RANS simulations, while retaining reasonable convergence properties, is demonstrated for typical aerospace applications.

1. Introduction

Computational fluid dynamics (CFD) is used in many engineering areas, where accurate predictions of flowfield over various vehicles such as aircrafts, automobiles, and ships, are highly required. Unstructured mesh methods are widely used because of its effectiveness in treating complicated vehicle geometries. For unstructured meshes involving various cell geometries, it is customary to employ finite volume method (FVM) for satisfying conservation laws rigorously. However, there is a serious problem that spatial accuracy of FVM is usually lower than that assumed in the formulation even for the second order case. The cause of this lower spatial accuracy comes from poor reconstruction of dependent variables especially when geometries of nearby cells are highly skewed. In addition, slope limiter can degrade the spatial accuracy of the reconstructed dependent variables at local minimum/maximum points. Furthermore, the choice of cell stencils for local reconstruction of dependent variables in two- or three-dimensional space is not unique. These issues make FVM for unstructured meshes not really effective for high order computing for flowfield of engineering interest.

Recently, Discontinuous Galerkin (DG) finite element method attracts attention because of its ability in achieving formal order of spatial accuracy even on unstructured meshes [1]. In this method, basis functions are introduced independently in each cell. The dependent variables are then reconstructed as a sum of these basis functions multiplied by the corresponding degree of freedoms (DOFs). Because DOFs are also independently obtained in each cell as solution unknowns, the resulting reconstructed dependent variables become discontinuous at cell interfaces. The numerical flux is then evaluated by using upwind schemes for convective flux function and by weakly continued viscous flux function. Note that high order spatial accuracy is achieved with information contained solely in each computational cell in DG method. This is in contrast with conventional high order FVMs which require information in nearby cells in reconstruction of dependent variables.

However, DG method has an obvious drawback that computational cost is very high due to increased number of equations for DOFs to be solved in each cell. In addition to implicit integration for achieving faster convergence for steady flow problems, parallel computation is definitely needed for DG method to reduce computing time. It is often pointed out in past works that DG method is quite suitable for parallel computation because of its compactness. Such compactness, however, has not been fully exploited in implicit DG method utilizing, for example, LU-SGS scheme [2]. Yasue et al. proposed the cellwise relaxation implicit DG scheme where implicit time integration was carried out in the same way as the point Jacobi scheme [3]. The second order accurate cellwise relaxation implicit DG scheme was shown to be stable with very large CFL number for solving the RANS equations. It was shown that the cellwise relaxation implicit DG scheme could be easily parallelized using MPI library to achieve more than 120 times speedup when 128 PEs were employed.

Efforts were then made to improve the spatial accuracy of the cellwise relaxation implicit DG scheme for the RANS equations, but a straightforward attempt to improve spatial accuracy was ended at emergence of numerical instability. In the work of Sawada and Yasue, a fourth order cellwise relaxation implicit DG scheme for computing wave propagation was shown to be stable [4]. Therefore, it seems likely that the numerical instability comes from the treatment of viscous terms. In order to clarify the cause of the numerical instability that appeared in the RANS solver, a scalar parabolic equation is instead considered to critically examine the implicit matrix of the cellwise relaxation implicit DG scheme.

In this paper, by examining a scalar parabolic equation, we will discuss rigorous treatment of implicit terms in cellwise relaxation scheme up to a fourth order case and how such treatment can be approximated while retaining numerical stability. We will then apply the rigorous and approximated treatments to third order accurate cellwise relaxation implicit DG method for the RANS equations. The outline of this paper is the following. In Section 2, the cellwise relaxation implicit DG method for a scalar parabolic equation is briefly described. In Section 3, we will show that the rigorous treatment of implicit matrix keeps the diagonal dominance for the third and fourth order accurate schemes, while an easier construction of the implicit matrix actually fails. We will then discuss some approximated construction of implicit matrix for reducing computational cost. In Section 4, the rigorous and approximated constructions of implicit matrix for the third order accurate cellwise relaxation implicit DG scheme are attempted for the RANS equations. The computed results for typical aerospace applications are shown. Finally, conclusions are given in Section 5.

2. Cellwise Relaxation Implicit DG Method for Scalar Parabolic Problem

2.1. Governing Equation

As a typical parabolic equation, we choose a scalar heat conduction equation given by 𝜕𝑢𝜕𝜕𝑡−𝛼²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²+𝜕²𝑢𝜕𝑧²=0, (1) where 𝑢 denotes the scalar variable and 𝛼 the scalar coefficient, assumed as 𝛼=1 for simplicity in the present study. In discretizing the above equation, we introduce the three different coordinate systems, namely, the physical coordinate (𝑥,𝑦,𝑧), the reference coordinate (𝜉,𝜂,𝜁), and the tensor coordinate (𝑟,𝑠,𝑡). Surface integrals of the above equation are evaluated in the reference coordinate using Gaussian quadrature formulae [5, 6], while volume integrals are in the tensor coordinate. The cell geometries considered in this study are tetrahedrons, prisms, pyramids, and hexahedrons. For these cell geometries, the abovementioned three coordinate systems are shown in Figure 1.

Figure 1 

Cell geometries in the physical coordinate (𝑥,𝑦,𝑧), the reference coordinate (𝜉,𝜂,𝜁), and the tensor coordinate (𝑟,𝑠,𝑡) [3].

2.2. Basis Function

In this study, the orthogonal bases functions in the tensor coordinate 𝜙_𝑗(𝑟,𝑠,𝑡) are chosen as described by Warburton [7]. For tetrahedral cells, the orthogonal bases functions are given by 𝜙_𝑗=𝑃_𝑙^0,0(𝑟)1−𝑠2^𝑙𝑃_𝑚^2𝑙+1,0(𝑠)1−𝑡2^𝑙+𝑚𝑃_𝑛^{2𝑙+2𝑚+2,0}(𝑡), (2) for hexahedral cells, 𝜙_𝑗=𝑃_𝑙^0,0(𝑟)𝑃_𝑚^0,0(𝑠)𝑃_𝑛^0,0(𝑡), (3) for prismatic cells, 𝜙_𝑗=𝑃_𝑙^0,0(𝑟)𝑃_𝑚^0,0(𝑠)(1−𝑡)^𝑙𝑃_𝑛^2𝑙+1,0(𝑡), (4) and for pyramidal cells, 𝜙_𝑗=𝑃_𝑙^0,0(𝑟)𝑃_𝑚^0,0(𝑠)(1−𝑡)^𝑙+𝑚𝑃_𝑛^{2𝑙+2𝑚+2,0}(𝑡). (5) In the above equations, 𝑃_𝛾^𝛼,𝛽 denotes the Jacobi polynomials given as 𝑃_𝛾^𝛼,𝛽=(−1)^𝛾2^𝛾()𝛾!1−𝑧^−𝛼(1+𝑧)^−𝛽𝑑^𝛾𝑑𝑧^𝛾(1−𝑧)^𝛼+𝛾(1+𝑧)^𝛽+𝛾. (6) Then subscript 𝑗 in (2) denotes the index of the DOFs defined by a set of subscripts 𝑙, 𝑚, and 𝑛. For the 3D second order case, (𝑙,𝑚,𝑛)=(0,0,0) for 𝑗=1, (𝑙,𝑚,𝑛)=(1,0,0) for 𝑗=2, (𝑙,𝑚,𝑛)=(0,1,0) for 𝑗=3, and (𝑙,𝑚,𝑛)=(0,0,1) for 𝑗=4. Thus, 4 DOFs are needed for the 3D second order case. Similarly, 10 DOFs are needed for the third order case and 20 DOFs for the fourth order case.

2.3. Discontinuous Galerkin Discretization

Let us consider DG discretization of the heat conduction equation. We first transform (1) in the reference coordinate as 𝜕𝑢−𝜕𝑡∇⋅𝜎=0, (7) where the transformed dependent variable 𝑢 and the flux functions 𝜎 are given by 𝑢=𝐽⁻¹𝑢,𝜎_𝜉=𝐽⁻¹𝜉_𝑥𝜕𝑢𝜕𝑥+𝜉_𝑦𝜕𝑢𝜕𝑦+𝜉_𝑧𝜕𝑢,𝜕𝑧𝜎_𝜂=𝐽⁻¹𝜂_𝑥𝜕𝑢𝜕𝑥+𝜂_𝑦𝜕𝑢𝜕𝑦+𝜂_𝑧𝜕𝑢,𝜕𝑧𝜎_𝜁=𝐽⁻¹𝜁_𝑥𝜕𝑢𝜕𝑥+𝜁_𝑦𝜕𝑢𝜕𝑦+𝜁_𝑧𝜕𝑢,𝜕𝑧 (8) and the transformation Jacobian 𝐽⁻¹ is given by 𝐽⁻¹=𝜕𝑥,𝑦,𝑧𝜕.𝜉,𝜂,𝜁 (9) The flux functions can further be rewritten as 𝜎=𝒜∇𝑢, (10) where 𝒜=𝒜_𝜉𝒜_𝜂𝒜_𝜁ℬ_𝜉ℬ_𝜂ℬ_𝜁𝒞_𝜉𝒞_𝜂𝒞_𝜁,𝒜_𝜉=𝐽⁻¹𝜕_𝑥𝜉𝜕_𝑥𝜉+𝜕_𝑦𝜉𝜕_𝑦𝜉+𝜕_𝑧𝜉𝜕_𝑧𝜉,𝒜_𝜂=𝐽⁻¹𝜕_𝑥𝜉𝜕_𝑥𝜂+𝜕_𝑦𝜉𝜕_𝑦𝜂+𝜕_𝑧𝜉𝜕_𝑧𝜂,𝒜_𝜁=𝐽⁻¹𝜕_𝑥𝜉𝜕_𝑥𝜁+𝜕_𝑦𝜉𝜕_𝑦𝜁+𝜕_𝑧𝜉𝜕_𝑧𝜁,ℬ_𝜉=𝐽⁻¹𝜕_𝑥𝜂𝜕_𝑥𝜉+𝜕_𝑦𝜂𝜕_𝑦𝜉+𝜕_𝑧𝜂𝜕_𝑧𝜉,ℬ_𝜂=𝐽⁻¹𝜕_𝑥𝜂𝜕_𝑥𝜂+𝜕_𝑦𝜂𝜕_𝑦𝜂+𝜕_𝑧𝜂𝜕_𝑧𝜂,ℬ_𝜁=𝐽⁻¹𝜕_𝑥𝜂𝜕_𝑥𝜁+𝜕_𝑦𝜂𝜕_𝑦𝜁+𝜕_𝑧𝜂𝜕_𝑧𝜁,𝒞_𝜉=𝐽⁻¹𝜕_𝑥𝜁𝜕_𝑥𝜉+𝜕_𝑦𝜁𝜕_𝑦𝜉+𝜕_𝑧𝜁𝜕_𝑧𝜉,𝒞_𝜂=𝐽⁻¹𝜕_𝑥𝜁𝜕_𝑥𝜂+𝜕_𝑦𝜁𝜕_𝑦𝜂+𝜕_𝑧𝜁𝜕_𝑧𝜂,𝒞_𝜁=𝐽⁻¹𝜕_𝑥𝜁𝜕_𝑥𝜁+𝜕_𝑦𝜁𝜕_𝑦𝜁+𝜕_𝑧𝜁𝜕_𝑧𝜁. (11)(12) For later use, we also write the flux function in the form of 𝜎_𝜉=𝒜_𝜉𝜕𝑢+𝜕𝜉𝒜_𝜂𝜕𝑢+𝜕𝜂𝒜_𝜁𝜕𝑢≡𝜕𝜁𝐴⋅∇𝑢,𝜎_𝜂=ℬ_𝜉𝜕𝑢+𝜕𝜉ℬ_𝜂𝜕𝑢+𝜕𝜂ℬ_𝜁𝜕𝑢≡𝜕𝜁𝐵⋅∇𝑢,𝜎_𝜁=𝒞_𝜉𝜕𝑢+𝜕𝜉𝒞_𝜂𝜕𝑢+𝜕𝜂𝒞_𝜁𝜕𝑢≡𝜕𝜁𝐶⋅∇𝑢. (13)

At first, we consider the BR2 formulation [8, 9] for a weak treatment of (10). Assuming that 𝜓_𝜉 represents a test function, we can obtain the following relation: _Ω𝜓_𝜉𝜎^∗_𝜉𝑑Ω=_𝜕Ω𝜓_𝜉𝑢^∗𝐴⋅𝑛𝑑𝜕Ω−_Ω∇⋅𝜓_𝜉𝐴𝑢𝑑Ω, (14) where 𝜕Ω denotes the interface of computational cell Ω in the reference coordinate, 𝜎^∗ the numerical flux, 𝑢^∗ the dependent variable at the interface, and 𝑛 the unit outward normal vector of 𝜕Ω. In the BR2 formulation, 𝑢^∗ is given by 𝑢^∗=𝑢^𝐿+𝑢^𝑅2on𝜕Ω_𝑖,𝑢_𝑏on𝜕Ω_𝑏, (15) where the superscripts 𝐿 and 𝑅 denote the inside value of the cell interface and that at the opposite side of the cell interface, respectively. On the other hand, the subscripts 𝑖 and 𝑏 denote the cell interface inside of the computational domain and that on the domain boundary, respectively. Therefore, the interface value in (15) for interior interface is simply evaluated by taking an average of dependent variables of both sides, while it is given as the boundary condition on the domain boundary. Then, (14) becomes _Ω𝜓_𝜉𝜎^∗_𝜉𝑑Ω=_𝜕Ω𝜓_𝜉𝑢^∗−𝑢^𝐿𝐴⋅𝑛𝑑𝜕Ω+_Ω𝜓_𝜉𝜎_𝜉𝑑Ω. (16) Therefore, the numerical flux function for 𝜉 component is weakly determined as 𝜎^∗_𝜉=𝜎_𝜉+𝛿_𝜉, (17) where 𝛿_𝜉 denotes the global lifting operator in 𝜉 direction in reference coordinate given by _Ω𝜓_𝜉𝛿_𝜉𝑑Ω=_𝜕Ω𝜓_𝜉𝑢^∗−𝑢^𝐿𝐴⋅𝑛𝑑𝜕Ω. (18) Similarly, the numerical flux functions for 𝜂 and 𝜁 component can be written as 𝜎^∗_𝜂=𝜎_𝜂+𝛿_𝜂,𝜎^∗_𝜁=𝜎_𝜁+𝛿_𝜁, (19) where _Ω𝜓_𝜂𝛿_𝜂𝑑Ω=_𝜕Ω𝜓_𝜂𝑢^∗−𝑢^𝐿𝐵⋅𝑛𝑑𝜕Ω,_Ω𝜓_𝜁𝛿_𝜁𝑑Ω=_𝜕Ω𝜓_𝜁𝑢^∗−𝑢^𝐿𝐶⋅𝑛𝑑𝜕Ω. (20)

Next, we consider the weak formulation of (7). Multiplying (7) by the basis function 𝜙_𝑖 as the test function and integrating over the computational cell, we have _Ω𝜙_𝑖𝜕𝑢𝑑𝜕𝑡Ω−_𝜕Ω𝜙_𝑖𝜎^∗⋅𝑛𝑑𝜕Ω+_Ω∇𝜙_𝑖⋅𝜎^∗𝑑Ω=0, (21) which then yields _Ω𝜙_𝑖𝜕𝑢𝑑𝜕𝑡Ω−_𝜕Ω𝜙_𝑖{𝜎⋅𝑛}𝑑𝜕Ω+_Ω∇𝜙_𝑖⋅𝜎𝑑Ω−𝜂_𝑒𝜕Ω𝜙_𝑖𝛿^𝑓⋅𝑛𝑑𝜕Ω+_Ω∇𝜙_𝑖⋅𝛿𝑑Ω=0, (22) where {⋅} denotes the averaged value of those inside and outside of the cell interface and 𝜂_𝑒 the stabilizing coefficient [9]. The surface integral of the global lifting term in (21) is replaced by the projected local lifting term in (22). It is assumed that the projected local lifting term can be decomposed by the local orthogonal basis as 𝛿^𝑓⋅𝑛=𝛿^𝑓_𝑛=_𝑗𝜙_𝑗𝛿^𝑓_𝑛,𝑗. (23) Therefore, by multiplying the local basis function and integrating over the computational cell, we have _Ω𝜙_𝑘𝛿^𝑓_𝑛𝑑Ω=_Ω𝜙_𝑘𝑗𝜙_𝑗𝛿^𝑓_𝑛,𝑗𝑑Ω=𝐼_𝑘𝑘𝛿^𝑓_𝑛,𝑘, (24) where we utilize the orthogonal property of basis function denoted by 𝐼_𝑖𝑗=_Ω𝜙_𝑖𝜙_𝑗𝑑Ω,0.𝑖=𝑗𝑖≠𝑗 (25) The coefficient in (23) is then obtained as 𝛿^𝑓_𝑛,𝑘=1𝐼_{𝑘𝑘𝜕Ω𝑓}𝜙_𝑘𝑛⋅𝒜𝑛𝑢^∗−𝑢^𝐿𝑑𝜕Ω, (26) which gives _𝜕Ω𝜙_𝑖𝛿^𝑓⋅𝑛𝑑𝜕Ω=_𝜕Ω𝜙_𝑖𝑗𝜙_𝑗𝛿^𝑓_𝑛,𝑗𝑑𝜕Ω. (27)

We note here that taking an average value in (27) is not straightforward because the local lifting term is defined in the reference space where the relation for coordinate transformation and thus the metric terms are not the same for those inside and outside of the interface. Therefore, the average is taken in the physical space and then transformed back to the reference space as shown below. From (22), we can write 𝜎^∗⋅𝑛=𝜎+𝛿^𝑓⋅𝑛=𝑃_𝑥𝜎_𝑥+𝛿^𝑓_𝑥+𝑃_𝑦𝜎_𝑦+𝛿^𝑓_𝑦+𝑃_𝑧𝜎_𝑧+𝛿^𝑓_𝑧, (28) where (𝑃_𝑥,𝑃_𝑦,𝑃_𝑧) can be found by substituting the following relations into (28) as 𝜎_𝜉+𝛿^𝑓_𝜉=𝐽⁻¹𝜉_𝑥𝜎_𝑥+𝛿^𝑓_𝑥+𝜉_𝑦𝜎_𝑦+𝛿^𝑓_𝑦+𝜉_𝑧𝜎_𝑧+𝛿^𝑓_𝑧,𝜎_𝜂+𝛿^𝑓_𝜂=𝐽⁻¹𝜂_𝑥𝜎_𝑥+𝛿^𝑓_𝑥+𝜂_𝑦𝜎_𝑦+𝛿^𝑓_𝑦+𝜂_𝑧𝜎_𝑧+𝛿^𝑓_𝑧,𝜎_𝜁+𝛿^𝑓_𝜁=𝐽⁻¹𝜁_𝑥𝜎_𝑥+𝛿^𝑓_𝑥+𝜁_𝑦𝜎_𝑦+𝛿^𝑓_𝑦+𝜁_𝑧𝜎_𝑧+𝛿^𝑓_𝑧. (29) Then, (28) can be further rewritten as 𝜎^∗⋅𝑛=𝑃²_𝑥+𝑃²_𝑦+𝑃²_𝑧×𝑚_𝑥𝜎_𝑥+𝛿^𝑓_𝑥+𝑚_𝑦𝜎_𝑦+𝛿^𝑓_𝑦+𝑚_𝑧𝜎_𝑧+𝛿^𝑓_𝑧, (30) where 𝑚_𝑥,𝑚_𝑦,𝑚_𝑧=1𝑃²_𝑥+𝑃²_𝑦+𝑃²_𝑧𝑃_𝑥,𝑃_𝑦,𝑃_𝑧. (31) Note that (𝜎_𝑥+𝛿^𝑓_𝑥,𝜎_𝑦+𝛿^𝑓_𝑦,𝜎_𝑧+𝛿^𝑓_𝑧) in (30) can be found either from (29) or by evaluating them directly in the physical space.

Finally, we write the scalar variable as a sum of basis functions multiplied by the related DOFs as 𝑢=_𝑗𝜙_𝑗𝑢_𝑗. (32) Therefore, (22) can be read as _𝑗𝑑𝑢_𝑗𝐼𝑑𝑡^𝑖𝑗−_𝜕Ω𝜙_𝑖{𝜎⋅𝑛}𝑑𝜕Ω+_Ω∇𝜙_𝑖⋅𝜎𝑑Ω−𝜂_𝑒𝜕Ω𝜙_𝑖𝛿^𝑓⋅𝑛𝑑𝜕Ω+_Ω∇𝜙_𝑖⋅𝛿𝑑Ω=0. (33)

2.4. Cellwise Relaxation Implicit Scheme

We first write the time linearization of the flux functions 𝜎 in the reference space as 𝜎^𝑛+1=𝜎^𝑛+𝜕𝜎𝜕𝑡Δ𝑡=𝜎^𝑛+𝜕𝜎𝜕Δ∇𝑢∇𝑢=𝜎^𝑛+𝒜_𝑗∇𝜙_𝑗Δ𝑢^𝑛_𝑗, (34) where superscript 𝑛 denotes time step. Therefore, the implicit treatment of the second term in (33) is given by _𝜕Ω𝜙_𝑖{𝜎⋅𝑛}^𝑛+1𝑑𝜕Ω≃_𝜕Ω𝜙_𝑖{𝜎⋅𝑛}^𝑛𝑑𝜕Ω+_𝑗𝜕Ω𝜙_𝑖𝒜∇𝜙_𝑗⋅𝑛Δ𝑢^𝑛_𝑗𝑑𝜕Ω. (35) Similarly, the third term in (33) becomes _Ω∇𝜙_𝑖⋅𝜎^𝑛+1𝑑Ω≃_Ω∇𝜙_𝑖⋅𝜎^𝑛𝑑Ω+_𝑗Ω∇𝜙_𝑖⋅𝒜∇𝜙_𝑗Δ𝑢^𝑛_𝑗𝑑Ω. (36) Next, the fourth term corresponding to the surface integral of lifting terms in (33) is approximated as _𝜕Ω𝜙_𝑖𝛿^𝑓⋅𝑛^𝑛+1𝑑𝜕Ω≃_𝜕Ω𝜙_𝑖𝛿^𝑓⋅𝑛^𝑛𝑑𝜕Ω+_𝜕Ω𝜙_𝑖𝑘𝜙_𝑘𝜕𝛿^𝑓_𝑛,𝑘𝜕𝑡Δ𝑡𝑑𝜕Ω=_𝜕Ω𝜙_𝑖𝛿^𝑓⋅𝑛^𝑛𝑑𝜕Ω+_𝑗𝜕Ω𝜙_𝑖𝑘𝜙_𝑘𝜕𝛿^𝑓_𝑛,𝑘𝜕𝑢𝜙_𝑗𝐽⁻¹Δ𝑢^𝑛_𝑗𝑑𝜕Ω, (37) where 𝜕𝛿^𝑓_𝑛,𝑘𝜕𝑢=1𝐼_{𝑘𝑘𝜕Ω𝑓}𝜙_𝑘𝑛⋅𝒜𝑛−12𝐽𝑑𝜕Ω_𝑓. (38) In the above equation, the surface integration is taken only on the local interface to maintain the compactness. On the other hand, the fifth term corresponding to the volume integral of lifting terms in (33) is approximated as _Ω∇𝜙_𝑖⋅𝛿^𝑛+1𝑑Ω≃_Ω∇𝜙_𝑖⋅𝛿^𝑛𝑑Ω+_𝜕Ω∇𝜙_𝑖⋅𝒜𝑛−12𝐽Δ𝑢𝑑𝜕Ω=_Ω∇𝜙_𝑖⋅𝛿^𝑛𝑑Ω+_𝑗𝜕Ω∇𝜙_𝑖⋅𝒜𝑛−12𝜙_𝑗Δ𝑢^𝑛_𝑗𝑑𝜕Ω. (39) Finally, we approximate the first term in (33) as _𝑗𝑑𝑢_𝑗𝐼𝑑𝑡^𝑖𝑗=1Δ𝑡_𝑗𝐼_𝑖𝑗𝐽⁻¹Δ𝑢^𝑛_𝑗. (40)

We now obtain an algebraic equation for temporal change of transformed DOFs, Δ𝑢^𝑛_𝑗, in each computational cell given by _𝑗𝑀_𝑖𝑗Δ𝑢^𝑛_𝑗=𝑅_𝑖, (41) where 𝑀_𝑖𝑗=1𝐼Δ𝑡^𝑖𝑗𝐽⁻¹−_𝜕Ω𝜙_𝑖𝒜∇𝜙_𝑗⋅𝑛𝑑Ω+_Ω∇𝜙_𝑖⋅𝒜∇𝜙_𝑗𝑑Ω−_𝜕Ω𝜙_𝑖𝑘𝜙_𝑘𝜕𝛿^𝑓_𝑛,𝑘𝜕𝑢𝜙_𝑗𝐽⁻¹𝑑𝜕Ω+_𝜕Ω∇𝜙_𝑖⋅𝒜𝑛−12𝜙_𝑗𝑑𝜕Ω,𝑅_𝑖=_𝜕Ω𝜙_𝑖{𝜎⋅𝑛}^𝑛𝑑𝜕Ω−_Ω∇𝜙_𝑖⋅𝜎^𝑛𝑑Ω+𝜂_𝑒𝜕Ω𝜙_𝑖𝛿^𝑓⋅𝑛^𝑛𝑑𝜕Ω−_Ω∇𝜙_𝑖⋅𝛿^𝑛𝑑Ω. (42)(43) The size of matrix 𝑀_𝑖𝑗 depends on the number of DOFs. For example, the size becomes 4×4 for 3D second order case and 10×10 for 3D third order case.

Here, we emphasize that 𝑀_𝑖𝑗 in (41) only needs the information of current cell. Therefore, (41) can be calculated in each cell independently as if it were an explicit scheme. We also note that a conventional implicit scheme such as LU-SGS scheme [2] probably can converge faster than the present cellwise one because the information over the entire computational domain is used for every time step through the implicit matrix inversion. However, the present implicit scheme can possibly have higher efficiency of parallelization than other conventional implicit schemes. We would like to see the consequence of choosing the present cellwise approach for DG scheme.

3. Computed Results for Scalar Parabolic Equation

In this chapter, we first show the computed results for the scalar parabolic equation using the second, third, and fourth order accurate cellwise relaxation implicit DG schemes. Spatial accuracies will be examined by comparing with the exact solution. Convergence properties will also be indicated. Then, we consider approximated construction of 𝑀_𝑖𝑗 in (41).

3.1. Numerical Conditions

We consider a scalar heat conduction equation. At first, a cubic domain is divided by Cartesian meshes having 5, 10, and 20 intervals along each coordinate line. The resulting hexahedral cell is further divided into 6 tetrahedral cells. Therefore, the total number of cells becomes 750, 6,000, and 48,000, respectively. We denote these unstructured meshes by Mesh-1a, -1b, and -1c. Then, these meshes are compressed in one coordinate direction until the cell aspect ratio becomes as large as 1,000. We denote the compressed unstructured meshes by Mesh-2a, -2b, and -2c. Figure 2 shows Mesh-1c and Mesh-2c. The boundary condition assumed on all the boundary surfaces of computational domain is given by the following analytic solution for steady state of a scalar heat conduction equation: 𝑢=(𝑥−0.5)²+𝑦−0.5²−2(𝑧−0.5)²+_𝑖,𝑗𝑒^−𝑖𝑖sin𝑗−𝑗cos𝑗, (44) where (𝑖,𝑗) is cyclically taken as 𝑖,𝑗=𝑥,𝑦,𝑦,𝑧,(𝑧,𝑥). (45) The steady analytic solution for the entire computational domain is also given by (44). The initial condition is to assume 𝑢=0 for all computational cell. The time step is so determined that CFL number becomes 10⁶ if a convection velocity of 1.0 is assumed.

(a) Mesh-1c

(b) Mesh-2c

Figure 2 

Unstructured meshes used for solving the heat conduction equation; (a) cell aspect ratio of 1 and (b) 1,000.

The stabilizing coefficient 𝜂_𝑒 is 1.0 for the case in which mesh aspect ratio is 1 and 1.5 for 1,000. All computations in this section are carried out using Intel(R) Core(TM) i7 CPU 3.47 GHz.

3.2. Computed Results

At first we examine the convergence characteristics and spatial accuracies of present schemes. In Table 1, the number of time steps needed for machine zero convergence and the CPU time needed for advancing one time step are summarized. For all spatial accuracies, the present implicit scheme gives machine zero convergence. However, the convergence becomes slower due to larger stabilizing coefficient 𝜂_𝑒 when a cell aspect ratio of 1,000 is chosen. The computational cost of the fourth order scheme for advancing one time step is about 20 times higher than that for the second order scheme. In addition, because the convergence rate for the fourth order scheme is slower than that of the second order scheme, the computational cost of the fourth order scheme for attaining machine zero convergence is about 55 times higher than that of the second order scheme. The error norms obtained for converged solutions are also shown in Table 1. One can evaluate the actual spatial accuracy from the slope of fitted line for error norms with different cell intervals. The attained actual spatial accuracies are also shown in Table 1. The present second, third, and fourth order DG schemes will attain the formal spatial accuracy not only for a cell aspect ratio of 1 but also for 1,000.

Table 1 

Computed results for the parabolic equation using 2nd, 3rd, and 4th order schemes.

Next, let us consider the possible reason why the second order cellwsie relaxation implicit DG scheme of Yasue et al. was numerically stable for solving the RANS equations while the third order scheme became unstable. In solving the RANS equations, Yasue et al. ignored the linearization of the derivative terms for conservative variables. Those second and third terms in (42) were not accounted for in their scheme. We denote the matrix where the second and third terms are ignored by reduction-matrix 𝑀_𝑖𝑗 while the matrix in (42) by the original-matrix 𝑀_𝑖𝑗. In Figures 3 and 4, comparisons of the convergence histories utilizing 𝑀_𝑖𝑗 and 𝑀_𝑖𝑗 are shown. In all calculations, Mesh-2b is employed. For the second order case, the converging rate with the reduction-matrix 𝑀_𝑖𝑗 is almost the same as that with the original-matrix 𝑀_𝑖𝑗, while the computational cost using 𝑀_𝑖𝑗 is reduced. However, for the third and fourth order cases, the computing with the reduction-matrix 𝑀_𝑖𝑗 becomes unstable.

(a) 2nd order

(b) 3rd order

(c) 4th order

Figure 3 

Convergence histories in terms of the number of iterations compared for original-matrix, reduction-matrix, and diagonal-matrix, for (a) 2nd order, (b) 3rd order, and (c) 4th order.

(a) 2nd order

(b) 3rd order

(c) 4th order

Figure 4 

Convergence histories in terms of CPU time compared for original-matrix, reduction-matrix, and diagonal-matrix, for (a) 2nd order, (b) 3rd order, and (c) 4th order.

It seems that the observed instability comes from loss of diagonal dominance of matrix 𝑀_𝑖𝑗. In order to see whether it is the case or not, let us consider for a sake of simplicity a cubic domain divided by Cartesian meshes and then compress only in 𝑧 direction to create a cuboid domain having hexahedral cells. Then, the matrix 𝒜 in (11) has a diagonal form as 𝒜=𝒜_𝜉000ℬ_𝜂000𝒞_𝜁, (46) where 𝒞_𝜁 can be much larger than 𝒜_𝜉 and ℬ_𝜂. When the second and third terms in (42) are ignored, the remaining terms containing 𝒞_𝜁 become −_𝜕Ω𝜙_𝑖𝑘𝜙_𝑘𝜕Ω𝜙_𝑘𝒞_𝜁𝑛²_𝜁−12𝐽𝑑𝜕Ω𝜙_𝑗𝐽⁻¹𝑑𝜕Ω+_𝜕Ω𝜕𝜙_𝑖𝜕𝜁𝒞_𝜁𝑛_𝜁−12𝜙_𝑗𝑑𝜕Ω. (47) Because the first term in the above equation contains 𝑛²_𝜁, the surface integration in the reference space becomes nonzero only on the interfaces 𝜁=−1 and 𝜁=1. Furthermore, on these surfaces, 𝜙_𝑖𝜙_𝑗 needs to be even function of 𝜉 and 𝜂. Otherwise, the surface integral on these surfaces vanishes. Similarly, the second term in (47) contains 𝑛_𝜁. Again, the surface integration can have a nonzero value only on the interfaces 𝜁=−1 and 𝜁=1, where (𝜕𝜙_𝑖/𝜕𝜁)𝜙_𝑗 needs to be even function of 𝜉 and 𝜂. For the second order case, the bases functions for a hexahedron are given by 𝜙₁=1,𝜙₂=𝜉,𝜙₃=𝜂,𝜙₄=𝜁, (48) which are the first order polynomials. Because 𝜙_𝑖𝜙_𝑗 cannot be an even function of 𝜉 and 𝜂, except for 𝜙₁𝜙₄, the first term in (47) vanishes for 𝑖≠𝑗. For 𝜙₁𝜙₄, the sum of surface integration on 𝜁=−1 and 𝜁=1 vanishes. Similarly, the second term in (47) vanishes for 𝑖≠𝑗 except for (𝜕𝜙₄/𝜕𝜁)𝜙₁. For (𝜕𝜙₄/𝜕𝜁)𝜙₁, the sum of surface integration on 𝜁=−1 and 𝜁=1 again vanishes. Therefore, the diagonal dominance is maintained for the second order case. For the third order case, the bases functions are given in addition to (48) as 𝜙₅1=−21−3𝜉²,𝜙₆=𝜉𝜂,𝜙₇1=−21−3𝜂²,𝜙₈𝜙=𝜂𝜁,₉1=−21−3𝜁²,𝜙₁₀=𝜉𝜁, (49) which are the second order polynomials. For these bases functions, (47) can have a nonzero value when 𝜙_𝑖𝜙_𝑗 or (𝜕𝜙_𝑖/𝜕𝜁)𝜙_𝑗 become even functions of 𝜉 and 𝜂. Although the diagonal term is also reinforced with these bases functions, we numerically observed that the resulting reduction-matrix 𝑀_𝑖𝑗 has no diagonal dominance. The situation is the same for the fourth order case. Although Yasue et al. employed the reduction-matrix 𝑀_𝑖𝑗 in their cellwise relaxation implicit DG scheme, the diagonal dominance was kept because the spatial accuracy of their code was second order.

An easier approach to maintain the diagonal dominance of the reduction-matrix 𝑀_𝑖𝑗 is to remove all the off-diagonal terms. We denote this matrix by the diagonal-matrix 𝑀_𝑖𝑗. The convergence histories for computing with 𝑀_𝑖𝑗 are also shown in Figures 3 and 4. For the third order case, the computing with 𝑀_𝑖𝑗 is found to be stable, but convergence is slower than that with the original-matrix 𝑀_𝑖𝑗. On the other hand, for the fourth order case, the computing with 𝑀_𝑖𝑗 is unstable. Therefore, we should take into account the terms arising from the linearization of the derivative terms for conservative variables in the third and fourth order cellwise relaxation implicit DG scheme.

Now that the linearization of the derivative terms for conservative variables should be included in implicit matrix for diagonal dominance; the computing cost to obtain the original-matrix 𝑀_𝑖𝑗 in each cell is very high. Therefore, it is favorable to examine a possibility of simplifying the original-matrix 𝑀_𝑖𝑗 in which some components of matrix 𝑀_𝑖𝑗 are removed for reducing computational cost while retaining the diagonal dominance. For the second order case, we employ the simplified form of 𝑀_𝑖𝑗 as 𝑀₁₁𝑀₁₂𝑀₁₃𝑀₁₄𝑀₂₁𝑀₂₂𝑀00₃₁0𝑀₃₃0𝑀₄₁00𝑀₄₄Δ𝑄₁Δ𝑄₂Δ𝑄₃Δ𝑄₄=𝑅₁𝑅₂𝑅₃𝑅₄, (50) while, for the third order case,𝑀₁₁𝑀₁₂𝑀₁₃𝑀₁₄𝑀₁₅𝑀₁₆𝑀₁₇𝑀₁₈𝑀₁₉𝑀₁₁₀𝑀₂₁𝑀₂₂00𝑀₂₅𝑀₂₆000𝑀₂₁₀𝑀₃₁0𝑀₃₃00𝑀₃₆𝑀₃₇𝑀₃₈𝑀00₄₁00𝑀₄₄000𝑀₄₈𝑀₄₉𝑀₄₁₀𝑀₅₁𝑀₅₂00𝑀₅₅𝑀00000₆₁𝑀₆₂𝑀₆₃00𝑀₆₆𝑀0000₇₁0𝑀₇₃000𝑀₇₇𝑀000₈₁0𝑀₈₃𝑀₈₄000𝑀₈₈𝑀00₉₁00𝑀₉₄0000𝑀₉₉0𝑀₁₀₁𝑀₁₀₂0𝑀₁₀₄00000𝑀₁₀₁₀Δ𝑄₁Δ𝑄₂Δ𝑄₃Δ𝑄₄Δ𝑄₅Δ𝑄₆Δ𝑄₇Δ𝑄₈Δ𝑄₉Δ𝑄₁₀=𝑅₁𝑅₂𝑅₃𝑅₄𝑅₅𝑅₆𝑅₇𝑅₈𝑅₉𝑅₁₀. (51)In (2), a set of (𝑙,𝑚,𝑛) determines the basis function. The simplified matrices in the above equations are derived by considering the relation of 𝑙, 𝑚, and 𝑛. For example, assume that (𝑙,𝑚,𝑛)=(1,1,0). Because the first order polynomials of 𝜉 and 𝜂 in the corresponding basis function can be derived from the combination of (𝑙,𝑚,𝑛)=(0,0,0), (1,0,0), and (0,1,0), we should retain these components of 𝑀₆₁, 𝑀₆₂, 𝑀₆₃, and 𝑀₆₆, while we all ignore 𝑀₆₄, 𝑀₆₅, 𝑀₆₇, 𝑀₆₈, 𝑀₆₉, and 𝑀₆₁₀. Similarly, a simplified form of the original-matrix 𝑀_𝑖𝑗 for the fourth order case can be obtained. It is found that the present approach of simplifying the original-matrix 𝑀_𝑖𝑗 is quite effective in reducing computational cost particularly for the higher order cases. In addition, we apply the successive overrelaxation (SOR) method in solving (41) because the size of matrices for higher order cases becomes very large. In the present SOR method, we first prepare the initial guess for the second order components as Δ𝑄₁=𝜔𝑀₁₁𝑅₁,Δ𝑄₂=𝜔𝑀₂₂𝑅₂−𝑀₂₁Δ𝑄₁,Δ𝑄₃=𝜔𝑀₃₃𝑅₃−𝑀₃₁Δ𝑄₁,Δ𝑄₄=𝜔𝑀₄₄𝑅₄−𝑀₄₁Δ𝑄₁. (52) Using the initial guess for the second order components, the forward sweep to obtain the third order components is carried out as Δ𝑄₅=(1−𝜔)Δ𝑄₅+𝜔𝑀₅₅𝑅₅−𝑀₅₁Δ𝑄₁+𝑀₅₂Δ𝑄₂,Δ𝑄₆=(1−𝜔)Δ𝑄₆+𝜔𝑀₆₆𝑅₆−𝑀₆₁Δ𝑄₁+𝑀₆₂Δ𝑄₂+𝑀₆₃Δ𝑄₃,Δ𝑄₇=(1−𝜔)Δ𝑄₇+𝜔𝑀₇₇𝑅₇−𝑀₇₁Δ𝑄₁+𝑀₇₃Δ𝑄₃,Δ𝑄₈=(1−𝜔)Δ𝑄₈+𝜔𝑀₈₈𝑅₈−𝑀₈₁Δ𝑄₁+𝑀₈₃Δ𝑄₃+𝑀₈₄Δ𝑄₄,Δ𝑄₉=(1−𝜔)Δ𝑄₉+𝜔𝑀₉₉𝑅₉−𝑀₉₁Δ𝑄₁+𝑀₉₄Δ𝑄₄,Δ𝑄₁₀=(1−𝜔)Δ𝑄₁₀+𝜔𝑀₁₀₁₀×𝑅₁₀−𝑀₁₀₁Δ𝑄₁+𝑀₁₀₂Δ𝑄₂+𝑀₁₀₄Δ𝑄₄. (53) Then the backward sweep for improving the second order components is followed as Δ𝑄₄=(1−𝜔)Δ𝑄₄+𝜔𝑀₄₄𝑅₄−𝑀₄₁Δ𝑄₁+𝑀₄₈Δ𝑄₈+𝑀₄₉Δ𝑄₉+𝑀₄₁₀Δ𝑄₁₀,Δ𝑄₃=(1−𝜔)Δ𝑄₃+𝜔𝑀₃₃𝑅₃−𝑀₃₁Δ𝑄₁+𝑀₃₆Δ𝑄₆+𝑀₃₇Δ𝑄₇+𝑀₃₈Δ𝑄₈,Δ𝑄₂=(1−𝜔)Δ𝑄₂+𝜔𝑀₂₂𝑅₂−𝑀₂₁Δ𝑄₁+𝑀₂₅Δ𝑄₅+𝑀₂₆Δ𝑄₆+𝑀₂₁₀Δ𝑄₁₀,Δ𝑄₁=(1−𝜔)Δ𝑄₁+𝜔𝑀₁₁𝑅₁−𝑀₁₂Δ𝑄₂+𝑀₁₃Δ𝑄₃+𝑀₁₄Δ𝑄₄+𝑀₁₅Δ𝑄₅+𝑀₁₆Δ𝑄₆+𝑀₁₇Δ𝑄₇+𝑀₁₈Δ𝑄₈+𝑀₁₉Δ𝑄₉+𝑀₁₁₀Δ𝑄₁₀. (54) In above equations, the relaxation coefficient is chosen as 𝜔=1.3 which can give fast convergence while retaining numerical stability. Solutions are iteratively improved by repeating (53) and (54) for 3 times in the present study. In Figures 5 and 6, the convergence histories for computing with the original-matrix and the matrix-simplification and that with SOR method are shown. One can find that the matrix-simplification is quite effective in reducing computational cost without affecting converging rate. In addition, the computational cost is further reduced when SOR method is used. As a result, we can successfully reduce 5%, 22%, and 39% of the total CPU times for the second, third, and fourth order cases, respectively, when we employ the matrix-simplification with the SOR method. We emphasize that the matrix-simplification with SOR method is more effective in reducing computational cost as the spatial accuracy becomes higher.

(a) 2nd order

(b) 3rd order

(c) 4th order

Figure 5 

Convergence histories in terms of the number of iterations compared for the original-matrix, matrix-simplification, and matrix-simplification with SOR method, for (a) 2nd order, (b) 3rd order, and (c) 4th order.

(a) 2nd order

(b) 3rd order

(c) 4th order

Figure 6 

Convergence histories in terms of CPU time compared for the original-matrix, matrix-simplification, and matrix-simplification with SOR method, for (a) 2nd order, (b) 3rd order, and (c) 4th order.

4. Computed Results for Navier-Stokes Equations or RANS Equations

In this section, we show the computed results for the Navier-Stokes equations using the second and third order accurate cellwise relaxation implicit DG schemes. The second order cellwise relaxation implicit DG scheme for RANS equations including the treatment of convection terms is thoroughly described in [3]. In this study, the second and third terms in (42) are included in the coefficient matrices, and those bases functions for high order terms, their associated DOFs, and sufficient number of Gaussian points are simply introduced to attain the third order spatial accuracy. In the following, we at first show the computed results for the laminar boundary layer flow and compare them with the Blasius solution. Then, we show the computed results for the turbulent boundary layer flow. Finally, the computed results for the vortical flowfield over a double-delta wing are shown.

4.1. Laminar Boundary Layer

The computational domain for the laminar boundary layer flow is set as −2.0≤𝑥≤1.0,0.0≤𝑦≤5.0,0.0≤𝑧≤0.05. (55) A uniform inflow is set at the inlet boundary (𝑥=−2.0) and the outflow condition with a prescribed pressure is set at the outlet boundary (𝑥=1.0) and at the top boundary (𝑦=2.0). The symmetric boundary is assumed not only at the side boundaries (𝑧=0.0,𝑧=0.05) but also at the bottom boundary of upstream region (−2.0≤𝑥≤0.0). At the bottom boundary of downstream region (0.0≤𝑥≤1.0), the nonslip boundary condition is set. A typical example of computational mesh is shown in Figure 7. The present computational mesh consists of 1,688 tetrahedral cells and 1,000 prismatic cells. The minimum gird spacing Δ𝑦_min for the prismatic layer on the wall is set to be 3.0×10⁻⁴ and 10 layers with a growth rate of 1.3 are given in the boundary layer region. The free-stream Mach number is chosen to be 0.2 and the Reynolds number is 1.69×10⁵ for a unit length. The computation is carried out by the cellwise relaxation implicit scheme in which the matrix-simplification with SOR method is employed. In the present case, we choose the relaxation coefficient in SOR method to be 0.9 and the number of iterations to be 10 in order to avoid emergence of numerical instability. A CFL number of 10⁶ is assumed in the implicit time integration. The stabilizing coefficient is assumed as 𝜂_𝑒=1.0. All computations in this subsection are carried out using Intel(R) Core(TM) i7 CPU 3.47 GHz.

Figure 7 

Unstructured hybrid mesh for computation of the laminar boundary layer over a flat plate. Those meshes in the prismatic layer are not shown.

In Figures 8 and 9, one can find that the computed velocity profiles obtained by the second and third order schemes agree well with Blasius solution at several streamwise locations. On the other hand, in Figures 10 and 11, the velocity profiles at 𝑥=0.9 are shown to duplicate each other for both the second and third order cases suggesting that the use of the matrix-simplification with/without SOR method does not alter the converged solutions. The convergence histories determined by density residuals are shown in Figure 12. It is evident that the use of the matrix-simplification with/without SOR method does not influence convergence property. In Table 2, the CPU time needed for advancing one time step is summarized. One can find that the use of the matrix-simplification reduces the needed CPU time by 6.75% and 35.58% for the second and third order schemes, respectively. When the matrix-simplification with SOR method is employed, the CPU time is further reduced by 9.65% and 42.51% for the second and third order schemes, respectively.

Table 2 

CPU time needed for advancing one time step for computing the laminar boundary layer flow.

(a) 2nd order

(b) 3rd order

Figure 8 

Velocity profiles in the laminar boundary layer over a flat plate in 𝑥 direction; (a) 2nd order and (b) 3rd order.

(a) 2nd order

(b) 3rd order

Figure 9 

Velocity profiles in the laminar boundary layer over a flat plate in 𝑦 direction; (a) 2nd order and (b) 3rd order.

(a) 2nd order

(b) 3rd order

Figure 10 

Influence of the matrix-simplification with/without SOR method for computing velocity profiles of the laminar boundary layer in 𝑥 direction at 𝑥=0.9; (a) 2nd order and (b) 3rd order.

(a) 2nd order

(b) 3rd order

Figure 11 

Influence of the matrix-simplification with/without SOR method for computing velocity profiles of the laminar boundary layer in 𝑦 direction at 𝑥=0.9; (a) 2nd order and (b) 3rd order.

(a) 2nd order

(b) 3rd order

Figure 12 

Convergence histories for computing the laminar boundary layer; (a) 2nd order and (b) 3rd order.

4.2. Turbulent Boundary Layer

Next, we consider the turbulent boundary layer flow over a flat plate. An example of the computational mesh is shown in Figure 13. The computational domain is set as −2.0≤𝑥≤1.0,0.0≤𝑦≤2.0,0≤𝑧≤0.05. (56) The boundary conditions are the same for those used in the calculation of the laminar boundary layer flow. The computational mesh consists of 1,496 tetrahedral cells and 2,100 prismatic cells. The minimum grid spacing Δ𝑦_min is set to be 4.0×10⁻⁶. As many as 15 prismatic layers are employed with a growth rate of 1.7 in the boundary layer region. The free-stream Mach number is 0.5, and the Reynolds number for a unit length is 1.05×10⁷. The cellwise relaxation implicit DG scheme is applied for time integration. A CFL number of 10⁶ is assumed. In order to reduce the computational cost, the matrix-simplification with SOR method is employed. The relaxation coefficient and the number of iterations are set to be 0.7 and 15 for this case. The stabilizing coefficient is assumed as 𝜂_𝑒=1.0. All computations are carried out again using Intel(R) Core(TM) i7 CPU 3.47 GHz.

Figure 13 

Unstructured hybrid mesh for computation of the turbulent boundary layer over a flat plate. Those meshes in the prismatic layer are not drawn.

The computed velocity profiles obtained by the second and third order schemes with asymptotic solutions of 𝑢⁺=𝑦⁺ and 𝑢⁺=5.75log₁₀𝑦⁺+5.5 are shown in Figure 14. For both the second and third order schemes, good agreements with the asymptotic solutions are obtained. In Figure 15, the computed results at 𝑥=0.9 obtained by the second and third order schemes are shown. As in the laminar boundary layer case, use of the matrix-simplification with/without SOR method does not influence the converged solutions. Figure 16 shows the velocity profiles of the converged solutions at 𝑥=0.9 given by the second and third order scheme. The velocity profiles obtained by the third order scheme better agree with those asymptotic solutions. The convergence histories determined by the density residuals are shown in Figure 17. Again, the use of the matrix-simplification with/without SOR method does not alter the convergence properties, although the density residual fluctuates in a slightly different way particularly in the third order case. The CPU time needed for advancing one time step for both of the spatial accuracies is summarized in Table 3. The use of the matrix-simplification reduces the needed CPU time by 4.30% and 35.90% for the second and third order schemes, while the matrix-simplification with SOR method by 9.95% and 43.18% for the second and third order scheme, respectively.

Table 3 

CPU time needed for advancing one time step for computing the turbulent boundary layer flow.

(a) 2nd order

(b) 3rd order

Figure 14 

Velocity profiles in the turbulent boundary layer over a flat plate; (a) 2nd order and (b) 3rd order.

(a) 2nd order

(b) 3rd order

Figure 15 

Influence of the matrix-simplification with/without SOR method for computing velocity profiles of the turbulent boundary layer at 𝑥=0.9; (a) 2nd order and (b) 3rd order.

Figure 16 

Velocity profiles at 𝑥=0.9 obtained by the second and third order schemes.

(a) 2nd order

(b) 3rd order

Figure 17 

Convergence histories for computing the turbulent boundary layer flow; (a) 2nd order and (b) 3rd order.

4.3. Vortical Flowfield over a Double-Delta Wing

Finally, the computed results for the flowfield over a double-delta wing are shown. The geometry of the double-delta wing is the same as the model used in the experiment conducted by Brennenstuhl and Hummel [10]. The root chord length is 𝑐_𝑖=0.5 [m]. The kink location is assumed to be the middle of the root chord. The swept-back angles are 80 degree in the upstream region of the kink position and 60 degree in the downstream region, respectively. The strake and main wing are flat plates having a thickness of 𝑑_𝑖=3.0×10⁻³ [m], while the leading edge is rounded by a radius of half wing thickness. An example of the computational mesh is shown in Figure 18. It has 305,619 tetrahedral cells, 343,469 prismatic cells, and 4,095 pyramidal cells. The minimum grid spacing is Δ𝑦_min=3.0×10⁻⁶ [m]. As many as 31 prismatic layers are employed with a growth rate of 1.25 in the boundary layer region. The outer boundary is assumed as a half sphere whose radius is 20 times of the root chord length. The free-stream Mach number is 0.3, and the Reynolds number is 1.3×10⁻⁶ defined by 𝑉Re=^∞𝑐_𝑖𝜈, (57) where 𝑉_∞ denotes the free-stream velocity and 𝜈 the kinematic viscosity. The angle of attack is 12.0 degree. The cellwise relaxation implicit DG scheme is employed. A CFL number of 10³ is chosen with the local time stepping for a fast convergence. The stabilizing coefficient is assumed as 𝜂_𝑒=1.0. For the present case, a parallel computation is performed using 16 processors of the Intel(R) Xeon(R) CPU 3.40 GHz. The METIS grid partitioning [11] and the Message Passing Interface (MPI) library are used in parallelization. In the third order computation, the converged solution of the second order scheme is used as the initial flowfield.

Figure 18 

Unstructured hybrid mesh for computation of a vortical flowfield over a double-delta wing.

The total pressure contour plots in several streamwise cross sections obtained by the second and third order accurate schemes are shown in Figures 19(a) and 19(b), respectively. In both results, two longitudinal vortices from the strake and the leading edge of the main wing are clearly captured which are almost merged at the downstream cross section. In addition, a secondary vortex can be found in this cross section. In Figures 20(a) and 20(b), the contour lines in 62.5% chordwise cross section computed by the second and third order accurate schemes, respectively, are compared. One can find that the outboard longitudinal vortex is better captured by the third order scheme. This is well indicated in the pressure coefficient distributions in 62.5% chordwise cross section shown in Figure 21. Apparently, two negative pressure peaks corresponding to these two longitudinal vortices are clearly observed in the third order computation resulting in better agreement with the experimental data.

(a) 2nd order

(b) 3rd order

Figure 19 

Total pressure contours plotted in streamwise cross section; (a) 2nd order and (b) 3rd order.

(a) 2nd order

(b) 3rd order

Figure 20 

Total pressure contours plotted in 62.5% chordwise cross section; (a) 2nd order and (b) 3rd order.

Figure 21 

Pressure coefficient distributions in 62.5% chordwise cross section.

In this case, the original-matrix is employed in the calculation. The CPU time needed for advancing one time step is 4.8 [s] and 32.3 [s] for the second and third order schemes, respectively, resulting in a ratio of 6.73. This value is close to the ratio of 6.62 obtained in the calculation of the turbulent boundary layer flow using the original-matrix without simplification. We expect to have a better computational efficiency when matrix-simplification and SOR method are employed in the simulation of vortical flowfield over the double-delta wing.

5. Conclusions

A third order accurate cellwise relaxation implicit DG scheme for RANS simulations using unstructured hybrid meshes is presented. In order to find a proper construction method of the implicit matrix for the viscous terms, we first consider a scalar parabolic equation. It is shown that those terms corresponding to linearization of the derivative terms for the conservative variables should be retained in the implicit matrix for the third and fourth order accurate cellwise relaxation implicit DG scheme. In addition, the matrix-simplification with SOR method is proposed to substantially reduce the computational cost. Utilizing the expertise learned in the study of solving the parabolic equation, the third order accurate cellwise relaxation implicit DG scheme for RANS equations is successfully developed. In the calculations of laminar and turbulent boundary layer flow and also vortical flowfield over a double-delta wing, not only the marked numerical stability but also the superior spatial accuracy of the present third order cellwise relaxation implicit DG scheme are well demonstrated.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.