Solving the Navier-Stokes Equations by the CESE Method
AIAA-2004-0075
Solving the Navier-Stokes Equations by the CESE Method Moujin Zhang1, S.-T. John Yu2 Mechanical Engineering Department Ohio State University, Columbus, OH 43210 Sin-Chung Chang3 NASA Glenn Research Center Cleveland, OH 44135
Abstract The CESE method is a new numerical paradigm for nonlinear hyperbolic conservation laws. In this paper, we extend the CESE method for solving viscous flows. Both viscous and inviscid fluxes are incorporated into the space-time integration of the CESE method to enforce flux conservation locally and globally. To overcome excessive numerical damping, incurred by clustered mesh for resolving the boundary layer, a CFL number insensitive scheme, recently developed by Chang, is employed. The capabilities of the present scheme are demonstrated by numerical solutions of a shock/ boundary layer interaction problem and a driven cavity flow. Without preconditioning the flow equations, the present approach can calculate flows at all speed.
1. Introduction The space-time conservation element and solution element (CESE) method is a new numerical framework for solving hyperbolic conservation laws. Contrast to modern upwind schemes, no reconstruction procedure or Riemann solver is needed and the computational logics and operational counts of the CESE method are much simpler and more efficient. Numerous one-, two-, and three-dimensional, steady and unsteady results, − Research Associate, Email: [email protected] 2 Associate Professor, Email: [email protected] 3 Senior Aerospace Engineer, Email: [email protected] 1
Copyright 2004 by the American Institute of Aeronautics and Astronatics, Inc, All right reserved
obtained by using the CESE method, have be reported in the cited references [1-8], including flows with complex shock systems, aero acoustics, flow dominated by vortices, The present paper aims at solving the Navier-stokes equations for viscous flows. In the setting of the space-time integral formation of the CESE method, the inviscid and viscous fluxes are incorporated in an equal footing manner to enforce local and global flux conservation. Moreover, to overcome the over-dissipation problem when CFL is smaller than 0.1, we also employ the newly developed CFL insensitive scheme by Chang and coworkers [9, 10, 11]. We note that numerical dissipation of the original CESE method tends to increase with decreasing local CFL number. Because the CFL numbers may vary significantly across the computational domain, usually due to non-uniform mesh, numerical solutions may become overly dissipated. To overcome this problem, Chang [9] proposed the CFL insensitive scheme for the Euler equations in one spatial dimension. The scheme was further extended for Euler equations in two spatial dimensions [10, 11]. Numerical results shows that the scheme is CFL number insensitive and could crisply resolve shocks as well as contact discontinuity. In the present paper, the new CFL insensitive methods are employed to resolve the boundary layer in solving the Navier-stokes equations. The present paper is organized as follows. In Section 2, we present the CESE method for the Navier-Stokes equations in two spatial dimensions. In section 3, we review the CFL insensitive scheme. Section 4 shows the numerical results for viscous flows. In particular, the same computer code has been used for a supersonic flow with shock/boundary layer interactions and incompressible cavity flow. We then offer concluding remarks and provide cited references.
1 American Institute of Aerospace and Astronautics
Let x1 = x , x 2 = y and x3 = t be the coordinates of a three-dimensional Eucilidean space E 3 . Equation (2.1) becomes the divergence free condition in E 3 ,
2. Numerical Methods 2.1 Space-Time Integration The dimensionless unsteady Navier-stokes equations in two-spatial-dimensions can be expressed as ∂u m ∂f m ∂g m ∂f vm ∂g vm + + − − =0, ∂t ∂x ∂y ∂x ∂y for m = 1, 2, 3, 4, where
(2.1)
(u1 , u 2 , u 3 , u 4 )
= ( ρ , ρu , ρv , e ) ,
( f 1 , f 2 , f 3 , f 4 )T
= ρu , ρu 2 + p, ρuv, (e + p )u
(g1 , g 2 , g 3 , g 4 )T
T
(2.2a)
(
)
(
)
T
, (2.2b)
= ρv, ρvu , ρv 2 + p, (e + p )v , (2.2c)
( fν 1 , fν 2 , fν 3 , fν 4 )T
T
(
= 0, τ xx , τ xy , uτ xx + vτ xy − q x
)T
(2.2d) and
(gν 1 , gν 2 , gν 3 , gν 4 )T
(
= 0, τ xy , τ yy , uτ xy + vτ yy − q y
)T
(2.2e) where ρ, p, u and v are density, pressure, velocity components in x- and y-direction. The specific total energy e is e = p (γ − 1) + ρ u 2 + v 2 2 . γ is the specific heat ratio. The stress components in Eq (2.2d-e) are as follows,
(
τ xx
)
1 4 ∂u 2 ∂v = − Re 3 ∂x 3 ∂y
τ yy =
∫∇⋅h
=
m dV
∫h
m
⋅ ds = 0 ,
(2.4)
S (V )
V
for m = 1, 2, 3, 4, where S(V) is the boundary of an arbitrary space-time region V in E 3 and ds = ndσ , where dσ and n are the area and the outward unit normal vector of a surface element on S(V).
2.2 CE and SE The two-dimensional spatial domain is divided into no-overlapped triangles. Refer to Fig. 1. Point G, the centroid of ∆BDF, is marked by a solid circle, and A, C and E, centorids of ∆FMB, ∆BJD and ∆DLF, and are marked by hollow circles. In Fig. 2, A, B, C, D, E and F form a hexagon ABCDEF. Point G , marked by a solid square is the centroid of hexagon ABCDEF, and it is the solution point of ∆BDF. Let S q1 , S q 2 and S q 3 be areas of quadrilaterals ABGF, BCDG and DEFG, and x q1 , y q1 , x q 2 , y q 2 and x q 3 , y q 3 be the spatial
(
) (
)
(
)
coordinates of their centroids. The spatial coordinates of point G are xG = yG =
1 4 ∂v 2 ∂u − Re 3 ∂y 3 ∂x
(2.3)
where h m = ( f m − fνm , g m − gνm , u m ) are the current density vectors in E 3 . By using Gauss’ divergence theorem in E 3 , we have
u m , f m , g m fνm and gνm
are flow variables, inviscid fluxes and viscous fluexs in x- and y- directions, respectively, as T
∇ ⋅ hm = 0 ,
S q1 x q1 + S q 2 x q 2 + S q 3 x q 3 S q1 + S q 2 + S q 3 S q1 y q1 + S q 2 y q 2 + S q 3 y q 3
.
(2.5)
S q1 + S q 2 + S q 3
Points A , C and E , marked by hollow squares in Fig. 2, and are the solution points of ∆FMB, ∆BJD and ∆DLF, respectively. In the space-time domain, A, B, C, D, E, F and G are at the time level n − 1 2 , and A' B ' ,
1 ∂u ∂v + Re ∂y ∂x 1 ∂T qx = 2 (γ − 1)M RePr ∂x 1 ∂T qy = 2 (γ − 1)M RePr ∂y
τ xy =
C ' , D ' , E ' , F ' and G ' are at the time level n . Points A" , B " , C " , D " , E " , F " and G " are at the time level n + 1 2 . Let j, k and n be indices for x, y and t,
where Re is the Reynolds number, Pr the Prandtl number and M the freestream Mach number.
respectively. Points G ' , A, C and G are marked by ( j, k , n ) , ( j1, k1, n − 1 2) , ( j 2, k 2, n − 1 2) and ( j 3, k 3, n − 1 2) , respectively. Shown in Fig. 3, the solution points G ' , A , C and E are placed in a staggered positions in E 3 , and their coordinates are
2 American Institute of Aerospace and Astronautics
marked by ( j, k , n ) , ( j1, k1, n − 1 2 ) , ( j 2, k 2, n − 1 2 ) and ( j 3, k 3, n − 1 2 ) . Note that, a triangle’s centroid G
'
4
( f mx )nj = ∑ ( f m,l )nj,k (ulx )nj,k ,
and the associated solution point, G ' have different spatial coordinates. In the calculation, flow variables are stored at the solution points.
( f my )nj = ∑ ( f m,l )nj,k (uly )nj,k ,
As presented in Fig. 4, the solution element SE ( j, k , n ) associated with point G ' ( j, k , n ) , is the
(g mx )nj,k
'
'
'
'
'
and SE ( j 3, k 3, n − 1 2 ) associated with points A, C and E, respectively. The surfaces of the four SEs form three CEs for point G ' . Refer to Fig. 3. Three CEs are ABGFA' B ' G ' F ' , quadrilateral cylinders ' ' ' ' ' ' ' ' CDGBC D G B and EFGDE F G D , and are referred to as CE1 ( j, k , n ) , CE 2 ( j, k , n ) and CE3 ( j, k , n ) , respectively. CE ( j, k , n ) is the union of CE1 ( j, k , n ) , CE 2 ( j, k , n ) and CE 3 ( j, k , n ) .
Inside SE ( j, k , n ) , the first-order Taylor series expansion is employed to descritize the flow variables and inviscid fluxes: u m*
( x, y , t ; j , k , n ) = ( )
u m nj ,k
+
+(
)
u mx nj ,k
u mt nj ,k
n
(
f m* (x, y, t; j, k , n ) = ( f m )nj ,k + ( f mx )nj ,k x − x j ,k + g m*
)
( ) (y − y j,k ) + ( ) (t − t ) , f my nj ,k
f mt nj ,k
(x, y , t; j, k , n ) = ( )
g m nj ,k
+(
)
g mx nj ,k
n
( )nj,k (y − y j,k ) + (g mt )nj,k (t − t n ) ,
( j, k , n ) , we let (umt )nj,k = −( f mx )nj,k − (g my )nj,k ,
(2.7)
(g my )nj,k = ∑ (g m,l )nj,k (uly )nj,k .
(2.14)
and 4
l =1
Aided by Eqs. (2.11-14), Eq. (2.9) can be recast to 4
(u mt )nj
=−
∑(f )
n m ,l j , k
l =1
(g mt )nj,k ( f mt )nj
4
for m, l = 1, 2, 3, 4. Aided by the chain rule, we have,
∑ (g ) (u ) l =1
∑(f )
=
n m ,l j , k
n ly j ,k
, (2.15)
(ult )nj,k
n m ,l j , k
l =1
4
4
∑(f ) ∑(f ) n m ,l j
4
−
n l , r j ,k
r =1
4
(u rx )nj,k
∑ ( f ) ∑ (g ) (u ) n m ,l j
4
(g mt )nj ,k
=
n l , r j ,k
r =1
∑ (g ) l =1
n m ,l j ,k
n ry j ,k
,
(2.16)
(ult )nj ,k
4
4
l =1
r =1
∑ (g m,l )nj,k ∑ ( f l ,r )nj,k (u rx )nj,k 4
∑ (g ) ∑ (g ) (u ) l =1
(2.10)
−
can be expressed as,
4
entries of Jacobian matrixes F and G, i.e.,
4
(u lx )nj,k
Aided by the chain rule and Eqs. (2.15), ( f mt )nj ,k and
−
for m = 1, 2, 3, 4. To proceed, let f m , l and g m , l be the ∂f m ∂g and g m,l = m , ∂ul ∂ul
(2.13)
(2.8)
(2.9)
(ulx )nj,k
n m ,l j ,k
,
l =1
=−
for m = 1, 2, 3, 4. At point
f m,l =
∑ (g )
l =1
(x − x j,k )
+ g my
4
=
l =1
(2.6)
(2.12)
l =1
=−
(x − x j,k )
( ) (y − y j,k ) + ( ) (t − t ) , u my nj ,k
4
'
union of four planes, i.e., the hexagon A B C D E F , the quadrilaterals B " BGG " , D " DGG " and F " FGG " , and their immediate neighborhood. Similarly, there are three SEs, i.e., SE ( j1, k1, n − 1 2 ) , SE ( j 2, k 2, n − 1 2 )
(2.11)
l =1
n m ,l j ,k
r =1
n l , r j ,k
n ry j ,k
,
(2.17)
for m = 1, 2, 3, 4. Aided by Eqs. (2.11-17), Eqs. (2.6-8) could fully specify the distribution of
u m* , f m* and g m*
inside SE ( j, k , n ) when values of (u m )nj ,k (u mx )nj ,k and
(umy )nj,k
are known.
The above expressions are all concerned with the distribution of flow variables and inviscid fluxes inside a SE. For expressing the viscous fluxes, fνm and gνm ,
3 American Institute of Aerospace and Astronautics
with u m u mx and umy , let’s express ∂u ∂x and ∂v ∂y with um , umx and umy as, ∂u 1 ∂ (ρu ) ρu ∂ρ u2 x u2u1x = − 2 = − 2 , ∂x ρ ∂x u1 ρ ∂x u1
(2.18a)
∂v 1 ∂ (ρv ) ρv ∂ρ u3 y u3u1 y = − 2 = − 2 . u1 ∂y ρ ∂y u1 ρ ∂y
(2.18b)
and EDD ' E ' belong to SE ( j3, k 3, n − 1 2) . Let S be the area of the surface. Over each area, let the outward normal vector be n, and the surface vector s = n S. The flux leaving a surface is equal to the scalar product between the vector h *m = f m* − fν*m , g m* − gν*m , u m* , evaluated at surface’s centroid, and the surface vector s. For hexahedron A' B ' C ' D ' E ' F ' in E3, its surface vector is
(
With the aid of Eq. (2.18), τ xx can be expressred as
τ xx =
uu u 1 uu u − 4 2 1x + 4 2 x + 2 3 1 y − 2 3 y . 3Re u12 u12 u1 u1
(
Similar expressions can be obtained for τ yy , τ xy , qx and q y in Eq. (2.2d-e). Moreover the viscous fluxes, fνm and gνm , can be expressred by um , umx and umy . Aided with Eq. (2.6), the distribution of flow variables inside SE ( j, k , n ) can be obtained. With flow variables, the distribution of viscous fluxes, fν*m and gν*m , can be
( )
n
n j, k
.
( )nj,k ,
With known values of (u m )nj ,k , (u mx )nj ,k and u my *
*
the distribution of flow variable u m , inviscid fluxes f m *
and g m , and viscous fluxes fν*m and gν*m , can be fully specified
(
h *m
f m*
= − approximated by
∫
SE ( j, k , n )
inside fν*m , g m*
−
gν*m , u m*
h *m ⋅ ds = 0 ,
),
)
)
s A' B 'C ' D ' E ' F ' = 0, 0, S q1 + S q 2 + S q 3 ,
(2.19)
further expressed by (u m ) j , k (u mx )nj ,k and u my
CDD 'C ' CBB 'C ' CDGB , and belong ( ) SE j 2 , k 2 , n − 1 2 . Quadrilaterals EFGD , EFF ' E ' to
.
Let
and Eq. (2.4) can be
S (V )
(
)
and the coordinates of its centroid G ' are xG , y G , t n . '
'
'
'
'
'
The flux leaving the surface A B C D E F is
(FLUX m ) A B C D E F '
'
'
'
'
'
(
)
= S q1 + S q 2 + S q 3 (u m )nj ,k . (2.22)
Let’s calculate the fluxes leaving surfaces belonging to SE ( j1, k1, n − 1 2 ) . For quadrilateral ABB ' A' , its surface vector is s ABB ' A' =
∆t ( yB − y A , x A − xB , 0 ) , 2
(2.23)
and the coordinates of its centroid are x A + x B y A + y B n −1 2 ∆t , ,t + . The flux leaving the 2 2 4 surface ABB ' A' is
(FLUX m ) ABB A '
(2.20)
(2.21)
'
=
∆t ( y B − y A )( f m )nj1−,1k12 + x A + x B − x A ( f mx )nj1−,1k12 2 2 y + yB + A − y A f my 2
( )nj1−,1k12 + ∆4t ( f mt )nj1−,1k12 − fνAB m
for m = 1, 2, 3, 4. 2.3 Time Marching for u
By imposing Eq. (2.20) over CE ( j, k , n ) , an algebraic equation can be obtained based on the global flux conservation and (u m )nj ,k can be obtained directly for
each m = 1, 2, 3, 4. To proceed, we calculate the flux leaving surfaces of CE ( j, k , n ) . Consider CE ( j, k , n ) , the hexahedral cylinder ABCDEFA' B 'C ' D ' E ' F ' referring to Fig. 3. The surfaces of CE ( j, k , n ) can be divided into four groups. Hexahedron A' B ' C ' D ' E ' F ' belongs to SE ( j, k , n ) . Quadrilaterals ABGF , ABB ' A' and AFF ' A' belong to SE ( j1, k1, n − 1 2 ) . Quadrilaterals
+
∆t (x A − x B )(g m )nj1−,1k12 + x A + x B − x A (g mx )nj1−,1k12 2 2
( )nj1−,1k12 + ∆4t (g mt )nj1−,1k12 − gνABm
y + yB + A − y A g my 2
where the viscous fluxes,
(2.24)
,
fνmAB and gνAB m , are obtained
from flow variables at the centroid of surface ABB ' A'
( )nj1−,1k12 .
and (u mx )nj1−,1k 12 and u my
For quadrilateral AFF ' A' , its surface vector is
4 American Institute of Aerospace and Astronautics
∆t ( y A − yF , xF − x A , 0 ) , 2
s AFF ' A' =
(2.25)
and the coordinates of its centroid are x A + x F y A + y F n −1 2 ∆t + . The flux leaving the , ,t 2 2 4 surface AFF ' A' is
(FLUX m ) AFF A '
=
'
∆t ( y A − y F )( f m )nj1−,1k12 + x F + x A − x A ( f mx )nj1−,1k12 2 2
( )nj1−,1k12 + ∆4t (g mt )nj1−,1k12 − gνAFm ,
y + yF + A − y A g my 2
(2.26) AF where the viscous fluxes, fνAF m and gνm , are obtained
from flow variables at the centroid of surface AFF ' A'
( )nj1−,1k12 .
and (u mx )nj1−,1k 12 and u my
The flow variables at
For quadrilateral ABGF, its surface vector is
)
s ABGF = 0, 0, − S q1 ,
(2.27)
(
)
and the coordinates of its centroid are x q1 , y q1 , t n −1 2 . The flux leaving the surface ABGF is
(
+ x q1 − x G
)( )
(
+ y q1 − y G
)( )
u my nj1−,1k12
].(2.28)
( fluxm )1n −1 2 = − S q1 [(u m )nj1−,1k12 + (x q1 − x A )(u mx )nj1−,1k12
(
+ y q1 − y A
)( )
]
∆t ( y B − y F )( f m )nj1−,1k12 + 2
∆t + [(x B + x A )( y B − y A ) + (x F + x A )( y A − y F ) 4
]
− 2 x A ( y B − y F ) ( f mx )nj1−,1k12
[
∆t 2 x F − x B2 − 2 x A (x F − x B ) (g mx )nj1−,1k 12 4
+
∆t [(x A − x B )( y A + y B ) 4
− y F )( f mt )nj1−,1k12
[
]
( )nj1−,1k12
+
.
−
(∆t )2 (x 8
F
− x B )(g mt )nj1−,1k12
[
∆t AB ( yB − y A ) fνAB m + (x A − x B )gνm 2
]
, (2.29) + ( y A − y F ) fνmAF + ( x F − x A )gνAF m
( )nj1−,1k12
where (u m )nj1−,1k12 , (u mx )nj1−,1k12 and u my
are values
stored at solution point A . Similarly, the fluxes leaving three quadrilaterals belonging to conservation over SE ( j 2, k 2, n − 1 2 ) are
[
(
)
− S q 2 (u m )nj 2−1,k 22 + x q 2 − xC (u mx )nj 2−1,k 22
( flux m )n2 −1 2 =
)( )nj 2−1,k 22 ]
(
+
∆t ( y D − y B )( f m )nj 2−1,k 22 2
+
∆t [(x D + xC )( y D − y C ) 4
]
The flux leaving three surfaces belonging to SE ( j1, k1, n − 1 2 ) is the sum of Eqs. (2.24), (2.26) and (2.28). We have
u my nj1−,1k12
+
B
+ (x B + x C )( y C − y B ) − 2 x C ( y D − y B ) ( f mx )nj 2−1,k 22
[
= − S q1 (u m )nj1-1,k21
u mx nj1−,1k12
∆t (x F − x B )(g m )nj1−,1k12 2
8
+ y q 2 − y C u my
centroid are obtained following Eq. (2.6).
(FLUX m ) ABGF
+
x +x ∆t (x F − x A )(g m )nj1−,1k12 + A F − x A (g mx )nj1−,1k12 + 2 2
(
(∆t )2 ( y
+ (x F − x A )( y A + y F ) − 2 y A (x F − x B )] g my
( )nj1−,1k12 + ∆4t ( f mt )nj1−,1k12 − fνAF m
y + yA + F − y A f my 2
+
]( )
∆t 2 + y B − y F2 − 2 y A ( y B − y F ) f my 4
n −1 2 j1,k 1
[
]( )
+
∆t 2 y D − y B2 − 2 y C ( y D − y B ) f my 4
+
(∆t )2 ( y
+
∆t (x B − x D )(g m )nj 2−1,k 22 2
+
∆t 2 x B − x D2 − 2 x C (x B − x D ) (g mx )nj 2−1,k 22 4
+
∆t [(x B − xC )( y B + yC ) 4
8
D
n −1 2 j 2,k 2
− y B )( f mt )nj 2−1,k 22
[
]
]( )nj1−,1k12
+ (xC − x D )( y D + y C ) − 2 y C (x B − x D ) g my
+
(∆t )2 (x 8
B
5 American Institute of Aerospace and Astronautics
− x D )(g mt )nj1−,1k 12
[
∆t CD ( y D − yC ) fνCD m + (xC − x D )gνm 2
−
(um )nj,k
]
CB + ( yC − y B ) fνCB m + ( x B − xC )gνm .
The fluxes leaving SE ( j 3, k 3, n − 1 2 ) are
the
three
[
( flux m )3n −1 2
)( )nj3−,1k 32 ]
2.4 Time Marching for ux and uy In this section, we illustrate the calculation of the spatial derivatives of the flow variables, i.e.,
)
(u mx )nj,k and (u my )nj,k . We note that point G '
∆t ( y F − y D )( f m )nj 3−,1k 32 + 2 ∆t [(x F + x E )( y F − y E ) 4
+ (x D + x E )( y E − y D ) − 2 x E ( y F − y D )]( f mx )nj 3−,1k 32
+
[
]( )
∆t 2 y F − y D2 − 2 y E ( y F − y D ) f my 4
(∆t )2 ( y + 8
− y D )(
F
n −1 2 j 3,k 3
+
∆t 2 x D − x F2 − 2 x E (x D − x F ) (g mx )nj 3−,1k 32 4
y E* = (3yG + 2 y E − y A − yC ) 3
[
]
( )nj3−,1k 32
− x F )(
)
]
,
,
fνEF m
,
,
gνEF m
,
fνED m
and
gνED m
fνCD m
,
, are
For flux conservation over CE ( j, k , n ) , we have '
'
'
'
'
'
(
(u )
(2.31)
AF obtained in the similar way of that for fνAF m and gνm .
(FLUX m ) A B C D E F
n −1 2 j1,k 1
(2.35)
.
(2.36)
n −1 2 mx j1,k 1 n j 2 ,k 2
+ ( flux m )1n −1 2
+ ( flux m )n2 −1 2 + ( flux m )3n −1 2 = 0 . (2.32)
' n m j 3,k 3
( )
∆t u mt 2
+
n −1 2 j1,k 1
)( )nj1−,1k12 . (2.37)
(
+ y A* − y A u my
( ) on (u )
and u m'
n j 3,k 3
, at C * and E * can
' n m j1,k 1
be obtained. Based
For each m =1, 2, 3, 4. The viscous fluxes, gνCB m
( ) )(u )
= um
( )
[
fνCB m
' n m j1,k 1
Similarly, u m'
∆t EF ( y F − y E ) fνEF − m + (x E − x F )gνm 2
gνCD m
(u )
+ x A* − x A
g mt nj 3−,1k 32
ED + ( y E − yD ) fνED m + ( xD − xE )gνm .
,
from Eq. (2.33). Flow variable at point A* is obtained by a first-order Taylor series, Eq. (2.6), as
+ (x E − x F )( y E + y F ) − 2 y E (x D − x F )] g my D
(2.34)
Flow variables at G ' , i.e., (u m )nj ,k , are calculated
∆t + [(x D − x E )( y D + y E ) 4
8
,
y A* = (3 y G + 2 y A − y C − y E ) 3
x E* = (3xG + 2x E − x A − xC ) 3
∆t (x D − x F )(g m )nj 3−,1k 32 2
(∆t )2 (x +
x A* = (3x G + 2 x A − x C − x E ) 3
y C * = (3 y G + 2 y C − y A − y E ) 3
)
is not the
centroid of ∆A 'C ' E ' unless a uniform mesh is used. As shown in Fig. 5, a triangle ∆A*C * E * , whose centroid is point G ' , is obtained by parallel moving ∆A 'C ' E ' in the spatial domain. The vertices’ coordinates of ∆A*C * E * are
xC * = (3xG + 2 xC − x A − x E ) 3
f mt nj 3−,1k 32
+
.
(2.33)
+ y q 3 − y E u my
+
S q1 + S q 2 + S q 3
of
= − S q 3 (u m )nj 3−,1k 32 + x q 3 − x E (u mx )nj 3−,1k 32
(
( fluxm )1n−1 2 + ( fluxm )n2−1 2 + ( flux m )3n−1 2
(2.30)
surfaces
(
=−
(u )
' n m j 2 ,k 2
,
and
on points A* , C * and E * , we apply central
( )nj,k
differencing to calculate (u mx )nj ,k and u my
at point
G ' , i.e.,
(u )
c n mx j ,k
(
=
1 2 S A*C *E *
)( )nj1,k1
(
y * − y * u' m E C
)( )nj 2,k 2 + (y A
+ y E * − y A* u m'
Aided by Eq. (2.22), an explicit expression for (u m )nj ,k can be obtained from above equation as
6 American Institute of Aerospace and Astronautics
*
)( )nj3,k 3 ,(2.38)
− y C * u m'
(u )
c n my j ,k
=
)( )nj1,k1
(
1 2 S A*C *E *
)( )nj 2,k 2 + (xC
(
level n = 0, 1, 2, ... , and hollow squares are at n = 1 2 , 3 2 , 5 2 , ... . According to the CESE method,
x * − x * u' m C E
+ x A* − x E * u m'
*
)( )nj3,k 3 ,(2.39)
the flow variables at solution point G ' (j, k, n) are determined by those at seven solution points A , B , C , D , E , F and G at the time level n-1. The hexagon A B C D E F is the numerical domain of dependence for the solution at G ' at time level n-1.
− x A* u m'
where S ∆A*C *E * is the area of ∆A*C * E * : S ∆A*C * E * =
(
1 x * y * + x C * y E * + x E * y A* 2 A C
)
− x A* y E * − xC * y A* − x E * y C * .
(2.40)
Similar central differencing can be applied to calculate
(u ( ) ) (u ( ) )
( )
1 n mx j ,k
(1) and u my
2 n my j ,k
for ∆A*G ' E *
n j ,k
( ) , and (u ( ) ) and (u ( ) ) (2 ) for ∆C * E *G ' , u mx 3 n mx j ,k
n j ,k
and
3 n my j ,k
for
∆A*C *G ' . Based on the re-weighting procedure in [4],
( )nj,k
we calculate (u mx )nj ,k and u my
(u mx )nj,k
( )
w = u mx
as,
( )nj,k = (umyw )nj,k ,
n j ,k
and u my
(2.41)
where
( )
w n umx j ,k
+ (θ m1θ m 3 )
α
( )
w n u my j ,k
( ) u (2 )
n mx j ,k
n j ,k
+ (θ m1θ m 2 )
α
( )
1 (1) = (θ m 2θ m 3 )α u my ω
+ (θ m1θ m 3 )
α
( ) u (2 )
n
my j ,k
( )
n mx j ,k
u (3)
, (2.42)
j ,k
+ (θ m1θ m 2 )
α
( ) u (3)
n
my j ,k
, (2.43)
θ mr
( ) ( ) 2
(r ) + u my
2
n
, j ,k
(2.44)
ω = (θ m1θ m 2 )α + (θ m 2θ m 3 )α + (θ m1θ m 3 )α .
(2.45)
The above CESE schemes are stable for CFL number < 1., while α≥0. Equations (2.33) and (2.41) form the commonly used CESE scheme, the a-α scheme. Hereafter we refer it as the original scheme.
3. CFL Insensitive Scheme 3.1 The CFL Condition The CFL number in two spatial dimensions is defined hereafter. The spatial projections of solution points are presented in Fig. 6, in which solid squares are at time
(3.1)
Let line segment G H be the distance from point G to boundary A B , we have
G H = 2 S ∆A B G
n
and
(r ) = u mx
G O , we have
G O = ∆t u 2 + v 2 , and θ 0 = arctan (v u ) .
( )
1 (1) = (θ m 2θ m 3 )α umx ω
Figure 7 shows the projection of a Mach cone on the spatial plane at t = (n-1)∆t, with point G ' being its vertex. The result is a circle with a radius of c∆t on the plane with O as the center of the circle. The boundary and interior of the circle form the domain of dependence for the flow solution at G ' . Thus, the CFL condition for the solution at G ' is defined such that if and only if the domain of dependence, i.e., the circle, lies in the interior of the hexagon A B C D E F . Let u, v and c be velocity components and sonic speed at solution point G ' at time level n-1. As shown in Fig. 7, for the velocity vector
θ1 = arctan
(x B − x A )2 + (y B − y A )2 , and
x A − xB yB − y A
.
(3.2)
As shown in Fig. 7, we choose a point P on the circle such that the line segment OP is parallel to line segment G H . Let R and S be the projection of O and P on G H . Obviously P is the closest point to the boundary A B on the circle. To keep the circle inside the hexagon, with A B as one of the boundary segments, we require that
v (1) = G S G H < 1 ,
(3.3)
where G S = ∆t c + u 2 + v 2 cos(θ 1 − θ 0 ) .
(3.4)
For the other five boundary line segments, we have similar conditions, i.e., v (2 ) for B C , v (3) for C D ,
v (4 ) for D E , v (5) for E F , and CFL condition is that,
{
v (6 ) for F A . The
}
v e = max v (1) , v (2 ) , v (3) , v (4 ) , v (5) , v (6 ) .
7 American Institute of Aerospace and Astronautics
(3.5)
ve < 1 .
(3.6)
Essentially, Eqs. (3.5-6) specify that the domain of dependence of the flow solution at G ′ must lies within its numerical domain of dependence, i.e., hexagon A B C D E F . In computation, ∆t is chosen to satisfy Eq. (3.5-6).
3.2 The CFL Insensitive Scheme According to [9-11], the CFL insensitive scheme is constructed such that (i) it would reduce to the original a scheme in the limit of CFL→0, and (ii) it would reduce to the original a-α scheme in the limit of CFL→1. As shown in Fig. 8, points Q1, Q2 and Q3 are centroids of quadrilaterals ABGF, BCDG and DEFG, points P1, P2 and P3 are defined within line segments A′Q1 , C ′Q2 and E ′Q3 as, x p1 = v e x A ′ + (1 − v e )x q1
y p1 = v e y A ′ + (1 − v e ) y q1
x p2 = v e xC ′ + (1 − v e )x q 2
y p3 = v e y E ′ + (1 − v e ) y q 3
(u )
' m p1*
( )
n −1 2 j1,k 1
= um
[(
,
(3.7)
,
(3.8)
,
(3.9)
Take Eq. (3.7) for instance, with CFL number, ve ,
1
(u )
' m p2*
( )
n −1 2 j 2,k 2
= um
[(
2
i.e., point G ′ , can be obtained by parallel translation of ∆A ' C ' E ' . The coordinates of the vertexes of ∆P1* P2* P3* are
x y
p1*
p2* p2*
(
)
(
3
p3*
)
1 − ve = v e xC * + 3xG ′ + 2 x q 2 − x q1 − x q 3 3 ,(3.11) 1 − ve = v e yC * + 3 y G ′ + 2 y q 2 − y q1 − y q 3 3
= ve y E *
(
(
∆t (u mt )nj 2−,1k 22 2
+
)
(u )
' m p3*
( )
n −1 2 j 3,k 3
= um
[(
+
(
)(
2
∆t (u mt )nj 3−,1k 32 2
)
(
)(
+ x p* − x E (u mx )nj 3−,1k 32 + y p* − y E u my 3
,
)nj 2−,1k 22 ]
3
)nj 3−,1k 32 ]
(3.15) P3* and G ′ known, the CFL insensitive scheme can be constructed. There are several versions of CFL insensitive scheme as stated in [9-11]. In the present paper, the Scheme II in [10] is applied. For completeness, the equations are provided in the following.
(u )
w n mx j ,k
≈
[
[
)
]( ) ](u ( ) ) + [1 + f (v )(s
1 n (1) 1 + f (v e )(s m1 ) j ,k u mx e
+ 1 + f (v e )(s m 2 )nj ,k
n j ,k
2 n mx j ,k
e
)
n m 3 j ,k
](u ( ) )
3 n mx j ,k
(3.16) w n my j , k
[
[
]( ) ](u ( ) ) + [1 + f (v )(s
1 (1) ≈ 1 + f (v e )(s m1 )nj ,k u my e 2 n my j , k
n
j ,k
e
)
n m3 j ,k
](u ( ) )
3 n my j , k
(3.17) where
[
]
e = 3 + f (ve )(s m1 )nj ,k + (s m 2 )nj ,k + (s m3 )nj ,k ,
(s mr )nj ,k
=
φ mr − 1, min(φ m1 , φ m 2 , φ m3 )
(3.18) (3.19)
φm1 = θ m 2θ m 3 , φm 2 = θ m1θ m 3 and φm 3 = θ m 2θ m1 , (3.20) and f (v e ) = 1 v e .
)
1 − ve 3 xG ′ + 2 x q 3 − x q1 − x q 2 3 .(3.12) 1 − ve + 3 y G ′′ + 2 y q 3 − y q1 − y q 2 3
x p* = v e x E * + y
)
1 − ve 3xG ′ + 2 x q1 − x q 2 − x q 3 = ve x A* + 3 ,(3.10) 1 − ve = ve y A* + 3 y G ′ + 2 y q1 − y q 2 − y q 3 3
(
1
)nj1−,1k12 ]
(3.14)
+ 1 + f (v e )(s m 2 )nj ,k
y
)(
+ x p* − x C (u mx )nj 2−1,k 22 + y p* − y C u my
Similar to that in section 2.4, the centroid of ∆P1* P2* P3* ,
1
(
(3.13)
(u )
(
)
+ x p* − x A (u mx )nj1−,1k12 + y p* − y A u my
changing from 0 to 1, point P1 moves from point Q1 to A ′ . Similar to the fact that the centroid of ∆A ' C ' E ' may not be G ′ , the centroid of ∆P1 P2 P3 may not be G ′ .
x p*
∆t (u mt )nj1−,1k12 2
+
With the values of flow variables at points P1* , P2* ,
y p2 = v e y C ′ + (1 − v e ) y q 2
x p3 = v e x E ′ + (1 − v e )x q 3
Flow variables at P1* , P2* and P3* are obtained by
)
8 American Institute of Aerospace and Astronautics
(3.21)
4. Numerical Results Two problems are solved by using the above CESE method for the Navier-Stokes equations.
4.1 Shock Wave Boundary Layer Interaction The shock wave boundary layer interaction problem in [12] is used as the test case. This problem is often used as a standard test for Navier-Stokes solver. As shown in Fig. 9, the computational domain is (x, y) ∈ [0, 2.4]×[0, 1.164]. The left boundary is defined as the inlet boundary where specified boundary condition is employed. The flow on the top boundary is specified to form an oblique shock, impinging on the wall. The right boundary is a supersonic outlet, where non-reflective boundary condition is used. The bottom boundary consists of a symmetric boundary and a solid wall, whose lengths are 0.8 and 1.6 respectively. For the solid wall, the no-slip boundary condition is employed. The incoming shock wave, emanating from the upper-left corner of the computational domain, impinges on the solid wall with an angle of 32.6o with respect to the wall. The flow Mach number on the left inlet boundary is 2.0. The flow condition on the top boundary is calculated based on the oblique shock condition to form the desired shock wave angle. We consider the viscous flow with Re=296,000. The computational domain is covered by 80,000 triangles. To resolve boundary layer, grid points are clustered to solid wall. The first grid point to the solid wall is 5×10-5 of the domain height. Due to the non-uniform mesh, the local CFL number varies significantly across the computational domain, the difference in CFL numbers is around 60 times. Figure 10 shows the pressure contours. There are 50 contours ranging equally from 0.01 to 0.25. The pressure contours obtained by CFL insensitive scheme are shown in Fig. 10a. Pressure contours by the original scheme are shown in Fig. 10b. The CFL insensitive scheme can crisply capture shocks, while the original scheme shows excessive numerical diffusion due to small CFL number. For viscous flow simulation, the CFL insensitive scheme plays an important role in capturing sharp shocks and detailed flow structure across the domain. Because the impinging shock is strong, boundary layer separation occurs at the shock impinging point at the wall. Figure 11 shows the velocity vectors in the recirculation zone. The change of boundary layer thickness and a boundary layer separation can be observed. Figures 12-13 show the pressure and the friction coefficient along solid wall, respectively. Curves are the calculated results and symbols are the
experimental data. The present results are in good agreement with the experimental data.
4.2 Driven Cavity Flow The driven cavity flow is also a benchmark problem for testing incompressible viscous solvers [13]. The CESE method is employed for this incompressible flow calculation. The computational domain is (x, y) ∈ [0, 1]×[0, 1.]. The top boundary is a moving solid wall, and the other three boundaries are stationary walls. We conducted computation with Re=1,000. 12,000 triangles are used to cover the computation domain. The distributions of velocity components along centerlines are plotted in Fig. 14 and 15. The velocity vector field is plotted in Fig. 16. The present results are in good agreement with the data reported in [13].
5. Concluding Remarks In this paper, the CESE method is extended to solve the Navier-stokes equations. In the setting of the space-time integration of the CESE method, the inviscid flux and viscous flux are incorporated to enforce local and global flux conservation. A CFL insensitive CESE scheme is employed to provide high resolution across the computational domain, for overcoming excessive damping incurred by small CFL numbers in the original CESE method. To demonstrate the capabilities of the new approach, we calculated a shock/boundary layer interaction problem and a driven cavity flow. Numerical results show that complex physical phenomena at a wide range of Mach numbers can be predicted accurately by the CESE method.
References 1.
2.
3.
4.
Chang, S.C. and To, W. M., “A New Numerical Framework for Solving Conservation Laws – The Method of Space-Time Conservation Element and Solution Element”, NASA TM 104498, 1991. Chang, S. C., “The Method of Space-time Conservation Element and Solution Element–A New Approach for Solving the Navier-Stokes and the Euler Equations”, Journal of Computational Physics, vol. 119, 1995, pp 295-324. Chang, S. C., Wang, X.Y. and Chow, C.Y., “The Space-Time Conservation Element and Solution Element Method: A New High-Resolution and Genuinely Multidimensional Paradigm for Solving Conservation Laws”, Journal of Computational Physics, 156, 1999, pp. 89-136. Wang, X. Y. and Chang, S. C., “A 2D Non-splitting Unstructured Triangular Mesh Euler Solver Based
9 American Institute of Aerospace and Astronautics
5.
6.
7.
8.
on the Space-time Conservation Element and Solution Element Method”, Computational Fluid Dynamics Journal, vol. 8, No. 2, 1999, pp.309-325. Chang, S. C., Wang, X.Y. and To, W.M., “Application of the Space-time Conservation Element and Solution Element Method to One-Dimensional Convection-Diffusion Problems”, Journal of Computational Physics, vol. 165, 2000, pp. 189-215. Zhang, Z.C., Yu, S. T. and Chang, S. C., “A Space-Time Conservation Element and Solution Element Method for Solving the Two- and Three-Dimensional Unsteady Euler Equations Using Quadrilateral and Hexahedral Meshes”, Journal of Computational Physics. Vol. 175, 2002, pp.168-199. Zhang, M.J., Lin, S. C., Yu, S. T., Chang, S. C. and Blankson, I., “Application of the Space-Time Conservation Element and Solution Element Method to the Ideal Magnetohydrodynamics Equations”, AIAA Paper 2002-3888, 2002. Zhang, M.J., Yu, S. T., Chang, S. C. and Blankson, I., “Direct Calculation of the Ideal MHD Equations by the CESE Method without Special Treatment of Constraint of ∇ ⋅ B = 0 ”, AIAA Paper 2003-0324, 2003.
9. 10.
11.
12.
13.
Chang, S. C., Courant Number Insensitive CESE Method, AIAA Paper 2002-3890, 2002. Zhang, M.J., Yu, S. T. and Chang, S. C., “CFL Number Insensitive CESE Schemes for the Two-Dimensional Euler Equations”, AIAA Paper 2003-3842, 2003. Chang, S. C. and X. Y, Wang, “Multi-dimensional Courant Number Insensitive CESE Euler Solvers for Applications Involving Highly Non-uniform Meshes”, AIAA Paper 2003-5285. Hakkinen, R. J., Greber, I. Trilling, L. and Abarbanel, S. S., “The Interaction of An Oblique Shock Wave with A Laminar Boundary Layer”, NASA MEMO 2-18-59W, 1959. Ghia, U., Ghia, K. N. and Shin, C. T., “High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method,” Journal of Computational Physics. Vol. 48, 1982, pp.387-411.
10 American Institute of Aerospace and Astronautics
Fig. 1: Spatial computational domain with a triangular mesh. Circles (solid or hollow) are triangles’ centroids.
Fig. 2: Definition of the solution points.
Fig. 3: Grid point arrangement in the space-time domain.
Fig. 4: Definition of Solution Element SE ( j, k , n ) associated with point G ' ( j, k , n ) .
Fig. 5: Parallel translation of ∆A 'C ' E ' and ∆A*C * E * . Solution point G ' is ∆A*C * E * ’s centroid.
Fig. 6: The numerical domain of dependence in the CESE method.
11
Fig. 7: Definition of the local CFL condition for two-dimensional problems.
Fig. 8: Definition of points Q1, Q2, Q3, P1, P2 and P3 for the CFL-insensitive schemes.
Fig. 9. Geometry of shock boundary layer interaction problem
10a.
10b.
Fig. 10. Pressure contours of the shock boundary layer interaction problem with Re=296000. (a). Pressure contours obtained by CFL insensitive scheme. (b). Pressure contours obtained by original CESE scheme.
11a.
11b.
Fig. 11. Velocity vectors around shock impinging point at wall. (a). Overall view. (b). Enlarged view.
12
Fig. 12. Pressure distribution along wall.
Fig. 13. Skin friction distribution along wall.
Fig. 14. Distribution of velocity component u along the vertical centerline.
Fig. 15. Distribution of velocity component v along the horizontal centerline.
Fig. 16. Velocity vectors of the cavity driven problem.
13
Solving the Navier-Stokes Equations by the CESE Method Moujin Zhang1, S.-T. John Yu2 Mechanical Engineering Department Ohio State University, Columbus, OH 43210 Sin-Chung Chang3 NASA Glenn Research Center Cleveland, OH 44135
Abstract The CESE method is a new numerical paradigm for nonlinear hyperbolic conservation laws. In this paper, we extend the CESE method for solving viscous flows. Both viscous and inviscid fluxes are incorporated into the space-time integration of the CESE method to enforce flux conservation locally and globally. To overcome excessive numerical damping, incurred by clustered mesh for resolving the boundary layer, a CFL number insensitive scheme, recently developed by Chang, is employed. The capabilities of the present scheme are demonstrated by numerical solutions of a shock/ boundary layer interaction problem and a driven cavity flow. Without preconditioning the flow equations, the present approach can calculate flows at all speed.
1. Introduction The space-time conservation element and solution element (CESE) method is a new numerical framework for solving hyperbolic conservation laws. Contrast to modern upwind schemes, no reconstruction procedure or Riemann solver is needed and the computational logics and operational counts of the CESE method are much simpler and more efficient. Numerous one-, two-, and three-dimensional, steady and unsteady results, − Research Associate, Email: [email protected] 2 Associate Professor, Email: [email protected] 3 Senior Aerospace Engineer, Email: [email protected] 1
Copyright 2004 by the American Institute of Aeronautics and Astronatics, Inc, All right reserved
obtained by using the CESE method, have be reported in the cited references [1-8], including flows with complex shock systems, aero acoustics, flow dominated by vortices, The present paper aims at solving the Navier-stokes equations for viscous flows. In the setting of the space-time integral formation of the CESE method, the inviscid and viscous fluxes are incorporated in an equal footing manner to enforce local and global flux conservation. Moreover, to overcome the over-dissipation problem when CFL is smaller than 0.1, we also employ the newly developed CFL insensitive scheme by Chang and coworkers [9, 10, 11]. We note that numerical dissipation of the original CESE method tends to increase with decreasing local CFL number. Because the CFL numbers may vary significantly across the computational domain, usually due to non-uniform mesh, numerical solutions may become overly dissipated. To overcome this problem, Chang [9] proposed the CFL insensitive scheme for the Euler equations in one spatial dimension. The scheme was further extended for Euler equations in two spatial dimensions [10, 11]. Numerical results shows that the scheme is CFL number insensitive and could crisply resolve shocks as well as contact discontinuity. In the present paper, the new CFL insensitive methods are employed to resolve the boundary layer in solving the Navier-stokes equations. The present paper is organized as follows. In Section 2, we present the CESE method for the Navier-Stokes equations in two spatial dimensions. In section 3, we review the CFL insensitive scheme. Section 4 shows the numerical results for viscous flows. In particular, the same computer code has been used for a supersonic flow with shock/boundary layer interactions and incompressible cavity flow. We then offer concluding remarks and provide cited references.
1 American Institute of Aerospace and Astronautics
Let x1 = x , x 2 = y and x3 = t be the coordinates of a three-dimensional Eucilidean space E 3 . Equation (2.1) becomes the divergence free condition in E 3 ,
2. Numerical Methods 2.1 Space-Time Integration The dimensionless unsteady Navier-stokes equations in two-spatial-dimensions can be expressed as ∂u m ∂f m ∂g m ∂f vm ∂g vm + + − − =0, ∂t ∂x ∂y ∂x ∂y for m = 1, 2, 3, 4, where
(2.1)
(u1 , u 2 , u 3 , u 4 )
= ( ρ , ρu , ρv , e ) ,
( f 1 , f 2 , f 3 , f 4 )T
= ρu , ρu 2 + p, ρuv, (e + p )u
(g1 , g 2 , g 3 , g 4 )T
T
(2.2a)
(
)
(
)
T
, (2.2b)
= ρv, ρvu , ρv 2 + p, (e + p )v , (2.2c)
( fν 1 , fν 2 , fν 3 , fν 4 )T
T
(
= 0, τ xx , τ xy , uτ xx + vτ xy − q x
)T
(2.2d) and
(gν 1 , gν 2 , gν 3 , gν 4 )T
(
= 0, τ xy , τ yy , uτ xy + vτ yy − q y
)T
(2.2e) where ρ, p, u and v are density, pressure, velocity components in x- and y-direction. The specific total energy e is e = p (γ − 1) + ρ u 2 + v 2 2 . γ is the specific heat ratio. The stress components in Eq (2.2d-e) are as follows,
(
τ xx
)
1 4 ∂u 2 ∂v = − Re 3 ∂x 3 ∂y
τ yy =
∫∇⋅h
=
m dV
∫h
m
⋅ ds = 0 ,
(2.4)
S (V )
V
for m = 1, 2, 3, 4, where S(V) is the boundary of an arbitrary space-time region V in E 3 and ds = ndσ , where dσ and n are the area and the outward unit normal vector of a surface element on S(V).
2.2 CE and SE The two-dimensional spatial domain is divided into no-overlapped triangles. Refer to Fig. 1. Point G, the centroid of ∆BDF, is marked by a solid circle, and A, C and E, centorids of ∆FMB, ∆BJD and ∆DLF, and are marked by hollow circles. In Fig. 2, A, B, C, D, E and F form a hexagon ABCDEF. Point G , marked by a solid square is the centroid of hexagon ABCDEF, and it is the solution point of ∆BDF. Let S q1 , S q 2 and S q 3 be areas of quadrilaterals ABGF, BCDG and DEFG, and x q1 , y q1 , x q 2 , y q 2 and x q 3 , y q 3 be the spatial
(
) (
)
(
)
coordinates of their centroids. The spatial coordinates of point G are xG = yG =
1 4 ∂v 2 ∂u − Re 3 ∂y 3 ∂x
(2.3)
where h m = ( f m − fνm , g m − gνm , u m ) are the current density vectors in E 3 . By using Gauss’ divergence theorem in E 3 , we have
u m , f m , g m fνm and gνm
are flow variables, inviscid fluxes and viscous fluexs in x- and y- directions, respectively, as T
∇ ⋅ hm = 0 ,
S q1 x q1 + S q 2 x q 2 + S q 3 x q 3 S q1 + S q 2 + S q 3 S q1 y q1 + S q 2 y q 2 + S q 3 y q 3
.
(2.5)
S q1 + S q 2 + S q 3
Points A , C and E , marked by hollow squares in Fig. 2, and are the solution points of ∆FMB, ∆BJD and ∆DLF, respectively. In the space-time domain, A, B, C, D, E, F and G are at the time level n − 1 2 , and A' B ' ,
1 ∂u ∂v + Re ∂y ∂x 1 ∂T qx = 2 (γ − 1)M RePr ∂x 1 ∂T qy = 2 (γ − 1)M RePr ∂y
τ xy =
C ' , D ' , E ' , F ' and G ' are at the time level n . Points A" , B " , C " , D " , E " , F " and G " are at the time level n + 1 2 . Let j, k and n be indices for x, y and t,
where Re is the Reynolds number, Pr the Prandtl number and M the freestream Mach number.
respectively. Points G ' , A, C and G are marked by ( j, k , n ) , ( j1, k1, n − 1 2) , ( j 2, k 2, n − 1 2) and ( j 3, k 3, n − 1 2) , respectively. Shown in Fig. 3, the solution points G ' , A , C and E are placed in a staggered positions in E 3 , and their coordinates are
2 American Institute of Aerospace and Astronautics
marked by ( j, k , n ) , ( j1, k1, n − 1 2 ) , ( j 2, k 2, n − 1 2 ) and ( j 3, k 3, n − 1 2 ) . Note that, a triangle’s centroid G
'
4
( f mx )nj = ∑ ( f m,l )nj,k (ulx )nj,k ,
and the associated solution point, G ' have different spatial coordinates. In the calculation, flow variables are stored at the solution points.
( f my )nj = ∑ ( f m,l )nj,k (uly )nj,k ,
As presented in Fig. 4, the solution element SE ( j, k , n ) associated with point G ' ( j, k , n ) , is the
(g mx )nj,k
'
'
'
'
'
and SE ( j 3, k 3, n − 1 2 ) associated with points A, C and E, respectively. The surfaces of the four SEs form three CEs for point G ' . Refer to Fig. 3. Three CEs are ABGFA' B ' G ' F ' , quadrilateral cylinders ' ' ' ' ' ' ' ' CDGBC D G B and EFGDE F G D , and are referred to as CE1 ( j, k , n ) , CE 2 ( j, k , n ) and CE3 ( j, k , n ) , respectively. CE ( j, k , n ) is the union of CE1 ( j, k , n ) , CE 2 ( j, k , n ) and CE 3 ( j, k , n ) .
Inside SE ( j, k , n ) , the first-order Taylor series expansion is employed to descritize the flow variables and inviscid fluxes: u m*
( x, y , t ; j , k , n ) = ( )
u m nj ,k
+
+(
)
u mx nj ,k
u mt nj ,k
n
(
f m* (x, y, t; j, k , n ) = ( f m )nj ,k + ( f mx )nj ,k x − x j ,k + g m*
)
( ) (y − y j,k ) + ( ) (t − t ) , f my nj ,k
f mt nj ,k
(x, y , t; j, k , n ) = ( )
g m nj ,k
+(
)
g mx nj ,k
n
( )nj,k (y − y j,k ) + (g mt )nj,k (t − t n ) ,
( j, k , n ) , we let (umt )nj,k = −( f mx )nj,k − (g my )nj,k ,
(2.7)
(g my )nj,k = ∑ (g m,l )nj,k (uly )nj,k .
(2.14)
and 4
l =1
Aided by Eqs. (2.11-14), Eq. (2.9) can be recast to 4
(u mt )nj
=−
∑(f )
n m ,l j , k
l =1
(g mt )nj,k ( f mt )nj
4
for m, l = 1, 2, 3, 4. Aided by the chain rule, we have,
∑ (g ) (u ) l =1
∑(f )
=
n m ,l j , k
n ly j ,k
, (2.15)
(ult )nj,k
n m ,l j , k
l =1
4
4
∑(f ) ∑(f ) n m ,l j
4
−
n l , r j ,k
r =1
4
(u rx )nj,k
∑ ( f ) ∑ (g ) (u ) n m ,l j
4
(g mt )nj ,k
=
n l , r j ,k
r =1
∑ (g ) l =1
n m ,l j ,k
n ry j ,k
,
(2.16)
(ult )nj ,k
4
4
l =1
r =1
∑ (g m,l )nj,k ∑ ( f l ,r )nj,k (u rx )nj,k 4
∑ (g ) ∑ (g ) (u ) l =1
(2.10)
−
can be expressed as,
4
entries of Jacobian matrixes F and G, i.e.,
4
(u lx )nj,k
Aided by the chain rule and Eqs. (2.15), ( f mt )nj ,k and
−
for m = 1, 2, 3, 4. To proceed, let f m , l and g m , l be the ∂f m ∂g and g m,l = m , ∂ul ∂ul
(2.13)
(2.8)
(2.9)
(ulx )nj,k
n m ,l j ,k
,
l =1
=−
for m = 1, 2, 3, 4. At point
f m,l =
∑ (g )
l =1
(x − x j,k )
+ g my
4
=
l =1
(2.6)
(2.12)
l =1
=−
(x − x j,k )
( ) (y − y j,k ) + ( ) (t − t ) , u my nj ,k
4
'
union of four planes, i.e., the hexagon A B C D E F , the quadrilaterals B " BGG " , D " DGG " and F " FGG " , and their immediate neighborhood. Similarly, there are three SEs, i.e., SE ( j1, k1, n − 1 2 ) , SE ( j 2, k 2, n − 1 2 )
(2.11)
l =1
n m ,l j ,k
r =1
n l , r j ,k
n ry j ,k
,
(2.17)
for m = 1, 2, 3, 4. Aided by Eqs. (2.11-17), Eqs. (2.6-8) could fully specify the distribution of
u m* , f m* and g m*
inside SE ( j, k , n ) when values of (u m )nj ,k (u mx )nj ,k and
(umy )nj,k
are known.
The above expressions are all concerned with the distribution of flow variables and inviscid fluxes inside a SE. For expressing the viscous fluxes, fνm and gνm ,
3 American Institute of Aerospace and Astronautics
with u m u mx and umy , let’s express ∂u ∂x and ∂v ∂y with um , umx and umy as, ∂u 1 ∂ (ρu ) ρu ∂ρ u2 x u2u1x = − 2 = − 2 , ∂x ρ ∂x u1 ρ ∂x u1
(2.18a)
∂v 1 ∂ (ρv ) ρv ∂ρ u3 y u3u1 y = − 2 = − 2 . u1 ∂y ρ ∂y u1 ρ ∂y
(2.18b)
and EDD ' E ' belong to SE ( j3, k 3, n − 1 2) . Let S be the area of the surface. Over each area, let the outward normal vector be n, and the surface vector s = n S. The flux leaving a surface is equal to the scalar product between the vector h *m = f m* − fν*m , g m* − gν*m , u m* , evaluated at surface’s centroid, and the surface vector s. For hexahedron A' B ' C ' D ' E ' F ' in E3, its surface vector is
(
With the aid of Eq. (2.18), τ xx can be expressred as
τ xx =
uu u 1 uu u − 4 2 1x + 4 2 x + 2 3 1 y − 2 3 y . 3Re u12 u12 u1 u1
(
Similar expressions can be obtained for τ yy , τ xy , qx and q y in Eq. (2.2d-e). Moreover the viscous fluxes, fνm and gνm , can be expressred by um , umx and umy . Aided with Eq. (2.6), the distribution of flow variables inside SE ( j, k , n ) can be obtained. With flow variables, the distribution of viscous fluxes, fν*m and gν*m , can be
( )
n
n j, k
.
( )nj,k ,
With known values of (u m )nj ,k , (u mx )nj ,k and u my *
*
the distribution of flow variable u m , inviscid fluxes f m *
and g m , and viscous fluxes fν*m and gν*m , can be fully specified
(
h *m
f m*
= − approximated by
∫
SE ( j, k , n )
inside fν*m , g m*
−
gν*m , u m*
h *m ⋅ ds = 0 ,
),
)
)
s A' B 'C ' D ' E ' F ' = 0, 0, S q1 + S q 2 + S q 3 ,
(2.19)
further expressed by (u m ) j , k (u mx )nj ,k and u my
CDD 'C ' CBB 'C ' CDGB , and belong ( ) SE j 2 , k 2 , n − 1 2 . Quadrilaterals EFGD , EFF ' E ' to
.
Let
and Eq. (2.4) can be
S (V )
(
)
and the coordinates of its centroid G ' are xG , y G , t n . '
'
'
'
'
'
The flux leaving the surface A B C D E F is
(FLUX m ) A B C D E F '
'
'
'
'
'
(
)
= S q1 + S q 2 + S q 3 (u m )nj ,k . (2.22)
Let’s calculate the fluxes leaving surfaces belonging to SE ( j1, k1, n − 1 2 ) . For quadrilateral ABB ' A' , its surface vector is s ABB ' A' =
∆t ( yB − y A , x A − xB , 0 ) , 2
(2.23)
and the coordinates of its centroid are x A + x B y A + y B n −1 2 ∆t , ,t + . The flux leaving the 2 2 4 surface ABB ' A' is
(FLUX m ) ABB A '
(2.20)
(2.21)
'
=
∆t ( y B − y A )( f m )nj1−,1k12 + x A + x B − x A ( f mx )nj1−,1k12 2 2 y + yB + A − y A f my 2
( )nj1−,1k12 + ∆4t ( f mt )nj1−,1k12 − fνAB m
for m = 1, 2, 3, 4. 2.3 Time Marching for u
By imposing Eq. (2.20) over CE ( j, k , n ) , an algebraic equation can be obtained based on the global flux conservation and (u m )nj ,k can be obtained directly for
each m = 1, 2, 3, 4. To proceed, we calculate the flux leaving surfaces of CE ( j, k , n ) . Consider CE ( j, k , n ) , the hexahedral cylinder ABCDEFA' B 'C ' D ' E ' F ' referring to Fig. 3. The surfaces of CE ( j, k , n ) can be divided into four groups. Hexahedron A' B ' C ' D ' E ' F ' belongs to SE ( j, k , n ) . Quadrilaterals ABGF , ABB ' A' and AFF ' A' belong to SE ( j1, k1, n − 1 2 ) . Quadrilaterals
+
∆t (x A − x B )(g m )nj1−,1k12 + x A + x B − x A (g mx )nj1−,1k12 2 2
( )nj1−,1k12 + ∆4t (g mt )nj1−,1k12 − gνABm
y + yB + A − y A g my 2
where the viscous fluxes,
(2.24)
,
fνmAB and gνAB m , are obtained
from flow variables at the centroid of surface ABB ' A'
( )nj1−,1k12 .
and (u mx )nj1−,1k 12 and u my
For quadrilateral AFF ' A' , its surface vector is
4 American Institute of Aerospace and Astronautics
∆t ( y A − yF , xF − x A , 0 ) , 2
s AFF ' A' =
(2.25)
and the coordinates of its centroid are x A + x F y A + y F n −1 2 ∆t + . The flux leaving the , ,t 2 2 4 surface AFF ' A' is
(FLUX m ) AFF A '
=
'
∆t ( y A − y F )( f m )nj1−,1k12 + x F + x A − x A ( f mx )nj1−,1k12 2 2
( )nj1−,1k12 + ∆4t (g mt )nj1−,1k12 − gνAFm ,
y + yF + A − y A g my 2
(2.26) AF where the viscous fluxes, fνAF m and gνm , are obtained
from flow variables at the centroid of surface AFF ' A'
( )nj1−,1k12 .
and (u mx )nj1−,1k 12 and u my
The flow variables at
For quadrilateral ABGF, its surface vector is
)
s ABGF = 0, 0, − S q1 ,
(2.27)
(
)
and the coordinates of its centroid are x q1 , y q1 , t n −1 2 . The flux leaving the surface ABGF is
(
+ x q1 − x G
)( )
(
+ y q1 − y G
)( )
u my nj1−,1k12
].(2.28)
( fluxm )1n −1 2 = − S q1 [(u m )nj1−,1k12 + (x q1 − x A )(u mx )nj1−,1k12
(
+ y q1 − y A
)( )
]
∆t ( y B − y F )( f m )nj1−,1k12 + 2
∆t + [(x B + x A )( y B − y A ) + (x F + x A )( y A − y F ) 4
]
− 2 x A ( y B − y F ) ( f mx )nj1−,1k12
[
∆t 2 x F − x B2 − 2 x A (x F − x B ) (g mx )nj1−,1k 12 4
+
∆t [(x A − x B )( y A + y B ) 4
− y F )( f mt )nj1−,1k12
[
]
( )nj1−,1k12
+
.
−
(∆t )2 (x 8
F
− x B )(g mt )nj1−,1k12
[
∆t AB ( yB − y A ) fνAB m + (x A − x B )gνm 2
]
, (2.29) + ( y A − y F ) fνmAF + ( x F − x A )gνAF m
( )nj1−,1k12
where (u m )nj1−,1k12 , (u mx )nj1−,1k12 and u my
are values
stored at solution point A . Similarly, the fluxes leaving three quadrilaterals belonging to conservation over SE ( j 2, k 2, n − 1 2 ) are
[
(
)
− S q 2 (u m )nj 2−1,k 22 + x q 2 − xC (u mx )nj 2−1,k 22
( flux m )n2 −1 2 =
)( )nj 2−1,k 22 ]
(
+
∆t ( y D − y B )( f m )nj 2−1,k 22 2
+
∆t [(x D + xC )( y D − y C ) 4
]
The flux leaving three surfaces belonging to SE ( j1, k1, n − 1 2 ) is the sum of Eqs. (2.24), (2.26) and (2.28). We have
u my nj1−,1k12
+
B
+ (x B + x C )( y C − y B ) − 2 x C ( y D − y B ) ( f mx )nj 2−1,k 22
[
= − S q1 (u m )nj1-1,k21
u mx nj1−,1k12
∆t (x F − x B )(g m )nj1−,1k12 2
8
+ y q 2 − y C u my
centroid are obtained following Eq. (2.6).
(FLUX m ) ABGF
+
x +x ∆t (x F − x A )(g m )nj1−,1k12 + A F − x A (g mx )nj1−,1k12 + 2 2
(
(∆t )2 ( y
+ (x F − x A )( y A + y F ) − 2 y A (x F − x B )] g my
( )nj1−,1k12 + ∆4t ( f mt )nj1−,1k12 − fνAF m
y + yA + F − y A f my 2
+
]( )
∆t 2 + y B − y F2 − 2 y A ( y B − y F ) f my 4
n −1 2 j1,k 1
[
]( )
+
∆t 2 y D − y B2 − 2 y C ( y D − y B ) f my 4
+
(∆t )2 ( y
+
∆t (x B − x D )(g m )nj 2−1,k 22 2
+
∆t 2 x B − x D2 − 2 x C (x B − x D ) (g mx )nj 2−1,k 22 4
+
∆t [(x B − xC )( y B + yC ) 4
8
D
n −1 2 j 2,k 2
− y B )( f mt )nj 2−1,k 22
[
]
]( )nj1−,1k12
+ (xC − x D )( y D + y C ) − 2 y C (x B − x D ) g my
+
(∆t )2 (x 8
B
5 American Institute of Aerospace and Astronautics
− x D )(g mt )nj1−,1k 12
[
∆t CD ( y D − yC ) fνCD m + (xC − x D )gνm 2
−
(um )nj,k
]
CB + ( yC − y B ) fνCB m + ( x B − xC )gνm .
The fluxes leaving SE ( j 3, k 3, n − 1 2 ) are
the
three
[
( flux m )3n −1 2
)( )nj3−,1k 32 ]
2.4 Time Marching for ux and uy In this section, we illustrate the calculation of the spatial derivatives of the flow variables, i.e.,
)
(u mx )nj,k and (u my )nj,k . We note that point G '
∆t ( y F − y D )( f m )nj 3−,1k 32 + 2 ∆t [(x F + x E )( y F − y E ) 4
+ (x D + x E )( y E − y D ) − 2 x E ( y F − y D )]( f mx )nj 3−,1k 32
+
[
]( )
∆t 2 y F − y D2 − 2 y E ( y F − y D ) f my 4
(∆t )2 ( y + 8
− y D )(
F
n −1 2 j 3,k 3
+
∆t 2 x D − x F2 − 2 x E (x D − x F ) (g mx )nj 3−,1k 32 4
y E* = (3yG + 2 y E − y A − yC ) 3
[
]
( )nj3−,1k 32
− x F )(
)
]
,
,
fνEF m
,
,
gνEF m
,
fνED m
and
gνED m
fνCD m
,
, are
For flux conservation over CE ( j, k , n ) , we have '
'
'
'
'
'
(
(u )
(2.31)
AF obtained in the similar way of that for fνAF m and gνm .
(FLUX m ) A B C D E F
n −1 2 j1,k 1
(2.35)
.
(2.36)
n −1 2 mx j1,k 1 n j 2 ,k 2
+ ( flux m )1n −1 2
+ ( flux m )n2 −1 2 + ( flux m )3n −1 2 = 0 . (2.32)
' n m j 3,k 3
( )
∆t u mt 2
+
n −1 2 j1,k 1
)( )nj1−,1k12 . (2.37)
(
+ y A* − y A u my
( ) on (u )
and u m'
n j 3,k 3
, at C * and E * can
' n m j1,k 1
be obtained. Based
For each m =1, 2, 3, 4. The viscous fluxes, gνCB m
( ) )(u )
= um
( )
[
fνCB m
' n m j1,k 1
Similarly, u m'
∆t EF ( y F − y E ) fνEF − m + (x E − x F )gνm 2
gνCD m
(u )
+ x A* − x A
g mt nj 3−,1k 32
ED + ( y E − yD ) fνED m + ( xD − xE )gνm .
,
from Eq. (2.33). Flow variable at point A* is obtained by a first-order Taylor series, Eq. (2.6), as
+ (x E − x F )( y E + y F ) − 2 y E (x D − x F )] g my D
(2.34)
Flow variables at G ' , i.e., (u m )nj ,k , are calculated
∆t + [(x D − x E )( y D + y E ) 4
8
,
y A* = (3 y G + 2 y A − y C − y E ) 3
x E* = (3xG + 2x E − x A − xC ) 3
∆t (x D − x F )(g m )nj 3−,1k 32 2
(∆t )2 (x +
x A* = (3x G + 2 x A − x C − x E ) 3
y C * = (3 y G + 2 y C − y A − y E ) 3
)
is not the
centroid of ∆A 'C ' E ' unless a uniform mesh is used. As shown in Fig. 5, a triangle ∆A*C * E * , whose centroid is point G ' , is obtained by parallel moving ∆A 'C ' E ' in the spatial domain. The vertices’ coordinates of ∆A*C * E * are
xC * = (3xG + 2 xC − x A − x E ) 3
f mt nj 3−,1k 32
+
.
(2.33)
+ y q 3 − y E u my
+
S q1 + S q 2 + S q 3
of
= − S q 3 (u m )nj 3−,1k 32 + x q 3 − x E (u mx )nj 3−,1k 32
(
( fluxm )1n−1 2 + ( fluxm )n2−1 2 + ( flux m )3n−1 2
(2.30)
surfaces
(
=−
(u )
' n m j 2 ,k 2
,
and
on points A* , C * and E * , we apply central
( )nj,k
differencing to calculate (u mx )nj ,k and u my
at point
G ' , i.e.,
(u )
c n mx j ,k
(
=
1 2 S A*C *E *
)( )nj1,k1
(
y * − y * u' m E C
)( )nj 2,k 2 + (y A
+ y E * − y A* u m'
Aided by Eq. (2.22), an explicit expression for (u m )nj ,k can be obtained from above equation as
6 American Institute of Aerospace and Astronautics
*
)( )nj3,k 3 ,(2.38)
− y C * u m'
(u )
c n my j ,k
=
)( )nj1,k1
(
1 2 S A*C *E *
)( )nj 2,k 2 + (xC
(
level n = 0, 1, 2, ... , and hollow squares are at n = 1 2 , 3 2 , 5 2 , ... . According to the CESE method,
x * − x * u' m C E
+ x A* − x E * u m'
*
)( )nj3,k 3 ,(2.39)
the flow variables at solution point G ' (j, k, n) are determined by those at seven solution points A , B , C , D , E , F and G at the time level n-1. The hexagon A B C D E F is the numerical domain of dependence for the solution at G ' at time level n-1.
− x A* u m'
where S ∆A*C *E * is the area of ∆A*C * E * : S ∆A*C * E * =
(
1 x * y * + x C * y E * + x E * y A* 2 A C
)
− x A* y E * − xC * y A* − x E * y C * .
(2.40)
Similar central differencing can be applied to calculate
(u ( ) ) (u ( ) )
( )
1 n mx j ,k
(1) and u my
2 n my j ,k
for ∆A*G ' E *
n j ,k
( ) , and (u ( ) ) and (u ( ) ) (2 ) for ∆C * E *G ' , u mx 3 n mx j ,k
n j ,k
and
3 n my j ,k
for
∆A*C *G ' . Based on the re-weighting procedure in [4],
( )nj,k
we calculate (u mx )nj ,k and u my
(u mx )nj,k
( )
w = u mx
as,
( )nj,k = (umyw )nj,k ,
n j ,k
and u my
(2.41)
where
( )
w n umx j ,k
+ (θ m1θ m 3 )
α
( )
w n u my j ,k
( ) u (2 )
n mx j ,k
n j ,k
+ (θ m1θ m 2 )
α
( )
1 (1) = (θ m 2θ m 3 )α u my ω
+ (θ m1θ m 3 )
α
( ) u (2 )
n
my j ,k
( )
n mx j ,k
u (3)
, (2.42)
j ,k
+ (θ m1θ m 2 )
α
( ) u (3)
n
my j ,k
, (2.43)
θ mr
( ) ( ) 2
(r ) + u my
2
n
, j ,k
(2.44)
ω = (θ m1θ m 2 )α + (θ m 2θ m 3 )α + (θ m1θ m 3 )α .
(2.45)
The above CESE schemes are stable for CFL number < 1., while α≥0. Equations (2.33) and (2.41) form the commonly used CESE scheme, the a-α scheme. Hereafter we refer it as the original scheme.
3. CFL Insensitive Scheme 3.1 The CFL Condition The CFL number in two spatial dimensions is defined hereafter. The spatial projections of solution points are presented in Fig. 6, in which solid squares are at time
(3.1)
Let line segment G H be the distance from point G to boundary A B , we have
G H = 2 S ∆A B G
n
and
(r ) = u mx
G O , we have
G O = ∆t u 2 + v 2 , and θ 0 = arctan (v u ) .
( )
1 (1) = (θ m 2θ m 3 )α umx ω
Figure 7 shows the projection of a Mach cone on the spatial plane at t = (n-1)∆t, with point G ' being its vertex. The result is a circle with a radius of c∆t on the plane with O as the center of the circle. The boundary and interior of the circle form the domain of dependence for the flow solution at G ' . Thus, the CFL condition for the solution at G ' is defined such that if and only if the domain of dependence, i.e., the circle, lies in the interior of the hexagon A B C D E F . Let u, v and c be velocity components and sonic speed at solution point G ' at time level n-1. As shown in Fig. 7, for the velocity vector
θ1 = arctan
(x B − x A )2 + (y B − y A )2 , and
x A − xB yB − y A
.
(3.2)
As shown in Fig. 7, we choose a point P on the circle such that the line segment OP is parallel to line segment G H . Let R and S be the projection of O and P on G H . Obviously P is the closest point to the boundary A B on the circle. To keep the circle inside the hexagon, with A B as one of the boundary segments, we require that
v (1) = G S G H < 1 ,
(3.3)
where G S = ∆t c + u 2 + v 2 cos(θ 1 − θ 0 ) .
(3.4)
For the other five boundary line segments, we have similar conditions, i.e., v (2 ) for B C , v (3) for C D ,
v (4 ) for D E , v (5) for E F , and CFL condition is that,
{
v (6 ) for F A . The
}
v e = max v (1) , v (2 ) , v (3) , v (4 ) , v (5) , v (6 ) .
7 American Institute of Aerospace and Astronautics
(3.5)
ve < 1 .
(3.6)
Essentially, Eqs. (3.5-6) specify that the domain of dependence of the flow solution at G ′ must lies within its numerical domain of dependence, i.e., hexagon A B C D E F . In computation, ∆t is chosen to satisfy Eq. (3.5-6).
3.2 The CFL Insensitive Scheme According to [9-11], the CFL insensitive scheme is constructed such that (i) it would reduce to the original a scheme in the limit of CFL→0, and (ii) it would reduce to the original a-α scheme in the limit of CFL→1. As shown in Fig. 8, points Q1, Q2 and Q3 are centroids of quadrilaterals ABGF, BCDG and DEFG, points P1, P2 and P3 are defined within line segments A′Q1 , C ′Q2 and E ′Q3 as, x p1 = v e x A ′ + (1 − v e )x q1
y p1 = v e y A ′ + (1 − v e ) y q1
x p2 = v e xC ′ + (1 − v e )x q 2
y p3 = v e y E ′ + (1 − v e ) y q 3
(u )
' m p1*
( )
n −1 2 j1,k 1
= um
[(
,
(3.7)
,
(3.8)
,
(3.9)
Take Eq. (3.7) for instance, with CFL number, ve ,
1
(u )
' m p2*
( )
n −1 2 j 2,k 2
= um
[(
2
i.e., point G ′ , can be obtained by parallel translation of ∆A ' C ' E ' . The coordinates of the vertexes of ∆P1* P2* P3* are
x y
p1*
p2* p2*
(
)
(
3
p3*
)
1 − ve = v e xC * + 3xG ′ + 2 x q 2 − x q1 − x q 3 3 ,(3.11) 1 − ve = v e yC * + 3 y G ′ + 2 y q 2 − y q1 − y q 3 3
= ve y E *
(
(
∆t (u mt )nj 2−,1k 22 2
+
)
(u )
' m p3*
( )
n −1 2 j 3,k 3
= um
[(
+
(
)(
2
∆t (u mt )nj 3−,1k 32 2
)
(
)(
+ x p* − x E (u mx )nj 3−,1k 32 + y p* − y E u my 3
,
)nj 2−,1k 22 ]
3
)nj 3−,1k 32 ]
(3.15) P3* and G ′ known, the CFL insensitive scheme can be constructed. There are several versions of CFL insensitive scheme as stated in [9-11]. In the present paper, the Scheme II in [10] is applied. For completeness, the equations are provided in the following.
(u )
w n mx j ,k
≈
[
[
)
]( ) ](u ( ) ) + [1 + f (v )(s
1 n (1) 1 + f (v e )(s m1 ) j ,k u mx e
+ 1 + f (v e )(s m 2 )nj ,k
n j ,k
2 n mx j ,k
e
)
n m 3 j ,k
](u ( ) )
3 n mx j ,k
(3.16) w n my j , k
[
[
]( ) ](u ( ) ) + [1 + f (v )(s
1 (1) ≈ 1 + f (v e )(s m1 )nj ,k u my e 2 n my j , k
n
j ,k
e
)
n m3 j ,k
](u ( ) )
3 n my j , k
(3.17) where
[
]
e = 3 + f (ve )(s m1 )nj ,k + (s m 2 )nj ,k + (s m3 )nj ,k ,
(s mr )nj ,k
=
φ mr − 1, min(φ m1 , φ m 2 , φ m3 )
(3.18) (3.19)
φm1 = θ m 2θ m 3 , φm 2 = θ m1θ m 3 and φm 3 = θ m 2θ m1 , (3.20) and f (v e ) = 1 v e .
)
1 − ve 3 xG ′ + 2 x q 3 − x q1 − x q 2 3 .(3.12) 1 − ve + 3 y G ′′ + 2 y q 3 − y q1 − y q 2 3
x p* = v e x E * + y
)
1 − ve 3xG ′ + 2 x q1 − x q 2 − x q 3 = ve x A* + 3 ,(3.10) 1 − ve = ve y A* + 3 y G ′ + 2 y q1 − y q 2 − y q 3 3
(
1
)nj1−,1k12 ]
(3.14)
+ 1 + f (v e )(s m 2 )nj ,k
y
)(
+ x p* − x C (u mx )nj 2−1,k 22 + y p* − y C u my
Similar to that in section 2.4, the centroid of ∆P1* P2* P3* ,
1
(
(3.13)
(u )
(
)
+ x p* − x A (u mx )nj1−,1k12 + y p* − y A u my
changing from 0 to 1, point P1 moves from point Q1 to A ′ . Similar to the fact that the centroid of ∆A ' C ' E ' may not be G ′ , the centroid of ∆P1 P2 P3 may not be G ′ .
x p*
∆t (u mt )nj1−,1k12 2
+
With the values of flow variables at points P1* , P2* ,
y p2 = v e y C ′ + (1 − v e ) y q 2
x p3 = v e x E ′ + (1 − v e )x q 3
Flow variables at P1* , P2* and P3* are obtained by
)
8 American Institute of Aerospace and Astronautics
(3.21)
4. Numerical Results Two problems are solved by using the above CESE method for the Navier-Stokes equations.
4.1 Shock Wave Boundary Layer Interaction The shock wave boundary layer interaction problem in [12] is used as the test case. This problem is often used as a standard test for Navier-Stokes solver. As shown in Fig. 9, the computational domain is (x, y) ∈ [0, 2.4]×[0, 1.164]. The left boundary is defined as the inlet boundary where specified boundary condition is employed. The flow on the top boundary is specified to form an oblique shock, impinging on the wall. The right boundary is a supersonic outlet, where non-reflective boundary condition is used. The bottom boundary consists of a symmetric boundary and a solid wall, whose lengths are 0.8 and 1.6 respectively. For the solid wall, the no-slip boundary condition is employed. The incoming shock wave, emanating from the upper-left corner of the computational domain, impinges on the solid wall with an angle of 32.6o with respect to the wall. The flow Mach number on the left inlet boundary is 2.0. The flow condition on the top boundary is calculated based on the oblique shock condition to form the desired shock wave angle. We consider the viscous flow with Re=296,000. The computational domain is covered by 80,000 triangles. To resolve boundary layer, grid points are clustered to solid wall. The first grid point to the solid wall is 5×10-5 of the domain height. Due to the non-uniform mesh, the local CFL number varies significantly across the computational domain, the difference in CFL numbers is around 60 times. Figure 10 shows the pressure contours. There are 50 contours ranging equally from 0.01 to 0.25. The pressure contours obtained by CFL insensitive scheme are shown in Fig. 10a. Pressure contours by the original scheme are shown in Fig. 10b. The CFL insensitive scheme can crisply capture shocks, while the original scheme shows excessive numerical diffusion due to small CFL number. For viscous flow simulation, the CFL insensitive scheme plays an important role in capturing sharp shocks and detailed flow structure across the domain. Because the impinging shock is strong, boundary layer separation occurs at the shock impinging point at the wall. Figure 11 shows the velocity vectors in the recirculation zone. The change of boundary layer thickness and a boundary layer separation can be observed. Figures 12-13 show the pressure and the friction coefficient along solid wall, respectively. Curves are the calculated results and symbols are the
experimental data. The present results are in good agreement with the experimental data.
4.2 Driven Cavity Flow The driven cavity flow is also a benchmark problem for testing incompressible viscous solvers [13]. The CESE method is employed for this incompressible flow calculation. The computational domain is (x, y) ∈ [0, 1]×[0, 1.]. The top boundary is a moving solid wall, and the other three boundaries are stationary walls. We conducted computation with Re=1,000. 12,000 triangles are used to cover the computation domain. The distributions of velocity components along centerlines are plotted in Fig. 14 and 15. The velocity vector field is plotted in Fig. 16. The present results are in good agreement with the data reported in [13].
5. Concluding Remarks In this paper, the CESE method is extended to solve the Navier-stokes equations. In the setting of the space-time integration of the CESE method, the inviscid flux and viscous flux are incorporated to enforce local and global flux conservation. A CFL insensitive CESE scheme is employed to provide high resolution across the computational domain, for overcoming excessive damping incurred by small CFL numbers in the original CESE method. To demonstrate the capabilities of the new approach, we calculated a shock/boundary layer interaction problem and a driven cavity flow. Numerical results show that complex physical phenomena at a wide range of Mach numbers can be predicted accurately by the CESE method.
References 1.
2.
3.
4.
Chang, S.C. and To, W. M., “A New Numerical Framework for Solving Conservation Laws – The Method of Space-Time Conservation Element and Solution Element”, NASA TM 104498, 1991. Chang, S. C., “The Method of Space-time Conservation Element and Solution Element–A New Approach for Solving the Navier-Stokes and the Euler Equations”, Journal of Computational Physics, vol. 119, 1995, pp 295-324. Chang, S. C., Wang, X.Y. and Chow, C.Y., “The Space-Time Conservation Element and Solution Element Method: A New High-Resolution and Genuinely Multidimensional Paradigm for Solving Conservation Laws”, Journal of Computational Physics, 156, 1999, pp. 89-136. Wang, X. Y. and Chang, S. C., “A 2D Non-splitting Unstructured Triangular Mesh Euler Solver Based
9 American Institute of Aerospace and Astronautics
5.
6.
7.
8.
on the Space-time Conservation Element and Solution Element Method”, Computational Fluid Dynamics Journal, vol. 8, No. 2, 1999, pp.309-325. Chang, S. C., Wang, X.Y. and To, W.M., “Application of the Space-time Conservation Element and Solution Element Method to One-Dimensional Convection-Diffusion Problems”, Journal of Computational Physics, vol. 165, 2000, pp. 189-215. Zhang, Z.C., Yu, S. T. and Chang, S. C., “A Space-Time Conservation Element and Solution Element Method for Solving the Two- and Three-Dimensional Unsteady Euler Equations Using Quadrilateral and Hexahedral Meshes”, Journal of Computational Physics. Vol. 175, 2002, pp.168-199. Zhang, M.J., Lin, S. C., Yu, S. T., Chang, S. C. and Blankson, I., “Application of the Space-Time Conservation Element and Solution Element Method to the Ideal Magnetohydrodynamics Equations”, AIAA Paper 2002-3888, 2002. Zhang, M.J., Yu, S. T., Chang, S. C. and Blankson, I., “Direct Calculation of the Ideal MHD Equations by the CESE Method without Special Treatment of Constraint of ∇ ⋅ B = 0 ”, AIAA Paper 2003-0324, 2003.
9. 10.
11.
12.
13.
Chang, S. C., Courant Number Insensitive CESE Method, AIAA Paper 2002-3890, 2002. Zhang, M.J., Yu, S. T. and Chang, S. C., “CFL Number Insensitive CESE Schemes for the Two-Dimensional Euler Equations”, AIAA Paper 2003-3842, 2003. Chang, S. C. and X. Y, Wang, “Multi-dimensional Courant Number Insensitive CESE Euler Solvers for Applications Involving Highly Non-uniform Meshes”, AIAA Paper 2003-5285. Hakkinen, R. J., Greber, I. Trilling, L. and Abarbanel, S. S., “The Interaction of An Oblique Shock Wave with A Laminar Boundary Layer”, NASA MEMO 2-18-59W, 1959. Ghia, U., Ghia, K. N. and Shin, C. T., “High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method,” Journal of Computational Physics. Vol. 48, 1982, pp.387-411.
10 American Institute of Aerospace and Astronautics
Fig. 1: Spatial computational domain with a triangular mesh. Circles (solid or hollow) are triangles’ centroids.
Fig. 2: Definition of the solution points.
Fig. 3: Grid point arrangement in the space-time domain.
Fig. 4: Definition of Solution Element SE ( j, k , n ) associated with point G ' ( j, k , n ) .
Fig. 5: Parallel translation of ∆A 'C ' E ' and ∆A*C * E * . Solution point G ' is ∆A*C * E * ’s centroid.
Fig. 6: The numerical domain of dependence in the CESE method.
11
Fig. 7: Definition of the local CFL condition for two-dimensional problems.
Fig. 8: Definition of points Q1, Q2, Q3, P1, P2 and P3 for the CFL-insensitive schemes.
Fig. 9. Geometry of shock boundary layer interaction problem
10a.
10b.
Fig. 10. Pressure contours of the shock boundary layer interaction problem with Re=296000. (a). Pressure contours obtained by CFL insensitive scheme. (b). Pressure contours obtained by original CESE scheme.
11a.
11b.
Fig. 11. Velocity vectors around shock impinging point at wall. (a). Overall view. (b). Enlarged view.
12
Fig. 12. Pressure distribution along wall.
Fig. 13. Skin friction distribution along wall.
Fig. 14. Distribution of velocity component u along the vertical centerline.
Fig. 15. Distribution of velocity component v along the horizontal centerline.
Fig. 16. Velocity vectors of the cavity driven problem.
13
