- Email: [email protected]
SUPG finite element computation of inviscid supersonic flows with YZβ shock-Capturing
Computers & Fluids 36 (2007) 147–159 www.elsevier.com/locate/compfluid
SUPG finite element computation of inviscid supersonic flows with YZb shock-Capturing Tayfun E. Tezduyar *, Masayoshi Senga Mechanical Engineering, Rice University, MS 321, 6100 Main Street, Houston, TX 77005, USA Received 27 May 2005; accepted 31 July 2005 Available online 29 November 2005
Abstract Stabilization and shock-capturing parameters introduced recently for the Streamline-Upwind/Petrov–Galerkin (SUPG) formulation of compressible flows based on conservation variables are assessed in test computations with inviscid supersonic flows and different types of finite element meshes. The new shock-capturing parameters, categorized as ‘‘YZb Shock-Capturing’’ in this paper, are compared to earlier parameters derived based on the entropy variables. In addition to being much simpler, the new shock-capturing parameters yield better shock quality in the test computations, with more substantial improvements seen for triangular elements. 2005 Elsevier Ltd. All rights reserved.
1. Introduction The Streamline-Upwind/Petrov–Galerkin (SUPG) formulation was introduced, for the advection–diffusion equation and incompressible flows, in an ASME paper [1], with further studies in a journal paper [2]. Within a couple of years, the SUPG formulation for compressible flows was introduced, in the context of conservation variables, in a NASA technical report [3]. A concise version of that was published as an AIAA paper [4], and a more thorough version with additional examples as a journal paper [5]. Following [3–5], the SUPG formulation for compressible flows was recast in entropy variables and supplemented with a shock-capturing (discontinuity-capturing) term [6]. It was shown in [7] that the SUPG formulation introduced in [3–5], when supplemented with a similar shock-capturing term, is very comparable in accuracy to the one that was recast in
*
Corresponding author. E-mail addresses: [email protected] (T.E. Tezduyar), msenga@ rice.edu (M. Senga). 0045-7930/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2005.07.009
entropy variables. The stabilized formulation introduced in [8] for advection–diffusion–reaction equations also included a discontinuity-capturing term, and precluded augmentation of the SUPG effect by the discontinuity-capturing effect when the advection and discontinuity directions coincide. The SUPG formulations include a stabilization parameter that is mostly known as ‘‘s’’. This parameter involves a measure of the local length scale (also known as ‘‘element length’’) and other parameters such as the element Reynolds and Courant numbers. Various ss were proposed starting with those in [1,2] and [3–5], followed by the one introduced in [8], and those proposed in the subsequently reported SUPG-based methods. Defining a separate s for each degree-of-freedom (i.e. for each equation), leading to a matrix form of the s, was proposed in [9], and generalized in [10]. We call the SUPG formulation introduced in [3–5] for compressible flows ‘‘(SUPG)82’’, and the set of ss introduced in conjunction with that ‘‘s82’’. The s used in [7] with (SUPG)82 is a slightly modified version of s82. A shock-capturing parameter, which we call ‘‘d91’’, was included in the shock-capturing term used in [7].
148
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
Calculation of d91 involves evaluation of the Jacobian matrix of the transformation from the entropy variables to the conservation variables. Subsequent minor modifications of s82 accounted for the interaction between the shock-capturing and the (SUPG)82 terms in a fashion similar to how it was done in [8] for advection–diffusion–reaction equations. All these slightly modified versions of s82 have always been used with the same d91, and we categorize them as ‘‘s82-MOD’’. The discontinuity-capturing directional dissipation (DCDD) stabilization was introduced in [11,12] for incompressible flows with sharp gradients. The DCDD takes effect where there is a sharp gradient in the velocity field and introduces dissipation in the direction of that gradient. The way the DCDD is added to the formulation precludes augmentation of the SUPG effect by the DCDD effect when the advection and discontinuity directions coincide. The DCDD involves a second element length scale, which was also introduced in [11,12] and is based on the solution gradient. This new element length scale is used together with the element length scales defined earlier in [8]. Recognizing this second element length as a diffusion length scale, new stabilization parameters for the diffusive limit were introduced for incompressible flows in [12,13]. Partly based on the ideas underlying the new ss for incompressible flows, new ways of calculating the ss for compressible flows were introduced in [13,14]. More significantly, new ways of calculating the shock-capturing parameters for compressible flows were also introduced in [13,14]. The objective was to have shockcapturing parameters that are simpler, and less costly to compute with, than d91. Some versions of these new shock-capturing parameters are based on ideas underlying the DCDD. Other versions, which we categorize as ‘‘YZb Shock-Capturing’’ in this paper, are based on scaled residuals and are defined with options for smoother or sharper shocks. A preliminary and limited set of test computations with inviscid supersonic flows were reported for YZb Shock-Capturing in [15]. In this paper, we carry out a more comprehensive set of test computations with inviscid supersonic flows, and investigate the performance of the new shock-capturing parameters in computations with different element types and mesh orientations.
oU oFi oEi þ R ¼ 0; ot oxi oxi
ð1Þ
where U = (q, qu1, qu2, qu3, qe) is the vector of conservation variables, and Fi and Ei are, respectively, the Euler and viscous flux vectors: 1
0
C B B ui qu1 þ di1 p C C B C B B Fi ¼ B ui qu2 þ di2 p C C; C B B u qu þ d p C i3 A @ i 3
B B B B Ei ¼ B B B B @
0
ui q
ui ðqe þ pÞ
0 T i1 T i2 T i3
1 C C C C C. C C C A
ð2Þ
qi þ T ik uk
Here dij are the components of the identity tensor I, qi are the components of the heat flux vector, and Tij are the components of the Newtonian viscous stress tensor: T ¼ kð$ uÞI þ 2leðuÞ;
ð3Þ
where k and l (=qm) are the viscosity coefficients, m is the kinematic viscosity, and e(u) is the strain-rate tensor: 1 T eðuÞ ¼ ðð$uÞ þ ð$uÞ Þ. 2
ð4Þ
It is assumed that k = 2l/3. The equation of state used here corresponds to the ideal gas assumption. The term R represents all other components that might enter the equations, including the external forces. Eq. (1) can further be written in the following form: oU oU o oU þ Ai Kij R ¼ 0; ð5Þ ot oxi oxi oxj where Ai ¼
oFi ; oU
Kij
oU ¼ Ei . oxj
ð6Þ
Appropriate sets of boundary and initial conditions are assumed to accompany Eq. (5).
3. SUPG formulations 2. Navier–Stokes equations of compressible flows 3.1. Semi-discrete Let X Rnsd be the spatial domain with boundary C, and (0, T) be the time domain. The symbols q, u, p and e will represent the density, velocity, pressure and the total energy, respectively. The Navier–Stokes equations of compressible flows can be written on X and "t 2 (0, T) as
Given Eq. (5), we form some suitably-defined finitedimensional trial solution and test function spaces ShU and VhU . The SUPG formulation of Eq. (5) can then be written as follows: find Uh 2 ShU such that 8Wh 2 VhU :
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
Z
h Z h h h oU oW h oU h oU þ Ai W dX þ Kij dX ot oxi oxi oxj X X Z Z Wh Hh dC Wh Rh dX h
CH
þ
nel Z X e¼1
X
h oW sSUPG e oxk X
h oU oUh o oUh þ Ahi Ahk Khij Rh dX ot oxi oxi oxj Z n h h el X oW oU mSHOC dX ¼ 0. þ oxi oxi Xe e¼1
ðn el Þn X e¼1
Ahk þ
ð7Þ
3.2. Space–time The space–time version of Eq. (7) can be written based on the Deforming-Spatial-Domain/Stabilized Space–Time (DSD/SST) formulation introduced in [16]. The finite element formulation of the governing equations is written over a sequence of N space–time slabs Qn, where Qn is the slice of the space–time domain between the time levels tn and tn+1. At each time step, the integrations involved in the finite element formulation are performed over Qn. The finite element interpolation functions are discontinuous across the space–time slabs. þ We use the notation ðÞ n and ðÞn to denote the values as tn is approached from below and above respectively. Each Qn is decomposed into space–time elements Qen , where e = 1, 2, . . . , (nel)n. The subscript n used with nel is to account for the general case in which the number of space–time elements may change from one space– time slab to another. For each slab Qn, we form some suitably-defined finite-dimensional trial solution and test function spaces ðShU Þn and ðVhU Þn . In the computations reported here, we use first-order polynomials as interpolation functions. The subscript n implies that corresponding to different space–time slabs we might have different discretizations. The DSD/SST formulation of Eq. (5) can then be written as follows: given ðUh Þn , find h h h h U 2 ðSU Þn such that 8W 2 ðVU Þn : h Z Z h h oU oW h oU h þ Ai W dQ þ ot oxi oxi Qn Qn Z Z h oU Wh Hh dP Wh dQ Khij oxj Pn Qn Z þ þ Rh dQ þ ðWh Þn ððUh Þn ðUh Þn Þ dX
Z sSUPG Qen
h oW oxk
oUh oUh o oUh þ Ahi Khij Rh dQ oxi ot oxi oxj
ðn el Þn X e¼1
Here Hh represents the natural boundary conditions associated with Eq. (5), and CH is the part of the boundary where such boundary conditions are specified. The SUPG stabilization and shock capturing parameters are denoted by sSUPG and mSHOC. They were discussed in Section 1 and will further be discussed in Section 4.
X
þ
149
Z mSHOC Qen
h h oW oU dQ ¼ 0. oxi oxi
ð8Þ
Here Pn is the lateral boundary of the space–time slab. The solution to Eq. (8) is obtained sequentially for all space–time slabs Q0, Q1, Q2, . . . , QN1, and the computah tions start with ðUh Þ 0 ¼ U0 , where U0 is the specified initial value of the vector U.
4. Calculation of the stabilization parameters for compressible flows and shock-capturing Various options for calculating the stabilization parameters and defining the shock-capturing terms in the context of the (SUPG)82 formulation were introduced in [13,14]. In this section we describe those options. For this purpose, we first define the acoustic speed as c, and define the unit vector j as j¼
$qh . jj$qh jj
ð9Þ
As the first option in computing sSUGN1 for each component of the test vector-function W, the stabilization parameters sqSUGN1 , suSUGN1 and seSUGN1 (associated with q, qu and qe, respectively) are defined by the following expression: sqSUGN1
¼
suSUGN1
¼
seSUGN1
¼
nen X
!1 h
ju $N a j
.
ð10Þ
a¼1
As the second option, they are defined as sqSUGN1 ¼ suSUGN1 ¼ seSUGN1 ¼
!1 nen X h ðcjj $N a j þ ju $N a jÞ .
ð11Þ
a¼1
In computing sSUGN2, the parameters sqSUGN2 , suSUGN2 and seSUGN2 are defined as follows: sqSUGN2 ¼ suSUGN2 ¼ seSUGN2 ¼
Dt ; 2
ð12Þ
where Dt is the time step. In computing sSUGN3, the parameter suSUGN3 is defined by using the expression suSUGN3 ¼
h2RGN ; 4m
ð13Þ
150
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
where hRGN ¼ 2
nen X
!1 jr $N a j
h
;
r¼
a¼1
$jju jj . jj$jjuh jjjj
ð14Þ
As the first option in defining the shock-capturing term, first the ‘‘shock-capturing viscosity’’ mSHOC is defined: 2
mSHOC ¼ sSHOC ðuint Þ ;
The parameter seSUGN3 is defined as ðheRGN Þ2 ; ð15Þ 4me where me is the ‘‘kinematic viscosity’’ for the energy equation, !1 nen X $hh e e hRGN ¼ 2 ; ð16Þ jr $N a j ; re ¼ jj$hh jj a¼1 seSUGN3 ¼
where sSHOC ¼
ðseSUPG ÞUGN ¼
1
ðseSUGN1 Þr
þ
1 ðseSUGN2 Þr
þ
1 ðseSUGN3 Þ
ð18Þ 1r . r
b hSHOC jj$qh jjhSHOC ; 2ucha qref
ð21Þ
line plots y
ðsqSUPG ÞUGN ,
and h is the temperature. The parameters ðsuSUPG ÞUGN and ðseSUPG ÞUGN are calculated from their components by using the ‘‘r-switch’’: 1r 1 1 ðsqSUPG ÞUGN ¼ þ ; ð17Þ r r ðsqSUGN1 Þ ðsqSUGN2 Þ 1r 1 1 1 þ þ ; ðsuSUPG ÞUGN ¼ r r r ðsuSUGN1 Þ ðsuSUGN2 Þ ðsuSUGN3 Þ
ð20Þ
Shock
M = 2.0
10˚
M = 1.64
29.3º
ð19Þ Typically, r = 2.
x
Fig. 1. Oblique shock. Problem description.
Fig. 2. Oblique shock. Meshes tested: quadrilateral elements (top) and triangular elements arranged in mesh types A (bottom left) and D (bottom right).
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
hSHOC ¼ hJGN ; hJGN ¼ 2
nen X
ð22Þ
!1 jj $N a j
ð23Þ
.
a¼1
Here ucha is a characteristic velocity (such as kuhk or acoustic speed c or a reference velocity uref), and uint is an intrinsic velocity (such as kuhk or c or ucha). Typically, uint = ucha = uref. The reference density qref is defined as bR =2 q qref ¼ qinf sca ; ð24Þ qinf where qinf is the density at the inflow and qsca is a scaling density (such as qinf or the density behind a normal shock corresponding to the inflow Mach number). The parameter bR is set as bR = 1 for smoother shocks and
151
bR = 2 for sharper shocks (in return for tolerating possible overshoots and undershoots). The parameter b is set as b = 1 for smoother shocks and b = 2 for sharper shocks. The compromise between the b = 1 and b = 2 selections is defined as the following averaged expression for sSHOC: 1 sSHOC ¼ ððsSHOC Þb¼1 þ ðsSHOC Þb¼2 Þ. 2
ð25Þ
As another option for calculating the shock-capturing parameter, mSHOC is defined as mSHOC
!b=21 b nsd X 1 oUh 2 hSHOC ¼ jjY Zjj ; Y ox 2 i i¼1 1
ð26Þ
1.6 CYZ12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y 1.6 CYZU12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y
Fig. 3. Oblique shock. Computed with the mesh made of quadrilateral elements. Density along x = 0.9, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
152
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
where Y is a diagonal scaling matrix constructed from the reference values of the components of U: 3 2 ðU 1 Þref 0 0 0 0 6 0 ðU 2 Þref 0 0 0 7 7 6 7 6 7; 0 0 ðU Þ 0 0 Y¼6 3 ref 7 6 7 6 0 0 ðU 4 Þref 0 5 4 0 0
0
0
0
and b = 1 or b = 2. In a variation of the expression given by Eq. (26), mSHOC is defined by the following expression:
mSHOC
nsd X 1 oUh 2 Y ¼ jjY Zjj ox jjY1 Uh jj1b
ð27Þ h
Z¼
h
oU oU þ Ahi ot oxi
ð28Þ
OR Z ¼ Ahi
oUh ; oxi
ð29Þ
i
i¼1
ðU 5 Þref
!b=21
1
hSHOC 2
b ð30Þ
.
The compromise between the b = 1 and b = 2 selections is defined as the following averaged expression for mSHOC: 1 ðmSHOC Þb¼1 þ ðmSHOC Þb¼2 . ð31Þ mSHOC ¼ 2
1.6 CYZ12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y 1.6 CYZU12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y
Fig. 4. Oblique shock. Computed with the mesh made of triangular elements arranged in mesh type A. Density along x = 0.9, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
153
1.6 CYZ12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y 1.6 CYZU12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y
Fig. 5. Oblique shock. Computed with the mesh made of triangular elements arranged in mesh type D. Density along x = 0.9, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
y
M = 2.3
78
29º 2 M = 2.9 line plots
M = 1.942
Shock
3
1 23º
x
Fig. 6. Reflected shock. Problem description.
We also propose versions of mSHOC that take into account the Mach number and shock intensity across a
shock. In doing that, we modify mSHOC given by Eqs. (26) and (30) as follows:
154
mSHOC
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
mSHOC
! b jj$qh jjhSHOC J 1=bM 1þ hM 1i . qref ð32Þ
Here M is the Mach number, ‘‘h i’’ is the Macaulay bracket: 0; x 6 y; hx yi ¼ ð33Þ x y; x > y; and the parameters bJ and bM can each be set to 1 for smoother shocks and 2 for sharper shocks. Based on Eq. (26), a separate mSHOC can be calculated for each component of the test vector-function W: !b=21 b nsd X 1 oUh 2 hSHOC 1 ðmSHOC ÞI ¼ jðY ZÞI j ; Y ox 2 i I i¼1 Similarly, a separate mSHOC for each component of W can be calculated based on Eq. (30): nsd X 1 oUh 2 Y ðmSHOC ÞI ¼ jðY ZÞI j oxi I i¼1 I ¼ 1;2;...nsd þ 2.
!b=21 jðY1 Uh ÞI j1b
b hSHOC ; 2
ð35Þ
Given mSHOC, the shock-capturing term is defined as nel Z X S SHOC ¼ $Wh : ðjSHOC $Uh Þ dX; ð36Þ e¼1
mSHOC
2
mSHOC switchðsSUPG ðj uÞ ;
sSUPG ðjj uj cÞ2 ; mSHOC Þ;
ð37Þ
where the ‘‘switch’’ function is defined as the ‘‘min’’ function or as the ‘‘r-switch’’ used earlier in this section. For viscous flows, the above modification would be made separately with each of sqSUPG , suSUPG and seSUPG , and this would result in mSHOC becoming a diagonal matrix even if the option given by Eq. (34) or Eq. (35) is not exercised.
ð34Þ
I ¼ 1; 2; . . . nsd þ 2.
1
option given by Eq. (34) or Eq. (35) is exercised, then mSHOC becomes an (nsd + 2) · (nsd + 2) diagonal matrix, and the matrix jSHOC becomes augmented from an nsd · nsd matrix to an (nsd · (nsd + 2)) · ((nsd + 2) · nsd) matrix. To preclude compounding, mSHOC can be modified as follows:
Xe
where jSHOC is defined as jSHOC = mSHOC I. As a possible alternative, it is defined as jSHOC = mSHOC jj. If the
5. Test computations The test computations were carried out by using the space–time SUPG formulation described in Section 3.2. In the option denoted by ‘‘CYZ12’’, we use Eq. (26) with Eq. (31), and in the option denoted by ‘‘CYZU12’’, Eq. (30) with Eq. (31). For smoother shocks (see discussion in Section 4, in the option denoted by ‘‘CYZ1’’, we use Eq. (26) with b = 1. This is equivalent to the option denoted by ‘‘CYZU1’’, where we use Eq. (30) with b = 1. In all these options, we use for Z the expression given by Eq. (29), and set jSHOC = mSHOC I. Also, in all these options, as stabiliza-
Fig. 7. Reflected shock. Meshes tested: quadrilateral elements (top) and triangular elements arranged in mesh types D1A2 (middle) and A1D2 (bottom).
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
tion parameters, we use Eq. (11), and in Eqs. (17)–(19) we do not include sSUGN2. The results are compared to those obtained with the s82-MOD and d91 combination. The version of s82-MOD used in this paper for comparison is similar to the one given in [17]: s82MOD ¼ maxð0; sa sd Þ; ð38Þ where h $jjUh jj hBGN d91 ; sa ¼ ; sd ¼ ; u ¼ c þ u cc 2 2ucc jj$jjUh jjjj ðucc Þ
hBGN
For completeness, we also provide here the expression for d91 [7]: oUh d91 ¼ Ahk ox k
, 1
A0 e
2 !12 nsd X onj h oU oxk ox 1 ; k A0 j¼1 e
ð41Þ
e 0 is the where njÕs are the element coordinates, and A Jacobian of the transformation from the entropy variables to the conservation variables. As test problems, we consider steady-state, inviscid cases. The solutions are obtained with time-marching to the steady-state. At each time step, the coupled nonlinear equations involved are solved with the Newton–Raphson method. An iterative technique with
ð39Þ
!1 nen h X $jjU jj ¼2 . jj$jjUh jjjj $N a
155
ð40Þ
a¼1
3 CYZ12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
2.5
3
3.5
4
x 3 CYZU12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2 x
Fig. 8. Reflected shock. Computed with the mesh made of quadrilateral elements. Density along y = 0.25, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
156
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
5.1. Oblique shock
nodal-block-diagonal preconditioner and GMRES update method [18] is employed for solving the linear equation system involved at each Newton–Raphson iteration. We compute two well-known test problems: ‘‘oblique shock’’ and ‘‘reflected shock’’. These were used in many earlier publications, and here we compute each of them with options CYZ12 and CYZU12 and with three different meshes. For each test problem, we use a mesh with quadrilateral elements and two triangular meshes where the elements are arranged in two different ways. We compute 100 time steps with a time-step size of 0.05. The number of Newton–Raphson iterations and the outer and inner GMRES iterations are 3, 1, and 5, respectively.
Fig. 1 shows the problem description. This is a Mach 2 uniform flow over a wedge at an angle of 10 with the horizontal wall. The solution involves an oblique shock at an angle of 29.3 emanating from the leading edge. The computational domain is a square with 0 6 x 6 1 and 0 6 y 6 1. The inflow conditions are given as M = 2.0, q = 1.0, u1 = cos 10, u2 = sin 10, and p = 0.179. This results in an exact solution with the following outflow data: M = 1.64, q = 1.46, u1 = 0.887, u2 = 0.0, and p = 0.305. All essential boundary conditions are imposed at the left and top boundaries, slip condition at the wall, and no boundary conditions at the right boundary. The mesh is essentially
3 CYZ12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
2.5
3
3.5
4
x 3 CYZU12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2 x
Fig. 9. Reflected shock. Computed with the mesh made of triangular elements arranged in mesh type D1A2. Density along y = 0.25, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
uniform and consists of 21 · 21 nodes. Fig. 2 shows the three different meshes tested. Figs. 3–5 show, for three different meshes, the density along x = 0.9, obtained with CYZ12 and CYZU12, compared with the solution obtained with the s82-MOD and d91 combination. In addition to being much simpler, the new shock-capturing parameters yield shocks with better quality. The improvement over the s82-MOD and d91 combination is even more substantial with the triangular meshes.
157
computational domain is a rectangle with 0 6 x 6 4.1 and 0 6 y 6 1. The inflow conditions in R1 are given as M = 2.9, q = 1.0, u1 = 2.9, u2 = 0.0, and p = 0.7143. Specifying these conditions and requiring the incident shock to be at an angle of 29 results in an exact solution with the following flow data: R2: M = 2.378, q = 1.7, u1 = 2.619, u2 = 0.506, and p = 1.528; R3: M = 1.942, q = 2.687, u1 = 2.401, u2 = 0.0, and p = 2.934. All essential boundary conditions are imposed at the left and top boundaries, slip condition at the wall, and no boundary conditions at the right boundary. The mesh is uniform and consists of 61 · 21 nodes. Fig. 7 shows the three different meshes tested. Figs. 8–10 show, for three different meshes, the density along y = 0.25, obtained with CYZ12 and CYZU12, compared with
5.2. Reflected shock Fig. 6 shows the problem description. This problem consists of three flow regions (R1, R2 and R3) separated by an oblique shock and its reflection from the wall. The
3 CYZ12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
2.5
3
3.5
4
x 3 CYZU12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2 x
Fig. 10. Reflected shock. Computed with the mesh made of triangular elements arranged in mesh type A1D2. Density along y = 0.25, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
158
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159 3 CYZ1 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
x
Fig. 11. Reflected shock. Computed with the mesh made of triangular elements arranged in mesh type A1D2. Density along y = 0.25, obtained with CYZ1 (equivalent to CYZU1), compared with the solution obtained with the s82-MOD and d91 combination.
the solution obtained with the s82-MOD and d91 combination. Again, the new, much simpler shock-capturing parameters yield shocks with better quality, and the improvement over the s82-MOD and d91 combination is even more substantial with the triangular meshes. The small undershoot seen in the solutions obtained with CYZ12 and CYZU12 and mesh type A1D2 can be eliminated by using CYZ1, which is equivalent to CYZU1. Fig. 11 shows, for mesh type A1D2, the density along y = 0.25, obtained with CYZ1, compared with the solution obtained with the s82-MOD and d91 combination. While the small overshoot has been eliminated at the expense of losing some of the sharpness of the shocks, they are still sharper compared to the solution obtained with the s82-MOD and d91 combination.
YZb Shock-Capturing and an earlier shock-capturing term that we call ‘‘d91’’. In addition to being much simpler, YZb Shock-Capturing yields substantially better shock quality, and we see even more substantial improvements for triangular elements.
Acknowledgements This work was supported by the US Army Natick Soldier Center (Contract No. DAAD16-03-C-0051) and NASA Johnson Space Center (Grant No. NAG91435).
References 6. Concluding remarks In the context of test computations with inviscid supersonic flows, we provided a comprehensive assessment of the stabilization and shock-capturing parameters that we introduced recently for the Streamline-Upwind/ Petrov–Galerkin (SUPG) formulation of compressible flows based on conservation variables. We focused on performance evaluation of the shock-capturing parameters that we categorized as ‘‘YZb Shock-Capturing’’ in this paper. We used two well-known test problems: ‘‘oblique shock’’ and ‘‘reflected shock’’. For each test problem, we used a mesh with quadrilateral elements and two triangular meshes where the elements are arranged in two different ways. With emphasis on the shock quality, we compared the solutions obtained with
[1] Hughes TJR, Brooks AN. A multi-dimensional upwind scheme with no crosswind diffusion. In: Hughes TJR, editor. Finite element methods for convection dominated flows, AMD-vol. 34. New York: ASME; 1979. p. 19–35. [2] Brooks AN, Hughes TJR. Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput Meth Appl Mech Eng 1982;32:199–259. [3] Tezduyar TE, Hughes TJR. Development of time-accurate finite element techniques for first-order hyperbolic systems with particular emphasis on the compressible Euler equations. NASA Technical Report NASA-CR-204772, NASA, 1982. Available from: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/ 19970023187_1997034954.pdf. [4] Tezduyar TE, Hughes TJR. Finite element formulations for convection dominated flows with particular emphasis on the compressible Euler equations. In: Proceedings of AIAA 21st aerospace sciences meeting, AIAA Paper 83-0125, Reno, Nevada, 1983.
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159 [5] Hughes TJR, Tezduyar TE. Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations. Comput Meth Appl Mech Eng 1984;45:217–84. [6] Hughes TJR, Franca LP, Mallet M. A new finite element formulation for computational fluid dynamics: VI. Convergence analysis of the generalized SUPG formulation for linear timedependent multi-dimensional advective–diffusive systems. Comput Meth Appl Mech Eng 1987;63:97–112. [7] Le Beau GJ, Tezduyar TE. Finite element computation of compressible flows with the SUPG formulation. Advances in finite element analysis in fluid dynamics, FED-vol. 123. New York: ASME; 1991. p. 21–7. [8] Tezduyar TE, Park YJ. Discontinuity capturing finite element formulations for nonlinear convection–diffusion–reaction equations. Comput Meth Appl Mech Eng 1986;59:307–25. [9] Hughes TJR, Mallet M. A new finite element formulation for computational fluid dynamics: III. The generalized streamline operator for multidimensional advective–diffusive systems. Comput Meth Appl Mech Eng 1986;58:305–28. [10] Shakib F, Hughes TJR, Johan Z. A new finite element formulation for computational fluid dynamics: X. The compressible Euler and Navier–Stokes equations. Comput Meth Appl Mech Eng 1991;89:141–219. [11] Tezduyar TE. Adaptive determination of the finite element stabilization parameters. In: Proceedings of the ECCOMAS computational fluid dynamics conference 2001 (CD-ROM), Swansea, Wales, United Kingdom, 2001.
159
[12] Tezduyar TE. Computation of moving boundaries and interfaces and stabilization parameters. Int J Numer Meth Fluids 2003;43:555–75. [13] Tezduyar TE. Finite element methods for fluid dynamics with moving boundaries and interfaces. In: Stein E, De Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics, vol. 3: Fluids. John Wiley & Sons; 2004. Chapter 17. [14] Tezduyar TE. Determination of the stabilization and shockcapturing parameters in SUPG formulation of compressible flows. In: Proceedings of the European congress on computational methods in applied sciences and engineering, ECCOMAS 2004 (CD-ROM), Jyvaskyla, Finland, 2004. [15] Tezduyar TE, Senga M. Stabilization and shock-capturing parameters in SUPG formulation of compressible flows. Comput Meth Appl Mech Eng, in press. [16] Tezduyar TE, Behr M, Liou J. A new strategy for finite element computations involving moving boundaries and interfaces–the deforming-spatial-domain/space–time procedure: I. The concept and the preliminary numerical tests. Comput Meth Appl Mech Eng 1992;94:339–51. [17] Aliabadi SK, Tezduyar TE. Parallel fluid dynamics computations in aerospace applications. Int J Numer Meth Fluids 1995;21: 783–805. [18] Saad Y, Schultz M. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Scient Statist Comput 1986;7:856–69.
SUPG finite element computation of inviscid supersonic flows with YZb shock-Capturing Tayfun E. Tezduyar *, Masayoshi Senga Mechanical Engineering, Rice University, MS 321, 6100 Main Street, Houston, TX 77005, USA Received 27 May 2005; accepted 31 July 2005 Available online 29 November 2005
Abstract Stabilization and shock-capturing parameters introduced recently for the Streamline-Upwind/Petrov–Galerkin (SUPG) formulation of compressible flows based on conservation variables are assessed in test computations with inviscid supersonic flows and different types of finite element meshes. The new shock-capturing parameters, categorized as ‘‘YZb Shock-Capturing’’ in this paper, are compared to earlier parameters derived based on the entropy variables. In addition to being much simpler, the new shock-capturing parameters yield better shock quality in the test computations, with more substantial improvements seen for triangular elements. 2005 Elsevier Ltd. All rights reserved.
1. Introduction The Streamline-Upwind/Petrov–Galerkin (SUPG) formulation was introduced, for the advection–diffusion equation and incompressible flows, in an ASME paper [1], with further studies in a journal paper [2]. Within a couple of years, the SUPG formulation for compressible flows was introduced, in the context of conservation variables, in a NASA technical report [3]. A concise version of that was published as an AIAA paper [4], and a more thorough version with additional examples as a journal paper [5]. Following [3–5], the SUPG formulation for compressible flows was recast in entropy variables and supplemented with a shock-capturing (discontinuity-capturing) term [6]. It was shown in [7] that the SUPG formulation introduced in [3–5], when supplemented with a similar shock-capturing term, is very comparable in accuracy to the one that was recast in
*
Corresponding author. E-mail addresses: [email protected] (T.E. Tezduyar), msenga@ rice.edu (M. Senga). 0045-7930/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2005.07.009
entropy variables. The stabilized formulation introduced in [8] for advection–diffusion–reaction equations also included a discontinuity-capturing term, and precluded augmentation of the SUPG effect by the discontinuity-capturing effect when the advection and discontinuity directions coincide. The SUPG formulations include a stabilization parameter that is mostly known as ‘‘s’’. This parameter involves a measure of the local length scale (also known as ‘‘element length’’) and other parameters such as the element Reynolds and Courant numbers. Various ss were proposed starting with those in [1,2] and [3–5], followed by the one introduced in [8], and those proposed in the subsequently reported SUPG-based methods. Defining a separate s for each degree-of-freedom (i.e. for each equation), leading to a matrix form of the s, was proposed in [9], and generalized in [10]. We call the SUPG formulation introduced in [3–5] for compressible flows ‘‘(SUPG)82’’, and the set of ss introduced in conjunction with that ‘‘s82’’. The s used in [7] with (SUPG)82 is a slightly modified version of s82. A shock-capturing parameter, which we call ‘‘d91’’, was included in the shock-capturing term used in [7].
148
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
Calculation of d91 involves evaluation of the Jacobian matrix of the transformation from the entropy variables to the conservation variables. Subsequent minor modifications of s82 accounted for the interaction between the shock-capturing and the (SUPG)82 terms in a fashion similar to how it was done in [8] for advection–diffusion–reaction equations. All these slightly modified versions of s82 have always been used with the same d91, and we categorize them as ‘‘s82-MOD’’. The discontinuity-capturing directional dissipation (DCDD) stabilization was introduced in [11,12] for incompressible flows with sharp gradients. The DCDD takes effect where there is a sharp gradient in the velocity field and introduces dissipation in the direction of that gradient. The way the DCDD is added to the formulation precludes augmentation of the SUPG effect by the DCDD effect when the advection and discontinuity directions coincide. The DCDD involves a second element length scale, which was also introduced in [11,12] and is based on the solution gradient. This new element length scale is used together with the element length scales defined earlier in [8]. Recognizing this second element length as a diffusion length scale, new stabilization parameters for the diffusive limit were introduced for incompressible flows in [12,13]. Partly based on the ideas underlying the new ss for incompressible flows, new ways of calculating the ss for compressible flows were introduced in [13,14]. More significantly, new ways of calculating the shock-capturing parameters for compressible flows were also introduced in [13,14]. The objective was to have shockcapturing parameters that are simpler, and less costly to compute with, than d91. Some versions of these new shock-capturing parameters are based on ideas underlying the DCDD. Other versions, which we categorize as ‘‘YZb Shock-Capturing’’ in this paper, are based on scaled residuals and are defined with options for smoother or sharper shocks. A preliminary and limited set of test computations with inviscid supersonic flows were reported for YZb Shock-Capturing in [15]. In this paper, we carry out a more comprehensive set of test computations with inviscid supersonic flows, and investigate the performance of the new shock-capturing parameters in computations with different element types and mesh orientations.
oU oFi oEi þ R ¼ 0; ot oxi oxi
ð1Þ
where U = (q, qu1, qu2, qu3, qe) is the vector of conservation variables, and Fi and Ei are, respectively, the Euler and viscous flux vectors: 1
0
C B B ui qu1 þ di1 p C C B C B B Fi ¼ B ui qu2 þ di2 p C C; C B B u qu þ d p C i3 A @ i 3
B B B B Ei ¼ B B B B @
0
ui q
ui ðqe þ pÞ
0 T i1 T i2 T i3
1 C C C C C. C C C A
ð2Þ
qi þ T ik uk
Here dij are the components of the identity tensor I, qi are the components of the heat flux vector, and Tij are the components of the Newtonian viscous stress tensor: T ¼ kð$ uÞI þ 2leðuÞ;
ð3Þ
where k and l (=qm) are the viscosity coefficients, m is the kinematic viscosity, and e(u) is the strain-rate tensor: 1 T eðuÞ ¼ ðð$uÞ þ ð$uÞ Þ. 2
ð4Þ
It is assumed that k = 2l/3. The equation of state used here corresponds to the ideal gas assumption. The term R represents all other components that might enter the equations, including the external forces. Eq. (1) can further be written in the following form: oU oU o oU þ Ai Kij R ¼ 0; ð5Þ ot oxi oxi oxj where Ai ¼
oFi ; oU
Kij
oU ¼ Ei . oxj
ð6Þ
Appropriate sets of boundary and initial conditions are assumed to accompany Eq. (5).
3. SUPG formulations 2. Navier–Stokes equations of compressible flows 3.1. Semi-discrete Let X Rnsd be the spatial domain with boundary C, and (0, T) be the time domain. The symbols q, u, p and e will represent the density, velocity, pressure and the total energy, respectively. The Navier–Stokes equations of compressible flows can be written on X and "t 2 (0, T) as
Given Eq. (5), we form some suitably-defined finitedimensional trial solution and test function spaces ShU and VhU . The SUPG formulation of Eq. (5) can then be written as follows: find Uh 2 ShU such that 8Wh 2 VhU :
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
Z
h Z h h h oU oW h oU h oU þ Ai W dX þ Kij dX ot oxi oxi oxj X X Z Z Wh Hh dC Wh Rh dX h
CH
þ
nel Z X e¼1
X
h oW sSUPG e oxk X
h oU oUh o oUh þ Ahi Ahk Khij Rh dX ot oxi oxi oxj Z n h h el X oW oU mSHOC dX ¼ 0. þ oxi oxi Xe e¼1
ðn el Þn X e¼1
Ahk þ
ð7Þ
3.2. Space–time The space–time version of Eq. (7) can be written based on the Deforming-Spatial-Domain/Stabilized Space–Time (DSD/SST) formulation introduced in [16]. The finite element formulation of the governing equations is written over a sequence of N space–time slabs Qn, where Qn is the slice of the space–time domain between the time levels tn and tn+1. At each time step, the integrations involved in the finite element formulation are performed over Qn. The finite element interpolation functions are discontinuous across the space–time slabs. þ We use the notation ðÞ n and ðÞn to denote the values as tn is approached from below and above respectively. Each Qn is decomposed into space–time elements Qen , where e = 1, 2, . . . , (nel)n. The subscript n used with nel is to account for the general case in which the number of space–time elements may change from one space– time slab to another. For each slab Qn, we form some suitably-defined finite-dimensional trial solution and test function spaces ðShU Þn and ðVhU Þn . In the computations reported here, we use first-order polynomials as interpolation functions. The subscript n implies that corresponding to different space–time slabs we might have different discretizations. The DSD/SST formulation of Eq. (5) can then be written as follows: given ðUh Þn , find h h h h U 2 ðSU Þn such that 8W 2 ðVU Þn : h Z Z h h oU oW h oU h þ Ai W dQ þ ot oxi oxi Qn Qn Z Z h oU Wh Hh dP Wh dQ Khij oxj Pn Qn Z þ þ Rh dQ þ ðWh Þn ððUh Þn ðUh Þn Þ dX
Z sSUPG Qen
h oW oxk
oUh oUh o oUh þ Ahi Khij Rh dQ oxi ot oxi oxj
ðn el Þn X e¼1
Here Hh represents the natural boundary conditions associated with Eq. (5), and CH is the part of the boundary where such boundary conditions are specified. The SUPG stabilization and shock capturing parameters are denoted by sSUPG and mSHOC. They were discussed in Section 1 and will further be discussed in Section 4.
X
þ
149
Z mSHOC Qen
h h oW oU dQ ¼ 0. oxi oxi
ð8Þ
Here Pn is the lateral boundary of the space–time slab. The solution to Eq. (8) is obtained sequentially for all space–time slabs Q0, Q1, Q2, . . . , QN1, and the computah tions start with ðUh Þ 0 ¼ U0 , where U0 is the specified initial value of the vector U.
4. Calculation of the stabilization parameters for compressible flows and shock-capturing Various options for calculating the stabilization parameters and defining the shock-capturing terms in the context of the (SUPG)82 formulation were introduced in [13,14]. In this section we describe those options. For this purpose, we first define the acoustic speed as c, and define the unit vector j as j¼
$qh . jj$qh jj
ð9Þ
As the first option in computing sSUGN1 for each component of the test vector-function W, the stabilization parameters sqSUGN1 , suSUGN1 and seSUGN1 (associated with q, qu and qe, respectively) are defined by the following expression: sqSUGN1
¼
suSUGN1
¼
seSUGN1
¼
nen X
!1 h
ju $N a j
.
ð10Þ
a¼1
As the second option, they are defined as sqSUGN1 ¼ suSUGN1 ¼ seSUGN1 ¼
!1 nen X h ðcjj $N a j þ ju $N a jÞ .
ð11Þ
a¼1
In computing sSUGN2, the parameters sqSUGN2 , suSUGN2 and seSUGN2 are defined as follows: sqSUGN2 ¼ suSUGN2 ¼ seSUGN2 ¼
Dt ; 2
ð12Þ
where Dt is the time step. In computing sSUGN3, the parameter suSUGN3 is defined by using the expression suSUGN3 ¼
h2RGN ; 4m
ð13Þ
150
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
where hRGN ¼ 2
nen X
!1 jr $N a j
h
;
r¼
a¼1
$jju jj . jj$jjuh jjjj
ð14Þ
As the first option in defining the shock-capturing term, first the ‘‘shock-capturing viscosity’’ mSHOC is defined: 2
mSHOC ¼ sSHOC ðuint Þ ;
The parameter seSUGN3 is defined as ðheRGN Þ2 ; ð15Þ 4me where me is the ‘‘kinematic viscosity’’ for the energy equation, !1 nen X $hh e e hRGN ¼ 2 ; ð16Þ jr $N a j ; re ¼ jj$hh jj a¼1 seSUGN3 ¼
where sSHOC ¼
ðseSUPG ÞUGN ¼
1
ðseSUGN1 Þr
þ
1 ðseSUGN2 Þr
þ
1 ðseSUGN3 Þ
ð18Þ 1r . r
b hSHOC jj$qh jjhSHOC ; 2ucha qref
ð21Þ
line plots y
ðsqSUPG ÞUGN ,
and h is the temperature. The parameters ðsuSUPG ÞUGN and ðseSUPG ÞUGN are calculated from their components by using the ‘‘r-switch’’: 1r 1 1 ðsqSUPG ÞUGN ¼ þ ; ð17Þ r r ðsqSUGN1 Þ ðsqSUGN2 Þ 1r 1 1 1 þ þ ; ðsuSUPG ÞUGN ¼ r r r ðsuSUGN1 Þ ðsuSUGN2 Þ ðsuSUGN3 Þ
ð20Þ
Shock
M = 2.0
10˚
M = 1.64
29.3º
ð19Þ Typically, r = 2.
x
Fig. 1. Oblique shock. Problem description.
Fig. 2. Oblique shock. Meshes tested: quadrilateral elements (top) and triangular elements arranged in mesh types A (bottom left) and D (bottom right).
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
hSHOC ¼ hJGN ; hJGN ¼ 2
nen X
ð22Þ
!1 jj $N a j
ð23Þ
.
a¼1
Here ucha is a characteristic velocity (such as kuhk or acoustic speed c or a reference velocity uref), and uint is an intrinsic velocity (such as kuhk or c or ucha). Typically, uint = ucha = uref. The reference density qref is defined as bR =2 q qref ¼ qinf sca ; ð24Þ qinf where qinf is the density at the inflow and qsca is a scaling density (such as qinf or the density behind a normal shock corresponding to the inflow Mach number). The parameter bR is set as bR = 1 for smoother shocks and
151
bR = 2 for sharper shocks (in return for tolerating possible overshoots and undershoots). The parameter b is set as b = 1 for smoother shocks and b = 2 for sharper shocks. The compromise between the b = 1 and b = 2 selections is defined as the following averaged expression for sSHOC: 1 sSHOC ¼ ððsSHOC Þb¼1 þ ðsSHOC Þb¼2 Þ. 2
ð25Þ
As another option for calculating the shock-capturing parameter, mSHOC is defined as mSHOC
!b=21 b nsd X 1 oUh 2 hSHOC ¼ jjY Zjj ; Y ox 2 i i¼1 1
ð26Þ
1.6 CYZ12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y 1.6 CYZU12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y
Fig. 3. Oblique shock. Computed with the mesh made of quadrilateral elements. Density along x = 0.9, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
152
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
where Y is a diagonal scaling matrix constructed from the reference values of the components of U: 3 2 ðU 1 Þref 0 0 0 0 6 0 ðU 2 Þref 0 0 0 7 7 6 7 6 7; 0 0 ðU Þ 0 0 Y¼6 3 ref 7 6 7 6 0 0 ðU 4 Þref 0 5 4 0 0
0
0
0
and b = 1 or b = 2. In a variation of the expression given by Eq. (26), mSHOC is defined by the following expression:
mSHOC
nsd X 1 oUh 2 Y ¼ jjY Zjj ox jjY1 Uh jj1b
ð27Þ h
Z¼
h
oU oU þ Ahi ot oxi
ð28Þ
OR Z ¼ Ahi
oUh ; oxi
ð29Þ
i
i¼1
ðU 5 Þref
!b=21
1
hSHOC 2
b ð30Þ
.
The compromise between the b = 1 and b = 2 selections is defined as the following averaged expression for mSHOC: 1 ðmSHOC Þb¼1 þ ðmSHOC Þb¼2 . ð31Þ mSHOC ¼ 2
1.6 CYZ12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y 1.6 CYZU12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y
Fig. 4. Oblique shock. Computed with the mesh made of triangular elements arranged in mesh type A. Density along x = 0.9, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
153
1.6 CYZ12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y 1.6 CYZU12 τ82-MOD with δ91 Exact solution 1.5
Density
1.4
1.3
1.2
1.1
1
0.9 0
0.2
0.4
0.6
0.8
1
y
Fig. 5. Oblique shock. Computed with the mesh made of triangular elements arranged in mesh type D. Density along x = 0.9, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
y
M = 2.3
78
29º 2 M = 2.9 line plots
M = 1.942
Shock
3
1 23º
x
Fig. 6. Reflected shock. Problem description.
We also propose versions of mSHOC that take into account the Mach number and shock intensity across a
shock. In doing that, we modify mSHOC given by Eqs. (26) and (30) as follows:
154
mSHOC
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
mSHOC
! b jj$qh jjhSHOC J 1=bM 1þ hM 1i . qref ð32Þ
Here M is the Mach number, ‘‘h i’’ is the Macaulay bracket: 0; x 6 y; hx yi ¼ ð33Þ x y; x > y; and the parameters bJ and bM can each be set to 1 for smoother shocks and 2 for sharper shocks. Based on Eq. (26), a separate mSHOC can be calculated for each component of the test vector-function W: !b=21 b nsd X 1 oUh 2 hSHOC 1 ðmSHOC ÞI ¼ jðY ZÞI j ; Y ox 2 i I i¼1 Similarly, a separate mSHOC for each component of W can be calculated based on Eq. (30): nsd X 1 oUh 2 Y ðmSHOC ÞI ¼ jðY ZÞI j oxi I i¼1 I ¼ 1;2;...nsd þ 2.
!b=21 jðY1 Uh ÞI j1b
b hSHOC ; 2
ð35Þ
Given mSHOC, the shock-capturing term is defined as nel Z X S SHOC ¼ $Wh : ðjSHOC $Uh Þ dX; ð36Þ e¼1
mSHOC
2
mSHOC switchðsSUPG ðj uÞ ;
sSUPG ðjj uj cÞ2 ; mSHOC Þ;
ð37Þ
where the ‘‘switch’’ function is defined as the ‘‘min’’ function or as the ‘‘r-switch’’ used earlier in this section. For viscous flows, the above modification would be made separately with each of sqSUPG , suSUPG and seSUPG , and this would result in mSHOC becoming a diagonal matrix even if the option given by Eq. (34) or Eq. (35) is not exercised.
ð34Þ
I ¼ 1; 2; . . . nsd þ 2.
1
option given by Eq. (34) or Eq. (35) is exercised, then mSHOC becomes an (nsd + 2) · (nsd + 2) diagonal matrix, and the matrix jSHOC becomes augmented from an nsd · nsd matrix to an (nsd · (nsd + 2)) · ((nsd + 2) · nsd) matrix. To preclude compounding, mSHOC can be modified as follows:
Xe
where jSHOC is defined as jSHOC = mSHOC I. As a possible alternative, it is defined as jSHOC = mSHOC jj. If the
5. Test computations The test computations were carried out by using the space–time SUPG formulation described in Section 3.2. In the option denoted by ‘‘CYZ12’’, we use Eq. (26) with Eq. (31), and in the option denoted by ‘‘CYZU12’’, Eq. (30) with Eq. (31). For smoother shocks (see discussion in Section 4, in the option denoted by ‘‘CYZ1’’, we use Eq. (26) with b = 1. This is equivalent to the option denoted by ‘‘CYZU1’’, where we use Eq. (30) with b = 1. In all these options, we use for Z the expression given by Eq. (29), and set jSHOC = mSHOC I. Also, in all these options, as stabiliza-
Fig. 7. Reflected shock. Meshes tested: quadrilateral elements (top) and triangular elements arranged in mesh types D1A2 (middle) and A1D2 (bottom).
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
tion parameters, we use Eq. (11), and in Eqs. (17)–(19) we do not include sSUGN2. The results are compared to those obtained with the s82-MOD and d91 combination. The version of s82-MOD used in this paper for comparison is similar to the one given in [17]: s82MOD ¼ maxð0; sa sd Þ; ð38Þ where h $jjUh jj hBGN d91 ; sa ¼ ; sd ¼ ; u ¼ c þ u cc 2 2ucc jj$jjUh jjjj ðucc Þ
hBGN
For completeness, we also provide here the expression for d91 [7]: oUh d91 ¼ Ahk ox k
, 1
A0 e
2 !12 nsd X onj h oU oxk ox 1 ; k A0 j¼1 e
ð41Þ
e 0 is the where njÕs are the element coordinates, and A Jacobian of the transformation from the entropy variables to the conservation variables. As test problems, we consider steady-state, inviscid cases. The solutions are obtained with time-marching to the steady-state. At each time step, the coupled nonlinear equations involved are solved with the Newton–Raphson method. An iterative technique with
ð39Þ
!1 nen h X $jjU jj ¼2 . jj$jjUh jjjj $N a
155
ð40Þ
a¼1
3 CYZ12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
2.5
3
3.5
4
x 3 CYZU12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2 x
Fig. 8. Reflected shock. Computed with the mesh made of quadrilateral elements. Density along y = 0.25, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
156
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
5.1. Oblique shock
nodal-block-diagonal preconditioner and GMRES update method [18] is employed for solving the linear equation system involved at each Newton–Raphson iteration. We compute two well-known test problems: ‘‘oblique shock’’ and ‘‘reflected shock’’. These were used in many earlier publications, and here we compute each of them with options CYZ12 and CYZU12 and with three different meshes. For each test problem, we use a mesh with quadrilateral elements and two triangular meshes where the elements are arranged in two different ways. We compute 100 time steps with a time-step size of 0.05. The number of Newton–Raphson iterations and the outer and inner GMRES iterations are 3, 1, and 5, respectively.
Fig. 1 shows the problem description. This is a Mach 2 uniform flow over a wedge at an angle of 10 with the horizontal wall. The solution involves an oblique shock at an angle of 29.3 emanating from the leading edge. The computational domain is a square with 0 6 x 6 1 and 0 6 y 6 1. The inflow conditions are given as M = 2.0, q = 1.0, u1 = cos 10, u2 = sin 10, and p = 0.179. This results in an exact solution with the following outflow data: M = 1.64, q = 1.46, u1 = 0.887, u2 = 0.0, and p = 0.305. All essential boundary conditions are imposed at the left and top boundaries, slip condition at the wall, and no boundary conditions at the right boundary. The mesh is essentially
3 CYZ12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
2.5
3
3.5
4
x 3 CYZU12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2 x
Fig. 9. Reflected shock. Computed with the mesh made of triangular elements arranged in mesh type D1A2. Density along y = 0.25, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159
uniform and consists of 21 · 21 nodes. Fig. 2 shows the three different meshes tested. Figs. 3–5 show, for three different meshes, the density along x = 0.9, obtained with CYZ12 and CYZU12, compared with the solution obtained with the s82-MOD and d91 combination. In addition to being much simpler, the new shock-capturing parameters yield shocks with better quality. The improvement over the s82-MOD and d91 combination is even more substantial with the triangular meshes.
157
computational domain is a rectangle with 0 6 x 6 4.1 and 0 6 y 6 1. The inflow conditions in R1 are given as M = 2.9, q = 1.0, u1 = 2.9, u2 = 0.0, and p = 0.7143. Specifying these conditions and requiring the incident shock to be at an angle of 29 results in an exact solution with the following flow data: R2: M = 2.378, q = 1.7, u1 = 2.619, u2 = 0.506, and p = 1.528; R3: M = 1.942, q = 2.687, u1 = 2.401, u2 = 0.0, and p = 2.934. All essential boundary conditions are imposed at the left and top boundaries, slip condition at the wall, and no boundary conditions at the right boundary. The mesh is uniform and consists of 61 · 21 nodes. Fig. 7 shows the three different meshes tested. Figs. 8–10 show, for three different meshes, the density along y = 0.25, obtained with CYZ12 and CYZU12, compared with
5.2. Reflected shock Fig. 6 shows the problem description. This problem consists of three flow regions (R1, R2 and R3) separated by an oblique shock and its reflection from the wall. The
3 CYZ12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
2.5
3
3.5
4
x 3 CYZU12 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2 x
Fig. 10. Reflected shock. Computed with the mesh made of triangular elements arranged in mesh type A1D2. Density along y = 0.25, obtained with CYZ12 (top) and CYZU12 (bottom), compared with the solution obtained with the s82-MOD and d91 combination.
158
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159 3 CYZ1 τ82-MOD with δ91 Exact solution 2.5
Density
2
1.5
1
0.5
0 0
0.5
1
1.5
2
2.5
3
3.5
4
x
Fig. 11. Reflected shock. Computed with the mesh made of triangular elements arranged in mesh type A1D2. Density along y = 0.25, obtained with CYZ1 (equivalent to CYZU1), compared with the solution obtained with the s82-MOD and d91 combination.
the solution obtained with the s82-MOD and d91 combination. Again, the new, much simpler shock-capturing parameters yield shocks with better quality, and the improvement over the s82-MOD and d91 combination is even more substantial with the triangular meshes. The small undershoot seen in the solutions obtained with CYZ12 and CYZU12 and mesh type A1D2 can be eliminated by using CYZ1, which is equivalent to CYZU1. Fig. 11 shows, for mesh type A1D2, the density along y = 0.25, obtained with CYZ1, compared with the solution obtained with the s82-MOD and d91 combination. While the small overshoot has been eliminated at the expense of losing some of the sharpness of the shocks, they are still sharper compared to the solution obtained with the s82-MOD and d91 combination.
YZb Shock-Capturing and an earlier shock-capturing term that we call ‘‘d91’’. In addition to being much simpler, YZb Shock-Capturing yields substantially better shock quality, and we see even more substantial improvements for triangular elements.
Acknowledgements This work was supported by the US Army Natick Soldier Center (Contract No. DAAD16-03-C-0051) and NASA Johnson Space Center (Grant No. NAG91435).
References 6. Concluding remarks In the context of test computations with inviscid supersonic flows, we provided a comprehensive assessment of the stabilization and shock-capturing parameters that we introduced recently for the Streamline-Upwind/ Petrov–Galerkin (SUPG) formulation of compressible flows based on conservation variables. We focused on performance evaluation of the shock-capturing parameters that we categorized as ‘‘YZb Shock-Capturing’’ in this paper. We used two well-known test problems: ‘‘oblique shock’’ and ‘‘reflected shock’’. For each test problem, we used a mesh with quadrilateral elements and two triangular meshes where the elements are arranged in two different ways. With emphasis on the shock quality, we compared the solutions obtained with
[1] Hughes TJR, Brooks AN. A multi-dimensional upwind scheme with no crosswind diffusion. In: Hughes TJR, editor. Finite element methods for convection dominated flows, AMD-vol. 34. New York: ASME; 1979. p. 19–35. [2] Brooks AN, Hughes TJR. Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput Meth Appl Mech Eng 1982;32:199–259. [3] Tezduyar TE, Hughes TJR. Development of time-accurate finite element techniques for first-order hyperbolic systems with particular emphasis on the compressible Euler equations. NASA Technical Report NASA-CR-204772, NASA, 1982. Available from: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/ 19970023187_1997034954.pdf. [4] Tezduyar TE, Hughes TJR. Finite element formulations for convection dominated flows with particular emphasis on the compressible Euler equations. In: Proceedings of AIAA 21st aerospace sciences meeting, AIAA Paper 83-0125, Reno, Nevada, 1983.
T.E. Tezduyar, M. Senga / Computers & Fluids 36 (2007) 147–159 [5] Hughes TJR, Tezduyar TE. Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations. Comput Meth Appl Mech Eng 1984;45:217–84. [6] Hughes TJR, Franca LP, Mallet M. A new finite element formulation for computational fluid dynamics: VI. Convergence analysis of the generalized SUPG formulation for linear timedependent multi-dimensional advective–diffusive systems. Comput Meth Appl Mech Eng 1987;63:97–112. [7] Le Beau GJ, Tezduyar TE. Finite element computation of compressible flows with the SUPG formulation. Advances in finite element analysis in fluid dynamics, FED-vol. 123. New York: ASME; 1991. p. 21–7. [8] Tezduyar TE, Park YJ. Discontinuity capturing finite element formulations for nonlinear convection–diffusion–reaction equations. Comput Meth Appl Mech Eng 1986;59:307–25. [9] Hughes TJR, Mallet M. A new finite element formulation for computational fluid dynamics: III. The generalized streamline operator for multidimensional advective–diffusive systems. Comput Meth Appl Mech Eng 1986;58:305–28. [10] Shakib F, Hughes TJR, Johan Z. A new finite element formulation for computational fluid dynamics: X. The compressible Euler and Navier–Stokes equations. Comput Meth Appl Mech Eng 1991;89:141–219. [11] Tezduyar TE. Adaptive determination of the finite element stabilization parameters. In: Proceedings of the ECCOMAS computational fluid dynamics conference 2001 (CD-ROM), Swansea, Wales, United Kingdom, 2001.
159
[12] Tezduyar TE. Computation of moving boundaries and interfaces and stabilization parameters. Int J Numer Meth Fluids 2003;43:555–75. [13] Tezduyar TE. Finite element methods for fluid dynamics with moving boundaries and interfaces. In: Stein E, De Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics, vol. 3: Fluids. John Wiley & Sons; 2004. Chapter 17. [14] Tezduyar TE. Determination of the stabilization and shockcapturing parameters in SUPG formulation of compressible flows. In: Proceedings of the European congress on computational methods in applied sciences and engineering, ECCOMAS 2004 (CD-ROM), Jyvaskyla, Finland, 2004. [15] Tezduyar TE, Senga M. Stabilization and shock-capturing parameters in SUPG formulation of compressible flows. Comput Meth Appl Mech Eng, in press. [16] Tezduyar TE, Behr M, Liou J. A new strategy for finite element computations involving moving boundaries and interfaces–the deforming-spatial-domain/space–time procedure: I. The concept and the preliminary numerical tests. Comput Meth Appl Mech Eng 1992;94:339–51. [17] Aliabadi SK, Tezduyar TE. Parallel fluid dynamics computations in aerospace applications. Int J Numer Meth Fluids 1995;21: 783–805. [18] Saad Y, Schultz M. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Scient Statist Comput 1986;7:856–69.