Dec 20, 2007 - from the Faculty of Science and Technology 374. Hybrid Methods for Unsteady Fluid. Flow Problems in Complex. Geometries. JING GONG.
Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 374
Hybrid Methods for Unsteady Fluid Flow Problems in Complex Geometries JING GONG
ACTA UNIVERSITATIS UPSALIENSIS UPPSALA 2007
ISSN 1651-6214 ISBN 978-91-554-7046-3 urn:nbn:se:uu:diva-8341
! "!
#" $ $$% &$'&( ) " ) ) " "* #" + , "* - .* $$%* / 0 " ) 1 1 + 2 - * 3 *
4%* 5
*
* 67 8 9%5:9&:((:%$:4* 6 " " )) " " +"" " " ) )) " ) " ) :
" * #" " " " )) ) " ) )) " " ) ) " ) " * ; ) " ) " " " " ) ) + * " " " :: ) + +"" * ; " )) = " * 3 " ": " + " ) * #" ) ? " >@ " " A + B "
" - " " * #" " + " )) " )) ) " * 3 )) " " ) " * #" " " "
" ) 8 >7 = * #"
" " " + " " ) " ) " 8 >7 = * " " ) )) " ) "
)) ) ! " # " $ % &&'"
" ()'*+,*
" C . - $$% 6778 &(&:& 67 8 9%5:9&:((:%$:4 ' ''' :54& D" 'EE * *E F G ' ''' :54&H
To Lihong, Chutian and Churen
List of Papers
This thesis is based on the following papers, which are referred to in the text by their Roman numerals. I
II
III
IV
V VI
VII
Svärd, M., Gong, J., Nordström, J. (2006) Stable artificial dissipation operators for finite volume schemes on unstructured grids. Applied Numerical Mathematics, 56(12):1481-1490 1 Nordström, J., Gong, J. (2005) A stable and efficient hybrid method for aeroacoustic sound generation and propagation. Comptes Rendus Mecanique, l333(9):713-718 2 Nordström, J., Gong, J. (2006) A stable hybrid method for hyperbolic problems. Journal of Computational Physics, 212(2):436453 3 Nordström, J., Ham, F., Shoeybi, M., van der Weide, E., Svärd, M., Mattsson, K., Iaccarino, Gong, J. A Hybrid Method for Unsteady Fluid Flow. Technical Report 2007-020, submitted Gong, J., Nordström, J. Interface procedures for viscous problems. Technical Report 2006-019, submitted Gong, J., Nordström, J. (2007) A stable and efficient hybrid scheme for viscous problems in complex geometries. Journal of Computational Physics, 226(2):1291-1309 4 Gong, J., Nordström, J., van der Weide, E., A Hybrid Method for the Unsteady Compressible Navier-Stokes Equations. Technical Report 2007-029, submitted
In addition to papers I through VII, the following related papers are not included in the thesis. I
Gong, J., Svärd, M., Nordström, J. (2004) Artificial Dissipation for Strictly Stable Finite Volume Methods on Unstructured Meshes. Sixth World Congress on Computational Mechanics in conjunction with the Second Asian-Pacific Congress on Computational Mechanics, Beijing, China
1 Reprints
were made with permission from Elsevier. were made with permission from Elsevier. 3 Reprints were made with permission from Elsevier. 4 Reprints were made with permission from Elsevier. 2 Reprints
v
II
III
vi
Efraimsson, G., Gong, J., Svärd, M., Nordström, J. (2006) An Investigation of the performance of a high-order accurate NavierStokes code. the European Conference on Computational Fluid Dynamics ECCOMAS CFD, TU Delft, the Netherlands Svärd, M., Gong, J., Nordström, J. (2007) An accuracy evaluation of unstructured node-centered finite volume methods. Applied Numerical Mathematics, in press
Contents
1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Summation-by-parts schemes . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The artifical dissipation operators . . . . . . . . . . . . . . . . . . . . . . 1.3 Stability and conservation conditions at the interface . . . . . . . . 1.4 Multi-dimensional extensions . . . . . . . . . . . . . . . . . . . . . . . . . 2 Summary of Papers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Paper I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Paper II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Paper III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Paper IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Paper V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Paper VI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Paper VII . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Summary in Swedish . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1 1 3 4 6 9 9 10 11 13 15 16 17 21 23 25
vii
1. Introduction
Finite difference and finite volume methods are two of the most widely used numerical methods employed to solve initial boundary value problems (IBVP). High-order finite difference methods (HOFDM) are very efficient for large problems in smooth geometries and the unstructured finite volume methods (UFVM) can more readily handle complex geometries. In computational physics, the computational domain is often, for efficiency and mesh generation reasons, divided into two parts, a small region where the geometry is complex and a large region where the the geometry is simple and smooth. It is natural to develop hybrid schemes which combine the strengths of both these methods [2, 5, 6, 10, 11, 16, 17, 22]. Although some of the previous attempts have been successfully implemented, others have been known to suffer from late-time instabilities on the interfaces between the structured and unstructured domains. In this thesis, hybrid schemes designed to combine the HOFDM and the UFVM in a stable way have been developed. We mainly focus on the stability of the coupling procedure at the interface between the HOFDM and UFVM. The stability of the hybrid method is obtained by using the energy method, summation-by-parts (SBP) operators, a weak form of interface conditions and a modification of the dual mesh in the UFVM.
1.1
Summation-by-parts schemes
Consider the scalar hyperbolic equation in one dimension, ut + aux = 0
u(1 t) = g(t)
x
[1 1]
t0 (1.1)
u(x 0) = f (x)
where a 0 is a constant. Multiplying (1.1) by u and integrating yields d u2 = a(g2 (t) u2 (1 t)) dt
(1.2)
Ê
where u2 = 1 1 u2 dx. As a result, an energy estimate can be obtained and (1.1) is well-posed [14]. A finite difference approximation of (1.1) is written as ([18, 20, 24, 25, 31]) vt + aD1 v = τ P 1 e0 (v0 g)
(1.3) 1
where v = [v0 v1 vN ]T is the discrete solution vector and e0 = [1 0 0]T . τ is a penalty parameter determined by stability requirements. Note that D1 = P 1 Q is an SBP operator of first derivative since P is a norm and Q + QT = B = diag( 1 0 0 1). Based on diagonal norms, SBP operators of different orders of accuracy have been constructed [31, 20]. Multiplying (1.3) by vT P and adding the transpose, we obtain d v2P = (a + 2τ ) v0 dt
2 τ g a + 2τ
τ2 2 g a + 2τ
av2N
(1.4)
where v2P = vT Pv. If τ a2, an energy estimate is obtained. When τ = a, equation (1.4) reduces to d v2P = a(v0 dt
g)2
av2N + ag2
(1.5)
which leads to an energy estimate and stability. Note that (1.5) is similar to the continuous estimate (1.2). Next we use the edge-based node-centered finite volume scheme to discretize equation (1.1) [26, 35, 33]. In a one dimensional computational domain, the unknown variables are stored at the nodes and the boundaries of control volumes are located in the midpoints between adjacent nodes, see Figure 1.1. Ω0 0
1 2
Ω1 1
3 2
2
Ω
1 2
+12
Ω 2
3 2
1
1
Ω 1 2
Figure 1.1: The control volume in one-dimensional domain.
Integration of equation (1.1) over a control volume Ω j , leads to
Ωj
ut dx +
Ωj
aux dx = 0
j=0 1
M
(1.6)
Let Ω j = x j+1 2 x j 1 2 (x j+1 2 = (x j+1 + x j )2) and u j+1 2 (u j+1 + u j )2. The finite volume semi-discretization of (1.1) can be written in matrix form as Put + aQu = ae0 (u0 2
g)
or
ut + aD1 u = aP 1 e0 (u0
g) (1.7)
where u = [u0
uM ]T is the finite volume discrete solution vector and
1 P= 2 1 Q= 2
(x1
x0 ) (x2 x0 ) .. . 1)
(x j+1 x j
..
. (xM xM 2 ) (xM xM
1 1
1 0 .. .
1 .. . 1
..
.
0 .. .
1 .. . 1
..
.
0 1
1
1)
1
Note that D1 = P 1 Q is an SBP operator since P is a diagonal matrix with the control volumes on the diagonal and Q + QT = B = diag( 1 0 0 1). Instead of imposing the boundary condition u0 = g strongly, we retain the boundary point u0 in the scheme (1.7). At the boundary point x0 , the relation x1
x0 2
(u0 )t + a
u0
u1 2
=
a(u0
g)
(1.8)
yields the term in the right hand side of (1.7). Obviously (1.7) is a stable SBP scheme with a penalty parameter τ = a, which leads to an energy estimate of the form (1.5). Notice that the scheme (1.7) is identical to the scheme (1.3) with the second order accurate SBP operator when the computational domain is discretized with same number of equidistant grid points.
1.2
The artifical dissipation operators
Since an SBP operator is essentially a centered difference scheme, added artifical dissipation is required for shock calculations and to reduce non-physical numerical oscillations (see [12, 21, 23, 33]). For the finite difference scheme, a semi-discrete approximation of (1.1) with an artifical dissipation operator DI2p can be written, v + aP 1 Qv = DI2p v + τ P 1 e0 (v0
g)v
(1.9)
Recall that without the artifical dissipation term DI2p the scheme (1.9) has been proven stable and energy bounded in the previous section. To maintain a non-growing energy of (1.9), DI2p should be a negative semi-definite matrix. 3
In [21] a higher order accurate artificial dissipation operator DI2p is given by the form (1.10) DI2p = P 1 DTp B p D p where D p is a consistent approximation of h p (d dx p ). B p is a diagonal matrix with positive numbers of O(h) on the diagonal. Inserting (1.10) into (1.9) and applying the energy method, we obtain d v2P + 2D pv2B p = (a + 2τ )(v0 dt
τ g)2 a + 2τ
av2N
τ2 2 g a + 2τ
(1.11)
where D p v2B p = (D p v)T B p (D p v). Comparing with (1.4), the artificial dissipation operator DI2p introduces a extra term D p v2B p which damps the energy. Analogous to the previous finite difference scheme (1.9), the finite volume scheme (1.7) with an added artifical dissipation operator P 1 (DTa A p Da ) can be written ut + aD1 u =
P 1 (DTa A p Da )u
aP 1 e0 (u0
g)u
(1.12)
where A p is symmetric positive semi-definite and based on the Laplacian operator introduced in [35]. The energy method leads to d u2P + 2Dau2A p = a(u0 dt
g)2
au2N + ag2
(1.13)
where Da u2A p = (Da u)T A p (Da u). As a result, an estimate for scheme (1.12) can be obtained.
1.3 Stability and conservation conditions at the interface
We consider a computational domain x [ 1 1] which consists of two subdomains. The unknown u = [u0 u1 uM ]T is discretized by using the finite volume scheme on the left subdomain [ 1 0] while the unknown v = [v0 v1 vN ]T is discretized by using the finite difference scheme on the right subdomain [0 1], see Figure 1.2. The superscripts L and R are introduced to identify the left and right subdomains. A semi-discretization of (1.1) which combines both the finite volume and finite difference approximation with a interface at x = 0 can be written as ut + aDL1 u = (PL ) 1 DTa A p Da u + σ L (PL ) 1 eLM (uM
v0 )
vt + aDR1 v = (PL ) 1 DTp B p D p u + σ R(PL ) 1 eR0 (v0
uM )
a(PL ) 1 eL0 (u0
g)
(1.14)
where eL0 = [1 0 0]T , eLM = [0 0 1]T and eR0 = [1 0 0]T . The penalty terms σ L and σ R for the interface conditions will be determined 4
Ù0 Ù1 Ù2
Ù
Ü= 1
Ú0 Ù 1 Ù Ú1 Ü=0
Ú2
Ú
Ú 1
Ú Ü=1
Interface Figure 1.2: The finite volume and finite difference grid points in a one-dimensional domain.
below. Though the two computational subdomains coincide on the interface x = 0, the variables uM and v0 located at the interface may have different values, see Figure 1.2. To obtain a conservation condition for the equation (1.14), let be a test L = ϕ R = ϕ . Multiply the function which is smooth at the interface, that is, ϕM I 0 L T L R T R two equations in (1.14) by ( ) P and ( ) P respectively and add them together. We obtain (L )T PL ut + (R )T PR vt = ( a + σ L + aϕ0L u0 aϕNR vN
Since the relation
σ R)ϕI (uM
v0 ) L L T aϕ0 (u0 g) + (Da ) A p (Da u) + (DL1 + (D p R )T B p (D p v) + (DR1 R )PR v
a+ σL
L
(1.15) )P u L
σ R = 0, or σR = σL
a
(1.16)
cancels the interface term in equation (1.15), then (1.16) is a conservation condition on the interface. Applying the energy method to (1.14), we have d dt
u2P + v2P L
= =
R
+ 2(Da u2A p + D p v2B p ) + a(u0
au2M + 2σ L uM (uM uM v0
T
a + 2σ L
σL
σR
v0 ) + av20 + 2σ R v0 (v0
σL σR a + 2σ R
uM
v0
g)2 + av2N uM )
ag2 (1.17)
= wT Mw
A sufficient condition that leads to M 0 in (1.17) is
σL
a 2
(1.18)
where we have also used (1.16). When the conditions (1.16) and (1.18) are satisfied, we can obtain an energy estimate for (1.14). The artifical dissipation operators do not influence the interface coupling. Notice that the conservation condition (1.16) is necessary for stability. 5
1.4
Multi-dimensional extensions
The SBP schemes can be generalized to system of equations in several space dimensions. To illustrate this, we consider the hyperbolic system 1x1 0y1
ut + (Au)x + (Bu)y = 0
(1.19)
with suitable initial and boundary conditions. A and B are constant symmetric matrices with k rows and columns. As an example, the computational domain is divided into two subdomains. The UFVM is used to discretize (1.19) on subdomain [ 1 0] [0 1] with an unstructured mesh while the HOFDM is used on subdomain [0 1] [0 1] with a structured mesh, see Figure 1.3. y=1 North
West
Interface
V
U
x=−1
East
x=1
y=0 South
Figure 1.3: The hybrid mesh on the computational domain.
The semi-discrete approximation of (1.19) can be written, ut + DLx Au + DLyBu = (N L ) 1 Σ¯L (uI vt + DRx Av + DRyBv = (N R ) 1 Σ¯R (vI
vI ) + BTL uI ) + BTR
(1.20)
with DLx = (PL ) 1 QLx Ik
DRx = (PxR ) 1 QRx IyR Ik
N L = PL Ik Σ¯L = (E L )T PyL ΣL
N R = PxR PyR Ik Σ¯R = (E R )T PyR ΣR
DLy = (PL ) 1 QLy Ik
DRy = IxR (PyR ) 1 QRy Ik
on both the subdomains. DLx and DLy are the finite volume approximations of the first derivatives in x-, and y-directions, respectively. Note that DLx and DLy are SBP operators since matrix PL is a positive diagonal matrix with the control volumes on the diagonal and QLx and QLy are almost skew symmetric 6
matrices, that is, QLx + (QLx )T = Y and QLy + (QLy )T = X with non-zero elements in Y and X corresponding to the boundary points, for more details see [26, 27]. Similarly DRx and DRy are the finite difference approximations of the first derivatives in x-, and y-directions, respectively. The operators DRx and DRy are also SBP operators since matrices PxR and PyR are symmetric and positive definite and the matrices QRx and QRy are nearly skew-symmetric, i.e., R R T QRx + (QRx )T = diag( 1 0 0 1) and (Qy ) + (Qy ) = diag( 1 0 0 1). L R L R ¯ ¯ Σ , Σ are the interface penalty terms. E and E are projection matrices that map u to uI and v to vI , respectively. Moreover, The sizes of the identity matrices Ix and Iy are same as the number of mesh grid points. BT L and BT R are the penalty terms that impose the outer boundary conditions weakly. Consider a simplified case here, PyL = PyR = Py
ΣL = (ΣL )T
ΣR = (ΣR )T
(1.21)
By applying the energy method to (1.20) we obtain, d dt
u2N + u2N L
R
+ SATL + SATR = [uI vI ]T ΠT (MI Py )Π[uI vI ] (1.22)
where SAT L and SAT R collect all outer boundary terms. Π is a permutation matrix and A + 2ΣL ΣL ΣR MI = ΣL ΣR A + 2ΣR To obtain an energy estimate the penalty terms should satisfy ΣR = ΣL A ΣL 2
A
(1.23) (1.24)
It can be proven that (1.23) is the conservation condition on the interface. Note that the stability conditions (1.23) and (1.24) in two-dimensions are similar to the stability conditions (1.16) and (1.18) in one-dimension. When the HOFDM with the higher order SBP operators is used on the right domain, we must modify the control volume for the UFVM at the points on the interface to guarantee PyL = PyR , see (1.21), for more details see Paper III.
7
2. Summary of Papers
2.1
Paper I
We construct strictly stable artifical dissipation operators for the node-based unstructured finite volume schemes (UFVM) on unstructured meshes. Three different orders (first, second and fourth) accurate artifical dissipation operators have been introduced. All the artifical dissipation operators are based on the discrete Laplacian considered in [35]. An approximation to the Laplacian operator is given as 1 un ui cni dsin ∑ Ωi n¾Ni rni
(2.1)
where Ωi is the control volume of the node ri . dsin is defined as the sum of the length of the “center of mass-midpoint-center of mass” lines passing over edge, Ni denotes the set of indices of points being neighbors to ri (see Figure 2.1), cni is a weight and c ji = ci j 0 guarantees stability. Moreover, rin = ri rn .
Ω
(Í )
Figure 2.1: The grid (solid lines) and the dual grid (dashed lines) in two dimensions.
A second order artificial dissipation operator based on (2.1) can be written as 1 1 D2 = (un Ωi Ωi n∑ ¾Ni
ui )rni dsin cni Δui Ωi
(2.2) 9
A fourth order artificial dissipation operator is obtained by applying (2.2) twice, that is, 1 1 1 (2.3) D4 = D2 D2 Δ(Δui )Ω2i Ωi Ωi Ωi With a different scaling in (2.2) (rni is removed), we obtain a first order dissipation operator, ui )dsin cni Δui Ωi
1 1 D1 = (un Ωi Ωi n∑ ¾Ni
1 2
2
2
0
0
−2
−2 No dissipation 1−st order dissipation 4−th order dissipation
−4 log10(||u−3||)
log10(||u−3||)
−4
−6
−8
−8
−10
−12
−12
−14
−14
0
0.1
0.2
0.3
0.4
0.5
time
(a) on a coarse mesh
0.6
0.7
No dissipation 1−st order dissipation 4−th order dissipation
−6
−10
−16
(2.4)
−16
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
time
(b) on a fine mesh
Figure 2.2: The deviation from the steady state solution as a function of t.
The first order dissipation can be used to handle the shock waves while the fourth-order dissipation can be used to reduce non-physical high frequency oscillations. As an example to show how the dissipation operators work, we consider a steady state solution for the Burgers’ equation. In the case, a shock with speed of 2 will move out of the computational domain around t = 0 35. Without dissipation, the numerical calculation becomes unstable. However, if we instead add the artificial dissipation, the shock propagates through the boundary without problem, see Figure 2.2.
2.2
Paper II
In this paper we discuss how to combine the node based unstructured finite volume method (UFVM) with high order finite difference method (HOFDM). Both the methods employ the SBP operators and impose the boundary conditions weakly. This specific character of the methods enables us to couple these different methods together in a stable way which is important for long time calculations on realistic can be obtained, see [8, 24, 25, 26]. The coupling procedure is based on energy estimates. Numerical experiments using finite difference methods that exemplify the theoretical discus10
sions are presented. The calculations show that increased efficiency and accuracy can be obtained using hybrid methods. It has also been shown that the key to accurate results are stable interface treatments. The stability conditions for our model problem ut + aux + buy = 0 (a, b 0) are ΣL a2 and ΣR = ΣL a. To indicate how important stability is, we used ΣL = a which doesn’t satisfy the stability conditions. The calculation explodes short after T = 025 and a huge error occurs along the interface, see Figure 2.3.
Figure 2.3: Sixth order accurate SBP operators with an unstable interface treatment on a mesh with 41 41 grid point in both the left and right domains.
2.3
Paper III
In this paper we show how to combine the UFVM described and analyzed in [26] with the HOFDM considered in [8, 24, 25] into a hybrid method for hyperbolic problems. The coupling procedure is based on energy estimates. Essentially, the whole procedure can be described as a way to modify the dual grid in the finite volume method in such a way that stability can be maintained at the interface. We consider the hyperbolic system (1.19). The typical computational domain is shown in Figure 1.3. The UFVM is used on left subdomain and HOFDM is used on right subdomain. In order to obtain an energy estimate, the diagonal norms PyL for UFVM and PyR for HOFDM in y-direction should be identical at in the interface, i.e. PyL = PyR . The specific SBP operators that are based on diagonal norms are given in [20, 31]. When we use the secondorder diagonal norm PyR = h diag[12 1 ] on the right subdomain, we do 11
Interface½
Interface
¿
(a) Fourth-order SBP
(b) Sixth-order SBP
Figure 2.4: The modified control volumes for the points on the interface.
not need to change the control volume since PyL = PyR . But the standard fourthand sixth-order diagonal norms are, 17 59 48 48 13649 PyR = h diag 43200
PyR = h diag
43 49 11 48 48 12013 2711 5359 7877 43801 11 8640 4320 4320 8640 43200
(2.5)
(2.6)
respectively. In both cases we need to modify the control volume for the UFVM at the points on the interface to guarantee PyL = PyR . The old dual grid for the points at the interface consists of the lines between the center of the triangles and the midpoints of the edges. In order to match PyL and PyR , the new lines will connect the center of the triangles and the points at the interface which correspond to the PyR , see Figure 2.4. To estimate the efficiency of the hybrid method we compare that calculation with a case where UFVM is used on the whole (l + 1 unit domains large) computational domain. We estimate the efficiency for large l as NHOFDM l NHOFDM + NUFV M = ∞ (l + 1) NUFV M NUFV M
Efficiency = lim l
(2.7)
where NHOFDM NUFV M denote the number of flops for the finite difference and finite volume calculation respectively. In Figures 2.5 we can see the hybrid method are more efficient than the UFVM method. The UFVM works well on unstructured grids in complex geometries. As an example, we exclude a part of the left subdomain shaped like a NACA0012 12
30
Efficiency (%)
25 2nd order 6th order
20 15 10 5 0 3
3.5
4
4.5 log10(Nodes)
5
5.5
Figure 2.5: Estimated efficiency for right-going waves, the “source to far field” case.
airfoil with length 0 2. A calculation for waves propagating from the lower-left corner is displayed in Figure 2.6. the airfoil shaped cut-out does not introduce a significant amount of error and the numerical experiment support that the interface treatment is truly stable. 1
1
0.8
0.5
0.6
y
0
0.4 −0.5
0.2 0 −1
−0.5
0
0.5
1
−1
x
(a) solution at T=1 0.04 0.02
0.6
0
0.4
−0.02
0.2
−0.04
y
1 0.8
0 −1
−0.5
0
0.5
1
x
(b) error at T = 1
Figure 2.6: Waves propagating from lower-left corner. log(L2 error) = 1 99 on the left domain with 3172 nodes and log(L2 error) = 2 27 on the right domain with 861 nodes.
2.4
Paper IV
A hybrid scheme which combines the HOFDM and the UFVM for hyperbolic systems of equations has been developed in Paper III. In this paper we applied 13
the theoretical results to the Euler equations by using two existing efficient codes and a third coupling code. A HOFDM code (NSSUS) and a UFVM code (CDP) [15, 32] are the initial building blocks for the hybrid scheme. Both the codes use SBP operators and penalty techniques for imposing the boundary and interface conditions. A third coupling code (Chimps-lite) will administer the coupling procedure and make it possible for the two solvers to communicate in an efficient and scalable way [3], see Figure 2.7.
Figure 2.7: Schematic interface communication.
We demonstrate the capability of the hybrid method in a calculation of the flow through a two-dimensional model of a coral. In this calculation we use a 6th order accurate version of NSSUS. The geometry and the corresponding mesh can be seen in Figure 2.8. The calculation proceeds as follows. First we
(a) Global view of mesh
(b) Zoom in at interface
Figure 2.8: Geometry and grid topology of hybrid calculation around coral.
compute a steady state solution. Next, we take the steady state solution and 14
add a vortex, that becomes our initial solution. At the beginning the vortex propagates through the coral and the shape of the vortex is modified by the coral. As time passes the vortex continuous to propagate via the interface and returns to its original form again, see 2.9. The results illustrate that the hybrid
(a) t = 0 0
(b) t = 2 3
(c) t = 3 6
(d) t = 4 4
Figure 2.9: Time sequence for the vortex-coral interaction.
method is an accurate, efficient and a practically useful computational tool that can handle complex geometries as well as wave propagation phenomena.
2.5
Paper V
Stable and accurate interface treatments for the finite difference methods applied to the linear advection–diffusion equation have been studied. Three stable interface procedures — the Baumann–Oden method [4], the “borrowing” method [8] and the local discontinuous Galerkin method [29, 9, 30] have been investigated. All the treatments are based on SBP operators and the simultaneous approximation term (SAT) technique [7, 8], which lead to an energy estimate and stability. Both compact and non-compact forms satisfied the SBP property [19] was used to approximate the second derivatives. By using the com15
pact form we can obtain one order higher accuracy than for the non-compact form. Moreover, the compact form introduces less reflection and oscillation than the non-compact form. Accurate high order calculation can be achieved in both single domain and multiple domains with interfaces by using all these schemes. Table 2.1 presents the results using a non-compact form of second derivatives. The convergence rate exactly coincide with the theoretical values.
Points left+right 41 + 11 81 + 21 161 + 41 321 + 81 641 + 161
BAN
BON
LDG
L2 -Err
q
L2 -Err
q
L2 -Err
q
1 45 2 62 3 81 5 05 6 26
1 34 2 65 3 86 5 06 6 26
1 33 2 18 3 27 4 47 5 68
3 97 4 04 4 13 4 04
4 37 4 03 3 97 3 99
2 83 3 62 3 98 4 03
Table 2.1: Grid convergence rate with non-compact form by using 6th order SBP operators in two subdomains of nonuniform mesh. Where BAN: Baumann-Oden method, BON: the “borrowing” method, LDG: the local discontinuous Galerkin method.
Since the size of the spectral radius influence the efficiency of the explicit time integration [34], the spectra of these schemes in two subdomain with different value of coefficients have been compared. It is shown that the spectral radius for the schemes depend on the chosen coefficients. Specially, when the centered fluxes is used in the local discontinuous Galerkin method, the minimal spectral radius has been obtained. Moreover, the interface procedures do not increase the spectral radius if suitable penalty parameters are chosen.
2.6
Paper VI
In Paper III, it was shown how to couple the UFVM and HOFDM in a stable way for hyperbolic problems. The energy method and a modification of the dual mesh in the UFVM lead to stability. In this paper we continue the study of stable interface treatment by considering hybrid schemes for viscous problems. We also add the additional complexity of a curvilinear mesh in the HOFDM region. We consider the model problem ut + aux + buy = ε (uxx + uyy ) u(x y 0) = f (x y) ∂u = g(x y t) αu + β ∂n
16
x y¾Ω
x y ¾ Ω x y ¾ ∂ Ω
0
t
(2.8a) (2.8b)
t
0
(2.8c)
The coefficients a, b and ε are constants. In general, the coefficients α and β depend on x, y, and t . The hybrid scheme is similar to that in Paper III which combines an UFVM on unstructured meshes and an HOFDM on (stretch) structured meshes. In order to obtain a stable procedure on a curved interface for the hybrid scheme, we must modify the control volume for the UFVM according to the diagonal norm and the shape of interfaces. As an example, a curved interface (see Figure 2.10) of the form x = x(0 η ) = 0 3 sin(2πη )
y = y(0 η ) = η
0
η
1
has been introduced. The interface is discretized using 11 points and each line segment has equal length. If we use a fourth-order accurate SBP operator to approximate the first derivative. the new modified dual grid points will be located as in Figure 2.10. 1
11
0.9
10
0.8
9
Interface 0.7 8
y
0.6
7 6
0.5
5
0.4 0.3
4
0.2 3
0.1 0 −0.5
2
1
−0.4
−0.3
−0.2
−0.1
0 x
0.1
0.2
0.3
0.4
0.5
Figure 2.10: The black ’o’ represents the midpoints of the interface and the red ’x’ represents the new modified dual grid points.
We exemplify our technique by computing the steady heat distribution around a set of rods. An unstructured mesh with 10300 nodes is used on the left-bottom subdomain. Structured meshes with 81 41, 41 81, 41 41 nodes are used on the left-top, right-bottom, right-top subdomains, respectively. The steady state solution is presented in Figure 2.11. One can clearly see the advantage with this technique. The near-field around the rods is captured and the far-field part is efficiently handled.
2.7
Paper VII
In this paper we continue the work in Paper VI and consider the full NavierStokes equations. In order to apply the energy method, we use a symmetric and conservative form of the time-dependent compressible viscous Navier– 17
(a) 2D
(b) 3D
Figure 2.11: The temperature distribution at t = 3 0 around four rods.
Stokes equations in two-dimensions (see [1, 28])
ut + (A1 u)x + (A2 u)y = ε (B11 ux + B12 uy )x + (B21 ux + B22 uy )y
(2.9)
where A1 , A2 , B11 , B12 , B21 and B22 are frozen symmetric coefficient matrices. Moreover, [B11 B12 ; B21 B22 ] is a positive semi-definite matrix. We focus on the stability of the coupling procedure at the interface and treat both the finite difference–finite difference coupling and the finite difference– finite volume coupling. Similar to the methods used in Paper III, the norms PyL and PyR in y-directions on both subdomains should be identical. This means that for the finite volume–finite difference coupling, we need to modify the control volume for the finite volume scheme on the interface to guarantee the condition PyL = PyR . We exemplify the whole procedure by using finite difference methods. A stationary viscous shock problem where the middle of shock is located at the interface is calculated. The Mach number in front of the shock in the reference frame of the shock is 2 0 and the angle of the shock relative to the Cartesian frame is 15Æ . The Reynolds number Re = 50 0 is based on the Mach number of shock. Since the errors for all variables (density, velocities and energy) are very similar, only the density errors are shown in our calculations. The accuracy by using different orders accurate HOFDMs is shown in Table 2.2. The convergence rates for the second-, fourth-, sixth- and eighth- order schemes are 2, 3, 4 and 5, respectively. The results are in agreement with the theory, see [13, 36]. A typical example can be seen in Figure 2.12(a) where density isolines using the 4th order discretization are shown. The corresponding cut at y = 0 can be found in Figure 2.12(b). The distribution of density close to the interface x = 0 is very smooth, which illustrates that the interface does not introduce large reflections and oscillations. In conclusion, a stable and con18
Points (left)+(right) 3333+1733 6565+3365 129129+65129 257257+129257 513513+257513
2nd order Err
q
1 43 2 02 2 65 3 26 3 88
1 94 2 09 2 05 2 05
4th order Err
2 35 3 23 4 13 5 02
q
2 59 2 91 2 98 2 98
6th order Err 1 30 2 06 3 02 4 09 5 22
q
8th order Err 1 39 2 04 3 14 4 48 5 92
2 52 3 21 3 53 3 76
q
3 68 4 43 4 80
Table 2.2: The Convergence rates of density by using the finite difference approximations of (2.9) on two non-uniform subdomains.
0.5 rho
Y
2.55 2.40 2.25 2.10 1.95 1.80
0
1.65 1.50 1.35 1.20 1.05
-0.5
-1
-0.5
0
X
0.5
1
(a) whole computational domain
2.75
2.5
2.25
rho
2
1.75
1.5
1.25
1 -1
-0.5
0
X
0.5
1
(b) at y = 0
Figure 2.12: The density isoline with the 4th order SBP operator. A mesh of 6565 grid points used in both subdomains.
19
servative hybrid method for solving the full Navier–Stokes equation has been developed. Stable and conservative finite difference–finite difference and finite difference-finite volume coupling at interfaces have been derived. The hybrid method makes it possible to combine the efficiency of the finite difference method and the flexibility of the finite volume schemes. The numerical experiments support the theoretical conclusions and show that the interface coupling is stable and converge at the correct order.
20
3. Summary in Swedish
I denna avhandling presenteras stabila och effektiva hybridmetoder som kombinerar högre ordningens finita differens-metoder med ostrukturerade finita volym-metoder för tidsberoende initial-randvärdesproblem. Hybridmetoderna kombinerar finita differens/metodernas effktivitet med finita volym-metodernas flexibilitet. Vi utför detaljerad stabilitets-analys av hybridmetoderna, med betoning på stabiliteten vid gränssnittet mellan de strukturerade och ostrukturerade blocken. Båda metoderna använder så kallade summation-by-parts (SBP) operatorer och påtvingar rand- och gränsskikts-villkor svagt, vilket leder till en energiuppskattning och stabilitet. Vi har konstruerat och analyserat första, andra och fjärde ordningens Laplacian-baserade artificiell dissipations operatorer för ostrukturerade finita volym-metoder. Den första ordningens artificiella dissipationen kan hantera stötvågor och fjärde ordningens artificiell dissipation eliminerar ofysikaliska numeriska oscillationer effektivt. Vi visar att en stabil hybridmetod kan konstrueras genom att modifiera finita volym-metodens duala nät nära gränssnittet. Hybridmetoden appliceras på Eulers ekvationer genom att koppla ihop två separata koder för strömmningsmekaniska beräkningar. Eftersom kopplingen hanteras av en tredje kod, tillåter hybridmetoden oberoende utveckling av de separata originalkoderna. Vi visar också att hybridmetoden är noggrann, effektiv och är ett praktiskt användbart verktyg för komplexa geometrier och vågutbrednings-fenomen. Stabil och noggrann hantering av gränssnittet för den linjära advektions-diffusions ekvationen har undersökts. Korrekta högre ordningens beräkningar har genomförts med flera block och gränssnitt. Tre olika stabila metoder för gränssnittshantering — Bauman–Oden, “borrowing” och lokal diskontinuerlig Galerkin — har studerats. Analys visar endast på små skillnader mellan metoderna. En stabil och konservativ hybridmetod för paraboliska problem har utvecklats. Hybridmetoden har tillämpats på Navier–Stokes ekvationner. De numeriska resultaten stödjer de teoretiska slutsatserna och visar att kopplingen är stabil och konvergerar med rätt nogrannhetsordning för Navier–Stokes ekvationer.
21
Introduktion Finita differens och finita volym-metoder är två av de mest använda typerna av numeriska metorder för att lösa initial-randvärdesproblem. Högre ordningens finita differens-metoder (HOFDM) är mycket effektiva för att lösa stora problem i enklare geometrier. Ostrukturerade finita volyms-metoder (UFVM) är mer lämpade för mer komplexa geometrier. Av prestanda- och nätgenererings-skäl, delar man ofta upp beräkningsdomänen i två delar; en mindre med komplex geomtri och en större med enklare gemetri. Det är beräkningstekniskt intressant att utveckla hybridmetoder, som kombinerar styrkorna hos de båda metoderna och undviker svagheterna. I denna avhandling presenteras hybridmetoder som har utvecklats för att stabilt koppla samman HOFDM och UFVM. De tekniker som har använts för att designa dessa metoder bygger på SBP operatorer och svagt tillämpade rand- och gränssnitts-villkor. Vi behandlar främst stabiliteten för kopplingen vid gränssnittet mellan HOFDM och UFVM. Stabiliteten hos hybridmetoden erhålls med hjälp av energimetoden, svagt påtvingade gränssnittsvillkor och modifierat dual-nät i UFVM.
22
4. Acknowledgments
This work was carried out at the division of Scientific Computing (TDB), department of Information Technology, Uppsala University, Sweden. The work was financially supported by the Graduate School in Mathematics and Computing (FMB). First I would like to express my gratitude to my advisor Jan Nordström for his excellent guidance and support. Without his full help, this work would not have been finished. I would also like to thank my co-advisor Gunilla Efraimsson for her contribution and discussion. Many thanks to Magnus Svärd for various help and fruitful collaboration. Many thanks to all people at TDB for making a nice working environment. Especially, I want to thank Tom Smedsaas and Carina Lindgren for the secretarial and administrative help as well as Jan Kärrman and Karolina Malm Holmgren for efficient computer advice. I would like to thank Malin Ljungberg for sharing her teaching experiences. Many thanks to Kenneth Duru and Sven-Erik Ekström for their kind proofreading this thesis. Jukka Komminaho and Jarmo Rantakokko are acknowledged for assistance concerning technical and implementational aspects in making our code run on the systems of UPPMAX. I would like to thank professor Parviz Moin for giving me the opportunity to visit the Stanford University several times. I would like to express my thankfulness for the following financial supports which allowed me to travel to many interesting places: Materielanslag from FMB, Stiftelsen G S Magnusons fond and Hierta–Retzius stipendiefond from Kungliga Vetenskapsakademien, EUA4X fellowships, Knut och Alice Wallenbergs stiftelses resefond from Svenska matematikersamfunder, Ångpanneföreningens Forskningsstiftelses resefond, Rektors Wallenbergmedel and C F Liljewalchs stiftelses resetipendiefond för vetenskaplig forskning from Uppsala University. I am deeply thankful to my parents, parents-in-law, sister and brother for their endless support. I wish to thank my two children, Chutian and Churen, who make everything meaningful. Last but not least I wish to thank my wife, Lihong Xu, for her love, patience and sacrifice. This work was funded by the Graduate School in Mathematics and Computing.
23
Bibliography
[1] S. Abarbanel and D. Gottlieb. Optimal time splitting for two- and threedimensional Navier-Stokes equations with mixed derivatives. Journal of Computational Physics, 41:1–43, 1981. [2] N. Adamas and K. Shariff. High-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems. Journal of Computational Physics, 127:27–51, 1996. [3] J. J. Alonso, S. Hahn, F. Ham, M. Herrmann, G. Iaccarino, G. Kalitzin, P. LeGresley, K. Mattsson, G. Medic, P. Moin, H. Pitsch, J. Schluter, M. Svärd, E. van der Weide, D. You, and X. Wu. Chimps: A high-performance scalable module for multi-physics simulations. AIAA paper, 2006-5274, 2006. [4] C.E. Baumann and J.T. Oden. A discontinuous hp finite element method for convection-diffusion problems. Computer Methods in Applied Mechanics and Engineering, 175:331–341, 1999. [5] L. Beilina, K. Samuelsson, and K. Åhlander. Efficiency of a hybrid method for the wave equation. In Proceedings of International Conference on Finite Element Methods: Three dimensional problems. GAKUTO international Series, Mathematical Sciences and Applications, Volume 15, September 2001. [6] C. H. Blanchard, G. Gutierrez, J. A. White, and R. B. Roemer. Hybrid finite element-finite difference method for thermal analysis of blood vessels. International Journal of Hyperthermia, 16(4):341–353, 2000. [7] M. H. Carpenter, D. Gottlieb, and S. Abarbanel. Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes. Journal of Computational Physics, 111(2):220–236, 1994. [8] M.H. Carpenter, J. Nordström, and D. Gottlieb. A stable and conservative interface treatment of arbitrary spatial accuracy. Journal of Computational Physics, 148:341–365, 1999. [9] B. Cockburn and C.-W. Shu. The local discontinuous Galerkin method for timedependent convection-diffusion systems. SIAM Journal on Numerical Analysis, 35:2440–2463, 1998.
25
[10] B. Costa and W. S. Don. High order hybrid central-WENO finite difference scheme for conservation laws. Preprint, 2006. [11] M.B. Giles. UNSFLO: A numerical method the calculation of unsteady flow in turbomachinery. Report No. 205, Gas Turbine Laboratory, MIT, 1991. [12] J. Gong, M. Svärd, and J. Nordström. Artificial dissipation for strictly stable finite volume methods on unstructured meshes. In Computational Mechanics, WCCM VI in conjunction with APCOM’04, Beijing, China, September 2004. [13] B. Gustafsson. The convergence rate for difference approximation to mixed initial boundary value problems. Mathematics of Computation, 29, 1975. [14] B. Gustafsson, H.-O. Kreiss, and J. Oliger. Time Dependent Problems and Difference Methods. John Wiley & Sons, Inc., 1995. [15] F. Ham, K. Mattsson, and G. Iaccarino. Accurate and stable finite volume operators for unstructured flow solvers. Annual Research Briefs, CTR, Stanford University, 2006. [16] H. Hefazi, V. Chin, and L.T. Chen. Two-dimensional zonal Euler method using composite structured and unstructured meshes. Journal of Aircraft, 31:651– 658, 1994. [17] D. Hill and D. Pullin. Hybrid tuned center-difference-WENO method for large eddy simulations in the presence of strong shocks. Journal of Computational Physics, 194:435–450, 2004. [18] H.-O. Kreiss and G. Scherer. Finite element and finite difference methods for
hyperbolic partial differential equations, in: C. De Boor (Ed.), Mathematical Aspects of Finite Elements in Partial Differential Equation. Academic Press, New York, 1974. [19] K. Mattsson. Boundary procedures for summation-by-parts operators. Journal of Scientific Computing, 18(1):133–153, 2003. [20] K. Mattsson and J. Nordström. Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics, 199:503–540, 2004. [21] K. Mattsson, M. Svärd, and J. Nordström. Stable and accurate artificial dissipation. Journal of Scientific Computing, 21(1):57–79, 2004. [22] K. Nakahashi and S. Obayashi. FDM-FVM zonal approach for viscous flow computations over multiple bodies. AIAA Paper 87-0604, 1987.
26
[23] J. Nordström. Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation. Journal of Scientific Computing, 29(3):375–404, 2006. [24] J. Nordström and M. H. Carpenter. Boundary and interface conditions for high order finite difference methods applied to the Euler and Navier-Stokes equations. Journal of Computational Physics, 148:621–645, 1999. [25] J. Nordström and M. H. Carpenter. High-order finite difference methods, multidimensional linear problems and curvilinear coordinates. Journal of Computational Physics, 173:149–174, 2001. [26] J. Nordström, K. Forsberg, C. Adamsson, and P. Eliasson. Finite volume methods, unstructured meshes and strict stability. Applied Numerical Mathematics, 45:453–473, 2003. [27] J. Nordström and J. Gong. A stable hybrid method for hyperbolic problems. Journal of Computational Physics, 212:436–453, 2006. [28] J. Nordström and M. Svärd. Well-posed boundary conditions for the NavierStokes equations. SIAM Journal on Numerical Analysis, 43(3):1231–1255, 2005. [29] W.H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. [30] C.-W. Shu. Different formulations of the discontinuous galerkin method for the viscous terms. In Z.-C. Shi, M. Mu, W. Xue, and J. Zou, editors, Advances in Scientific Computing. Science Press, 2001. [31] B. Strand. Summation by parts for finite difference approximation for d/dx. Journal of Computational Physics, 110(1):47–67, 1994. [32] M. Svärd and E. Van der Weide. Stable and high-order accurate finite difference schemes on singular grids. Annual Research Briefs, CTR, Stanford University, 2006. [33] M. Svärd, J. Gong, and J. Nordström. Stable artificial dissipation operators for finite volume schemes on unstructured grids. Applied Numerical Mathematics, 56(12):1481–1490, 2006. [34] M. Svärd, K. Mattsson, and J. Nordström. Steady-stable computations using summation-by-parts operators. Journal of Scientific Computing, 24(1):79– 95, 2005. [35] M. Svärd and J. Nordström. Stability of finite volume approximations for the Laplacian operator on quadrilateral and triangular grids. Applied Numerical Mathematics, 51:101–125, 2004.
27
[36] M. Svärd and J. Nordström. On the order of accuracy for difference approximations of initial-boundary value problems. Journal of Computational Physics, 218(1):333–352, 2006.
28
Acta Universitatis Upsaliensis Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 374 Editor: The Dean of the Faculty of Science and Technology A doctoral dissertation from the Faculty of Science and Technology, Uppsala University, is usually a summary of a number of papers. A few copies of the complete dissertation are kept at major Swedish research libraries, while the summary alone is distributed internationally through the series Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology. (Prior to January, 2005, the series was published under the title “Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology”.)
Distribution: publications.uu.se urn:nbn:se:uu:diva-8341
ACTA UNIVERSITATIS UPSALIENSIS UPPSALA 2007
