- Email: [email protected]
Discontinuous Galerkin solution of compressible flow in time-dependent domains
Available online at www.sciencedirect.com
Mathematics and Computers in Simulation 80 (2010) 1612–1623
Discontinuous Galerkin solution of compressible flow in time-dependent domains M. Feistauer ∗ , V. Kuˇcera, J. Prokopová Charles University Prague, Faculty of Mathematics and Physics, Sokolovská 83, 186 75 Praha, Czech Republic Received 30 September 2008; received in revised form 19 January 2009; accepted 29 January 2009 Available online 10 February 2009
Abstract This work is concerned with the simulation of inviscid compressible flow in time-dependent domains. We present an arbitrary Lagrangian–Eulerian (ALE) formulation of the Euler equations describing compressible flow, discretize them in space by the discontinous Galerkin method and introduce a semi-implicit linearized time stepping for the numerical solution of the complete problem. Special attention is paid to the treatment of boundary conditions and the limiting procedure avoiding the Gibbs phenomenon in the vicinity of discontinuities. The presented computational results show the applicability of the developed method. © 2009 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Compressible flow; Time-dependent domain; Fluid–structure interaction; ALE method; Discontinuous Galerkin method; Shock capturing; Flow in channels with moving walls
1. Introduction The interaction of fluid flow with vibrating bodies plays a significant role in many areas of science and technology. We mention, for example, development of airplanes (vibrations of wings) or turbines (blade vibrations), some problems from civil engineering (interaction of wind with constructions as bridges, TV towers or cooling towers of power stations), car industry (vibration of various elements of a carosery), but also medicine (hemodynamics or flow in vocal folds). In a number of these examples the moving medium is a gas, i.e. compressible flow. For low Mach number flows incompressible models are used (as e.g. in [13]), but in some cases compressibility plays an important role. The solution of fluid–structure interaction requires the coupling of the solution of equations describing the fluid flow with equations describing the structural behaviour. Due to the deformation and/or vibrations of structures, the computational domain is time dependent. There exist several techniques of the solution of incompressible flow in time-dependent domains. See, e.g. [13] and references therein. The numerical simulation of compressible flow is much more difficult, particularly in time-dependent domains. It is necessary to overcome difficulties caused by nonlinear convection dominating over diffusion, which leads to boundary layers and wakes for large Reynolds numbers and shock waves and contact discontinuities for high Mach numbers and instabilities caused by acoustic effects for low Mach numbers.
∗
Corresponding author. E-mail addresses: [email protected] (M. Feistauer), [email protected] (V. Kuˇcera), [email protected] (J. Prokopová).
0378-4754/$36.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2009.01.020
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1613
It appears that a suitable numerical method for the solution of compressible flow is the discontinuous Galerkin finite element method (DGFEM). It employs piecewise polynomial approximations without any requirement on the continuity on interfaces between neighbouring elements. The DGFEM was used for the numerical simulation of the compressible Euler equations, for example, by Bassi and Rebay in [1], where the space DG discretization was combined with explicit Runge-Kutta time discretization. In [2] Baumann and Oden describe an hp version of the space DG discretization with explicit time stepping to compressible flow. Van der Vegt and van der Ven apply space-time discontinuous Galerkin method to the solution of the Euler equations in [14], where the discrete problem is solved with the aid of a multigrid accelerated pseudo-time-integration. The papers [6,8] are concerned with a semi-implicit DGFEM technique for the solution of inviscid compressible flow, which is unconditionally stable. In [10], this method was extended so that the resulting scheme is robust with respect to the magnitude of the Mach number. For some DG techniques see [3]. The goal of our research is the numerical simulation of interaction of compressible flow with structures as, e.g. flow-induced airfoil vibrations or the flow past an elastic wall vibrating due to the influence of a moving gas. The present paper represents the first step to the solution of this problem. We are concerned with the generalization of the method from [10] to the solution of inviscid compressible flow in time-dependent domains. The main ingredients of the method is the discontinuous Galerkin space semidiscretization of the Euler equations written in the arbitrary Lagrangian–Eulerian (ALE) form, semi-implicit time discretization, suitable treatment of boundary conditions so that they are transparent for acoustic waves at the inlet and outlet and the shock capturing avoiding Gibbs phenomenon at discontinuities. Numerical experiments prove the applicability of the method. 2. Continuous problem We are concerned with inviscid compressible flow in a bounded domain t ⊂ R2 depending on time t ∈ [0, T ]. Because of simplicity we consider two-dimensional flow, but the method can be applied to 3D flow as well. We assume that the boundary of t consists of three disjoint parts I , O , Wt : ∂t = I ∪ O ∪ Wt , where I and O represent the time-independent inlet and outlet, respectively, and Wt represents moving impermeable walls. The flow is described by the Euler equations written in the conservative form ∂w ∂f s (w) + =0 ∂t ∂xs 2
in t ,
t ∈ (0, T ),
(1)
s=1
where w = (w1 , . . . , w4 )T = (ρ, ρv1 , ρv2 , E)T is the so-called state vector and f s (w) = (ρvs , ρvs v1 + δs1 p, ρvs v2 + δs2 p, (E + p) vs )T i = 1, 2, are the inviscid (Euler) fluxes of the quantity w in the directions xs , s = 1, 2. We use the following notation— ρ: density, p: pressure, E: total energy, v = (v1 , v2 ): velocity vector, γ > 1: Poisson adiabatic constant, x = (x1 , x2 ): points in t and X = (X1 , X2 ): points of the reference configuration 0 . System (1) is completed by the thermodynamical relation 1 2 p = (γ − 1) E − ρ|v| , (2) 2 boundary conditions, which are treated in Section 5, and the initial condition w(x, 0) = w0 (x),
x ∈ 0 .
(3)
We define the flux of the quantity w in the direction n = (n1 , n2 ) ∈ R2 , by F (w, n) =
2
f s (w)ns
(4)
s=1
and its Jacobi matrix DF (w, n) As (w)ns , = Dw 2
P(w, n) =
s=1
(5)
1614
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
where Df s (w) , s = 1, 2, (6) Dw are the Jacobi matrices of the mappings f s . It is possible to show that f s , s = 1, 2, are homogeneous mappings of order one, which implies that As (w) =
f s (w) = As (w)w,
s = 1, 2.
(7)
The matrix P is diagonalizable (cf. [9], Section 3.1.5), i.e. there exist matrices T = T(w, n) and = diag(λ1 , . . . , λ4 ) such that P = TT−1 .
(8)
3. ALE formulation The time dependence of the domain is taken into account with the aid of a regular one-to-one ALE mapping (cf. [12]) ¯0 → ¯ t , i.e. X ∈ ¯ 0 → x = x(X, t) = At (X) ∈ ¯ t. At :
(9)
We define the ALE velocity: ∂ At (X), ∂t z(x, t) = z˜ (A−1 t (x), t),
z˜ (X, t) =
t ∈ [0, T ], X ∈ 0 ,
(10)
t ∈ [0, T ], x ∈ t
and the ALE derivative of a function f = f (x, t) defined for x ∈ t and t ∈ [0, T ]: DA ∂f˜ f (x, t) = (X, t), Dt ∂t
(11)
where f˜ (X, t) = f (At (X), t),
X ∈ 0 , x = At (X).
(12)
By the chain rule, DA f ∂f ∂f = + z · ∇f = + div (zf ) − f div z. Dt ∂t ∂t This leads us to the ALE form of the Euler equations: DA w ∂gs (w) + + w div z = 0, Dt ∂xs
(13)
2
(14)
s=1
where gs , s = 1, 2, are modified inviscid fluxes gs (w):=f s (w) − zs w.
(15)
4. Discretization of the problem in the time-dependent domain In this section, we shall describe the discretization of the initial-boundary value problem for the Euler equations written in the ALE form (14). 4.1. Space semidiscretization The space semidiscretization will be carried out by the discontinuous Galerkin finite element method (DGFEM) based on piecewise polynomial approximations constructed over a triangular mesh without any requirement on the continuity on interfaces between neighbouring elements.
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1615
¯ ht consisting of closed Let ht be a polygonal approximation of t at time t. By Tht we denote a partition of triangles Ki ∈ Tht , i ∈ It , with disjoint interiors. (It ⊂ Z+ = {0, 1, 2, . . .} is a suitable index set.) The mesh Tht can be nonconforming with hanging nodes. By ij we denote a common edge between two neighbouring elements Ki and Kj . We set st (i) = {j ∈ It ; Kj is a neighbour of Ki }. The boundary ∂ht is formed by a finite number of faces of elements Ki adjacent to ∂ht . We denote all these boundary faces by Sj , where j ∈ Ibt ⊂ Z− = {−1, −2, . . .}. Now we set γt (i) = {j ∈ Ibt ; Sj is a face of Ki ∈ Tht } and ij = Sj for Ki ∈ Tht such that Sj ⊂ ∂Ki , j ∈ Ibt . For Ki not containing any boundary face Sj we set γt (i) = ∅. Obviously, st (i) ∩ γt (i) = ∅ for all i ∈ It . If we write St (i) = st (i) ∪ γt (i), we have ∂Ki = ij , ∂Ki ∩ ∂ht = ij . (16) j ∈ St (i)
j ∈ γt (i)
The symbol nij = ((nij )1 , (nij )2 ) will denote the unit outer normal to ∂Ki on the face ij . The approximate solution will be sought at each time instant t as an element of the finite-dimensional space S ht = S r,−1 (ht , Tht ) = {ϕh ; ϕh |K ∈ P r (K)
∀K ∈ Tht }4 ,
(17)
where r ≥ 0 is an integer and P r (K) denotes the space of all polynomials on K of degree ≤ r. Function ϕ ∈ S ht are in general discontinuous on interfaces ij . By ϕ|ij and ϕ|ji we denote the values of ϕ on ij considered from the interior and the exterior of Ki , respectively. The symbols 1 (ϕ|ij + ϕ|ji ), 2
ϕij =
[ϕ]ij = ϕ|ij − ϕ|ji
(18)
denote the avarege and jump of a function ϕ on ij = ij . In order to derive the discrete problem, we assume that w is a sufficiently regular solution of system (14), multiply (14) by a test function ϕ ∈ S ht , integrate over any element Ki , i ∈ It , apply Green’s theorem and sum over all i ∈ It . We get the relation Ki ∈ Tht
Ki
2 2 DA w(t) ∂ϕ gs (w(t)) · dx − gs (w(t))(nij )s · ϕ dS · ϕ dx = Dt ∂xs Ki ∈ Tht Ki s=1 Ki ∈ Tht j ∈ St (i) ij s=1 − div z(w · ϕ) dx. (19) Ki ∈ Tht
Ki
Now the exact solution w(t) is approximated by an element wh (t) ∈ S ht and the fluxes through the faces ij are approximated with the aid of a numerical flux H = H(u, w, n). It means that
2
ij s=1
gs (w(t)) nij
s
· ϕ dS ≈ ij
H g (wh (t)|ij , wh (t)|ji , nij ) · ϕ dS.
(20)
Here H g is an analogy to the Vijayasundaram numerical flux consistent with the fluxes gs , s = 1, 2. Taking into account (15), (5) and (6), we have Dgs (w) Df s (w) = − zs I = As − zs I, Dw Dw
(21)
and can write ˜ P(w, n) =
2 Dgs (w) s=1
Dw
ns =
2
(As ns − zs ns I) = P(w, n) − (z · n)I.
(22)
s=1
This and (8) imply that ˜ −1 , P˜ = TT
˜ = diag(λ˜ 1 , . . . , λ˜ 4 ),
λ˜ i = λi − z · n,
i = 1, . . . , 4.
(23)
1616
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
Now we define the “positive” and “negative” parts of the matrix P˜ as ˜ ± = diag(λ˜ ± , . . . , λ˜ ± ), 1 4
P˜ ± = T± T−1 ,
(24)
where λ+ = max(λ, 0), λ− = min(λ, 0). Using the above concepts, we introduce the modified Vijayasundaram numerical flux (cf. [15] or [9]) as wL + w R wL + w R H g (wL , wR , n) = P˜ + (25) , n wL + P˜ − , n wR . 2 2 If we define the forms (ψh , ϕh )h = ψh · ϕh dx,
(26)
ht
b˜ h (wh , ϕh ) = −
Ki ∈ Tht
dh (wh , ϕh ) = −
2
Ki s=1
Ki ∈ Tht
Ki
gs (wh ) ·
∂ϕh dx + ∂xs
Ki ∈ Tht j ∈ St (i) ij
H g (wh |ij , wh |ji , nij ) · ϕh dS,
divz(wh · ϕh ) dx,
(27)
(28)
we arrive at the semi-implicit discrete problem. We define an approximate solution as a function wh = wh (t) satisfying the conditions (a) (b)
wh (t) ∈ S ht , ∀t ∈ [0, T ],
DA wh (t) + b˜ h (wh (t), ϕh ) − dh (wh (t), ϕh ) = 0 , ϕh Dt
∀ϕh ∈ S ht ,
∀t ∈ (0, T ),
(29)
h
(c)
wh (0) = h w0 ,
where h w0 is the L2 (0 )-projection of w0 from the initial condition (3) on the space S h0 . 4.2. Time discretization The discrete problem (29) is equivalent to a large system of ordinary differential equations. For solving this system one can apply several numerical schemes like Runge-Kutta methods. However, they are conditionally stable and the time step is strongly limited by the CFL-stability condition. For low Mach number flows these methods fail. Therefore, we apply here a semi-implicit scheme, which is a generalization of the technique proposed in [10]. We introduce the partition 0 = t0 < t1 < . . . of the time interval [0, T ] and set τj = tj+1 − tj . The function wh (·, tj ) j will be approximated by wj , defined in tj . Let us assume that the approximate solution wh has already been computed for j = 0, . . . , k. We are interested in the computation of the approximate solution wk+1 at time instant tk+1 . If we set h ˆ h (x) = wh (Atj (A−1 w tk+1 )(x)), j
j
x ∈ htk+1 ,
(30)
then, on the basis of (11), we can approximate the ALE derivative using the first order backward difference:
ˆ kh −w wk+1 DA wh (t) h , ϕh ≈ , ϕh . Dt τk h t=tk+1
h
(31)
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1617
Further, on the basis of (7) and the definition of the modified Vijayasundaram numerical flux we define a partial linearization bh of the form b˜ h : 2 ∂ϕ (x) ˆ kh , wk+1 ˆ k (x)) − zs (x)I)wk+1 (x)) · h , ϕ bh (w ) = − ((As (w dx h h ∂xs Ki Ki ∈ Thtk+1
+
s=1
Ki ∈ Thtk+1 j ∈ Stk+1 (i) ij
˜ − ˆ k , nij )wk+1 | ) · ϕh dS. ˆ kh ij , nij )wk+1 (P˜ + (w ji h ij h |ij + P (w h (32)
The term dh linear with respect to wh will be treated implicitly. These considerations lead us to the following semi-implicit scheme: For k = 0.1, . . . find wk+1 such that h (a) (b)
∈ S htk+1 , wk+1 h
k+1 ˆ kh wh − w k+1 ˆ kh , wk+1 , ϕh + bh (w h , ϕh ) − dh (wh , ϕh ) = 0 τk
∀ϕh ∈ S htk+1 ,
(33)
h
(c)
w0h = h w0 .
This is a first order accurate scheme in time. In the solution of nonstationary flows, it is necessary to apply a scheme, which is sufficiently accurate in space as well as in time. One possibility is to apply the following two step second order time discretization based on the backward difference formula (BDF): DA wh (x, t) 2τk + τk−1 k+1 τk + τk−1 k τk ˆ h (x) + ˆ k−1 (x), x ∈ htk+1 , ≈ w (x) − w w Dt τk (τk + τk−1 ) h τk τk−1 τk−1 (τk + τk−1 ) h t=tk+1 (34) for the approximation of the ALE derivative of the solution at tk+1 and the second order approximation ˜ k+1 w = h
τk + τk−1 k τk k−1 ˆh − ˆ w w , τk−1 τk−1 h
(35)
ˆ kh in the form bh . As a result we get the of wh (tk+1 ) obtained with the aid of extrapolation, which replaces the state w ∈ S htk+1 such that following two-step second-order scheme: For each k ≥ 1 find wk+1 h 2τk + τk−1 τk + τk−1 k τk k+1 ˆ k−1 , ϕh )h ˜ k+1 ˆ h , ϕ h )h − (wk+1 (w (w h , ϕh )h + bh (w h , wh , ϕh ) = τk (τk + τk−1 ) τk τk−1 τk−1 (τk + τk−1 ) h × ∀ϕh ∈ S htk+1 ,
k = 0, 1, . . . ,
(36)
where w0h = h w0 , w1h obtained by the Runge-Kutta method. In order to guarantee the stability of the scheme, we use the CFL condition 1 τk maxKi ∈ Thtk (maxj ∈ S(i) |ij |λmax ) ≤ CFL, (37) ˆ kh |ij ,nij ) P˜ (w |Ki | where |Ki | denotes the area of Ki , |ij | is the length of the edge ij , CFL is a given constant and λmax ˆ k| P˜ (w
h ij ,nij )
is the
˜ w ˆ kh |ij , nij ), here the maximum is taken over ij . This condition is a heuristic maximal eigenvalue of the matrix P( extension of a similar condition applied in the finite volume method (see, e.g. [9]). Numerical experiments show that the presented method is practically unconditionally stable and the CFL number can be chosen large—cf. [10]. 5. Boundary conditions If ij ⊂ ∂ht , i.e. j ∈ γt (i), it is necessary to specify the boundary state w|ji appearing in the numerical flux H g in the definition of the inviscid form bh . Similarly as in [10], the choice of the boundary conditions is governed by the ˜ sign of the eigenvalues λ˜ i , i = 1, . . . , 4, of the matrix P(w, n).
1618
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
On the inlet and outlet, which are assumed fixed, we proceed in the same way as in [10], Section 4, where we prescribe the state w|ji in such a way that the locally linearized Euler equations are well posed. Using the rotational invariance, we transform the Euler equations to the coordinates x˜ 1 , parallel with the normal direction n = (n1 , n2 ) = nij to the boundary, and x˜ 2 , tangential to the boundary, neglect the derivative with respect to x˜ 2 and linearize the system around the state qij = Q(n)w|ij , where ⎛ ⎞ 1, 0, 0, 0 ⎜ 0, n , n , 0 ⎟ 1 2 ⎜ ⎟ Q(n) = ⎜ (38) ⎟ ⎝ 0, −n2 , n1 , 0 ⎠ 0,
0,
0,
1
is the rotational matrix. Then we obtain the linear system ∂q ∂q = 0, + A1 (qij ) ∂t ∂˜x1
(39)
for the transformed vector-valued function q = Q(n)w, considered in the set (−∞, 0) × (0, ∞) and equipped with the initial and boundary conditions q(˜x1 , 0) = qij , x˜ 1 < 0,
and q(0, t) = qji , t > 0.
(40)
The goal is to choose qji in such a way that this initial-boundary value problem is well posed, i.e. has a unique solution. The method of characteristics leads to the following process: Let us put q∗ji = Q(n)w∗ji , where w∗ji is a given boundary state at the inlet or outlet. We calculate the eigenvectors rs corresponding to the eigenvalues λs , s = 1, . . . , 4, of the matrix A1 (qij ), arrange them as columns in the matrix T and calculate T−1 . Now we set α = T−1 qij ,
β = T−1 q∗ji .
and define the state qji by the relations 4 αs , λs ≥ 0, qji := γs rs , γs = βs , λs < 0. s=1
(41)
(42)
Finally, the sought boundary state w|ji is defined as w|ji = wji = Q−1 (n)qji .
(43)
On the impermeable moving wall we prescribe the normal component of the velocity v · n = z · n,
(44)
˜ n) vanish, where n is the unit outer normal to Wt and z is the wall velocity. This means that two eigenvalues of P(w, one is positive and one is negative. Then, in analogy to [9], Section 3.3.6, we should prescribe one quantity, namely v · n, and extrapolate three quantities—tangential velocity, density and pressure. However, here we define the numerical flux on Wt as the physical flux through the boundary with the assumption (44) taken into account. Thus, on Wt we write 2
gs (w)ns = (v · n − z · n)w + p (0, n1 , n2 , v · n)T = p (0, n1 , n2 , z · n)T =: H g .
(45)
s=1
As was shown in [1] or [9], Section 4.6.8, the piecewise linear approximation of a curved boundary can lead to a spurious nonphysical entropy production near vertices. In order to avoid this drawback, isoparametric curved elements should be used. In our code we follow [9], Section 4.6.8, where the reference triangle is transformed by a bilinear mapping onto the approximation of a curved triangle adjacent to the boundary ∂.
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1619
Fig. 1. Triangulations at time instants t = 0.0, 6.1392 and 10.0592.
6. Limiting procedure at discontinuities For high speed flows with shock waves and contact discontinuities it is necessary to avoid the Gibbs phenomenon manifested by spurious overshoots and undershoots in computed quantities near discontinuities. In order to avoid the Gibbs phenomenon, we apply the limiting procedure from [10] based on the discontinuity indicator 2 [ρˆ hk ] dS gk (i) = , Ki ∈ Thtk+1 , (46) 3/4 ∂Ki hKi |Ki | ˆ kh . Then we define the discrete introduced in [7]. The density ρˆ hk represents the first component of the state vector w discontinuity indicator Gk (i) = 0,
if gk (i) < 1,
Gk (i) = 1,
if gk (i) ≥ 1,
and the artificial viscosity forms k ˆ kh , wk+1 βh (w , ϕ ) = ν h G (i) 1 Ki h h i ∈ Itk+1
and ˆ kh , wk+1 Jh (w h , ϕh ) = ν2
i ∈ I j ∈ stk+1 (i)
Ki
Ki ∈ Thtk+1 ,
∇wk+1 · ∇ϕh dx h
1 k (G (i) + Gk (j)) 2
(47)
(48)
ij
[wk+1 h ] · [ϕh ] dS,
(49)
with ν1 , ν2 = O(1). Then the resulting scheme obtained by limiting of (33)(b) has the following form: we seek wk+1 ∈ S htk+1 such that h
ˆ kh −w wk+1 h ˆ kh , wk+1 ˆ kh , wk+1 ˆ kh , wk+1 , ϕh + bh (w h , ϕh ) + βh (w h , ϕh ) + Jh (w h , ϕh ) = 0, ∀ ϕh ∈ S htk+1 , τk h
k = 0, 1, . . . , (Similarly we obtain a stabilized version of scheme (36).)
1620
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
Fig. 2. Velocity (left) and pressure (right) isolines at time instants t = 0.3992, 2.2192, 3.1992, 4.1792, 6.1392 and 7.1192.
In practical computations, the integrals appearing in the definition of the approximated solution are evaluated with the aid of quadrature formulae. Linear algebraic systems are solved either with the aid of the direct solver UMFPACK [4] or by the GMRES method with block diagonal preconditioning. 7. Numerical experiments We consider inviscid compressible flow in the rectangular channel with the initial shape 0 = [−2, 2] × [0, 1], where the lower wall of the channel is moving. The ALE mapping is equal to identity in the sets [−2, −1] × [0, 1] and
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1621
Fig. 3. Velocity (left) and pressure (right) isolines at time instants t = 8.0992, 10.0592, 12.0192, 13.9792, 14.9592 and 15.9392.
[1, 2] × [0, 1]. In the rest of the channel we construct the ALE mapping so that the lower wall is represented at time t by the graph of the function 0.45 sin(0.4t) (cos(πX1 ) + 1),
X1 ∈ (−1, 1).
(50)
This movement is interpolated to the rest of the domain resulting in the ALE mapping At . The computation was carried out with the dimensionless form of the Euler equations, using the dimensionless conservation variables, on a triangulation with 1160 elements. Fig. 1 shows the triangulation at time t = 0.0, 6.1392 and 10.0592. The triangular mesh was constructed by the technique from [5]. The inlet Mach number Minlet = 0.067,
1622
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
the dimensionless inlet density ρinlet = 1.0. At the outlet, the dimensionless pressure poutlet = 159.12. In the artificial viscosity forms (48) and (49) the values ν1 = ν2 = 0.2 were used. Figs. 2 and 3 show velocity and pressure isolines at different time instants t = 0.3992, 2.2192, 3.1992, 4.1792, 6.1392, 7.1192, 8.0992, 10.0592, 12.0192, 13.9792, 14.9592 and 15.9392 during one period. In the solution we can observe a vortex formation, when the lower wall starts to descend. This vortex is convected through the domain. Moreover, we see that a contact discontinuity is developed, when the channel becomes narrow. The contact discontinuity is characterized by the discontinuity in the velocity, whereas the pressure remains continuous. In Fig. 2, one can see some “artefacts” in pressure isolines near the moving wall at time t = 6.1392, which might be caused by the realization of the boundary conditions on the curved boundary. This indicates that the treatment of boundary conditions should be further analyzed testing also other approaches as, e.g. from [11]. 8. Conclusion We have presented an efficient higher-order numerical scheme for the solution of the compressible Euler equations in time-dependent domains. It is based on several important ingredients: • • • • • •
the ALE method applied to the compressible Euler equations, the application of the discontinuous Galerkin method for the space discretization, semi-implicit time discretization, special treatment of boundary conditions, artificial viscosity applied in the vicinity of discontinuities, the use of superparametric elements near curved parts of the boundary.
The presented method behaves as unconditionally stable and appears to be robust with respect to the magnitude of the Mach number. Future work will be concentrated on the following topics: • • • •
further analysis of various treatments of boundary conditions, realization of a remeshing in case of closing the channel, coupling of the developed method with the solution of equations describing flow-induced vibrations of structures, extension of the method to the solution of compressible viscous flow.
Acknowledgments This work is a part of the research project No. MSM 0021620839 of the Ministry of Education of the Czech Republic. The research of M. Feistauer was partly supported by the Grant No. 201/08/0012 of the Czech Grant Agency. The research of J. Prokopová was supported by the Grant No. 46807 of the Grant Agency of the Charles University. The research of V. Kuˇcera was partly supported by the project LC06052 of the Ministry of Education of the Czech Republic (Jindˇrich Neˇcas Center for Mathematical Modelling). The authors acknowledge the support of these institutions. References [1] F. Bassi, S. Rebay, High-order accurate discontinuous finite element solution of the 2D Euler equations, J. Comput. Phys. 138 (1997) 251–285. [2] C.E. Baumann, J.T. Oden, A discontinuous hp finite element method for the Euler and Navier–Stokes equations, Int. J. Numer. Methods Fluids 31 (1999) 79–95. [3] B. Cockburn, G.E. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods, Lecture Notes in Computational Science and Engineering 11, Springer, Berlin, 2000. [4] T.A. Davis, I.S. Duff, A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Trans. Math. Softw. 25 (1999) 1–19. [5] V. Dolejˇsí, Anisotropic mesh adaptation for finite volume and finite element methods on triangular meshes, Comput. Vis. Sci. 1 (3) (1998) 165–178. [6] V. Dolejˇsí, M. Feistauer, A semi-implicit discontinuous Galerkin finite element method for the numerical solution of inviscid compressible flow, J. Comput. Phys. 198 (2004) 727–746. [7] V. Dolejˇsí, M. Feistauer, C. Schwab, On some aspects of the discontinuous Galerkin finite element method for conservation laws, Math. Comput. Simul. 61 (2003) 333–346.
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1623
[8] M. Feistauer, V. Dolejˇsí, V. Kuˇcera, On the discontinuous Galerkin method for the simulation of compressible flow with wide range of Mach numbers, Comput. Vis. Sci. 10 (2007) 17–27. [9] M. Feistauer, J. Felcman, I. Straˇskaraba, Mathematical and Computational Methods for Compressible Flow, Clarendon Press, Oxford, 2003. [10] M. Feistauer, V. Kuˇcera, On a robust discontinuous Galerkin technique for the solution of compressible flow, J. Comput. Phys. 224 (2007) 208–221. [11] L. Krivodonova, M. Berger, High-order accurate implementation of solid wall boundary conditions in curved geometries, J. Comput. Phys. 211 (2006) 492–512. [12] T. Nomura, T.J.R. Hughes, An arbitrary Lagrangian–Eulerian finite element method for interaction of fluid and a rigid body, Comput. Methods Appl. Mech. Eng. 95 (1992) 115–138. [13] P. Sváˇcek, M. Feistauer, J. Horáˇcek, Numerical simulation of flow induced airfoil vibrations with large amplitudes, J. Fluids Struct. 23 (2007) 391–411. [14] J.J.W. van der Vegt, H. van der Ven, Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flow, J. Comput. Phys. 182 (2002) 546–585. [15] G. Vijayasundaram, Transonic flow simulation using upstream centered scheme of Godunov type in finite elements, J. Comput. Phys. 63 (1986) 416–433.
Mathematics and Computers in Simulation 80 (2010) 1612–1623
Discontinuous Galerkin solution of compressible flow in time-dependent domains M. Feistauer ∗ , V. Kuˇcera, J. Prokopová Charles University Prague, Faculty of Mathematics and Physics, Sokolovská 83, 186 75 Praha, Czech Republic Received 30 September 2008; received in revised form 19 January 2009; accepted 29 January 2009 Available online 10 February 2009
Abstract This work is concerned with the simulation of inviscid compressible flow in time-dependent domains. We present an arbitrary Lagrangian–Eulerian (ALE) formulation of the Euler equations describing compressible flow, discretize them in space by the discontinous Galerkin method and introduce a semi-implicit linearized time stepping for the numerical solution of the complete problem. Special attention is paid to the treatment of boundary conditions and the limiting procedure avoiding the Gibbs phenomenon in the vicinity of discontinuities. The presented computational results show the applicability of the developed method. © 2009 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Compressible flow; Time-dependent domain; Fluid–structure interaction; ALE method; Discontinuous Galerkin method; Shock capturing; Flow in channels with moving walls
1. Introduction The interaction of fluid flow with vibrating bodies plays a significant role in many areas of science and technology. We mention, for example, development of airplanes (vibrations of wings) or turbines (blade vibrations), some problems from civil engineering (interaction of wind with constructions as bridges, TV towers or cooling towers of power stations), car industry (vibration of various elements of a carosery), but also medicine (hemodynamics or flow in vocal folds). In a number of these examples the moving medium is a gas, i.e. compressible flow. For low Mach number flows incompressible models are used (as e.g. in [13]), but in some cases compressibility plays an important role. The solution of fluid–structure interaction requires the coupling of the solution of equations describing the fluid flow with equations describing the structural behaviour. Due to the deformation and/or vibrations of structures, the computational domain is time dependent. There exist several techniques of the solution of incompressible flow in time-dependent domains. See, e.g. [13] and references therein. The numerical simulation of compressible flow is much more difficult, particularly in time-dependent domains. It is necessary to overcome difficulties caused by nonlinear convection dominating over diffusion, which leads to boundary layers and wakes for large Reynolds numbers and shock waves and contact discontinuities for high Mach numbers and instabilities caused by acoustic effects for low Mach numbers.
∗
Corresponding author. E-mail addresses: [email protected] (M. Feistauer), [email protected] (V. Kuˇcera), [email protected] (J. Prokopová).
0378-4754/$36.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2009.01.020
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1613
It appears that a suitable numerical method for the solution of compressible flow is the discontinuous Galerkin finite element method (DGFEM). It employs piecewise polynomial approximations without any requirement on the continuity on interfaces between neighbouring elements. The DGFEM was used for the numerical simulation of the compressible Euler equations, for example, by Bassi and Rebay in [1], where the space DG discretization was combined with explicit Runge-Kutta time discretization. In [2] Baumann and Oden describe an hp version of the space DG discretization with explicit time stepping to compressible flow. Van der Vegt and van der Ven apply space-time discontinuous Galerkin method to the solution of the Euler equations in [14], where the discrete problem is solved with the aid of a multigrid accelerated pseudo-time-integration. The papers [6,8] are concerned with a semi-implicit DGFEM technique for the solution of inviscid compressible flow, which is unconditionally stable. In [10], this method was extended so that the resulting scheme is robust with respect to the magnitude of the Mach number. For some DG techniques see [3]. The goal of our research is the numerical simulation of interaction of compressible flow with structures as, e.g. flow-induced airfoil vibrations or the flow past an elastic wall vibrating due to the influence of a moving gas. The present paper represents the first step to the solution of this problem. We are concerned with the generalization of the method from [10] to the solution of inviscid compressible flow in time-dependent domains. The main ingredients of the method is the discontinuous Galerkin space semidiscretization of the Euler equations written in the arbitrary Lagrangian–Eulerian (ALE) form, semi-implicit time discretization, suitable treatment of boundary conditions so that they are transparent for acoustic waves at the inlet and outlet and the shock capturing avoiding Gibbs phenomenon at discontinuities. Numerical experiments prove the applicability of the method. 2. Continuous problem We are concerned with inviscid compressible flow in a bounded domain t ⊂ R2 depending on time t ∈ [0, T ]. Because of simplicity we consider two-dimensional flow, but the method can be applied to 3D flow as well. We assume that the boundary of t consists of three disjoint parts I , O , Wt : ∂t = I ∪ O ∪ Wt , where I and O represent the time-independent inlet and outlet, respectively, and Wt represents moving impermeable walls. The flow is described by the Euler equations written in the conservative form ∂w ∂f s (w) + =0 ∂t ∂xs 2
in t ,
t ∈ (0, T ),
(1)
s=1
where w = (w1 , . . . , w4 )T = (ρ, ρv1 , ρv2 , E)T is the so-called state vector and f s (w) = (ρvs , ρvs v1 + δs1 p, ρvs v2 + δs2 p, (E + p) vs )T i = 1, 2, are the inviscid (Euler) fluxes of the quantity w in the directions xs , s = 1, 2. We use the following notation— ρ: density, p: pressure, E: total energy, v = (v1 , v2 ): velocity vector, γ > 1: Poisson adiabatic constant, x = (x1 , x2 ): points in t and X = (X1 , X2 ): points of the reference configuration 0 . System (1) is completed by the thermodynamical relation 1 2 p = (γ − 1) E − ρ|v| , (2) 2 boundary conditions, which are treated in Section 5, and the initial condition w(x, 0) = w0 (x),
x ∈ 0 .
(3)
We define the flux of the quantity w in the direction n = (n1 , n2 ) ∈ R2 , by F (w, n) =
2
f s (w)ns
(4)
s=1
and its Jacobi matrix DF (w, n) As (w)ns , = Dw 2
P(w, n) =
s=1
(5)
1614
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
where Df s (w) , s = 1, 2, (6) Dw are the Jacobi matrices of the mappings f s . It is possible to show that f s , s = 1, 2, are homogeneous mappings of order one, which implies that As (w) =
f s (w) = As (w)w,
s = 1, 2.
(7)
The matrix P is diagonalizable (cf. [9], Section 3.1.5), i.e. there exist matrices T = T(w, n) and = diag(λ1 , . . . , λ4 ) such that P = TT−1 .
(8)
3. ALE formulation The time dependence of the domain is taken into account with the aid of a regular one-to-one ALE mapping (cf. [12]) ¯0 → ¯ t , i.e. X ∈ ¯ 0 → x = x(X, t) = At (X) ∈ ¯ t. At :
(9)
We define the ALE velocity: ∂ At (X), ∂t z(x, t) = z˜ (A−1 t (x), t),
z˜ (X, t) =
t ∈ [0, T ], X ∈ 0 ,
(10)
t ∈ [0, T ], x ∈ t
and the ALE derivative of a function f = f (x, t) defined for x ∈ t and t ∈ [0, T ]: DA ∂f˜ f (x, t) = (X, t), Dt ∂t
(11)
where f˜ (X, t) = f (At (X), t),
X ∈ 0 , x = At (X).
(12)
By the chain rule, DA f ∂f ∂f = + z · ∇f = + div (zf ) − f div z. Dt ∂t ∂t This leads us to the ALE form of the Euler equations: DA w ∂gs (w) + + w div z = 0, Dt ∂xs
(13)
2
(14)
s=1
where gs , s = 1, 2, are modified inviscid fluxes gs (w):=f s (w) − zs w.
(15)
4. Discretization of the problem in the time-dependent domain In this section, we shall describe the discretization of the initial-boundary value problem for the Euler equations written in the ALE form (14). 4.1. Space semidiscretization The space semidiscretization will be carried out by the discontinuous Galerkin finite element method (DGFEM) based on piecewise polynomial approximations constructed over a triangular mesh without any requirement on the continuity on interfaces between neighbouring elements.
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1615
¯ ht consisting of closed Let ht be a polygonal approximation of t at time t. By Tht we denote a partition of triangles Ki ∈ Tht , i ∈ It , with disjoint interiors. (It ⊂ Z+ = {0, 1, 2, . . .} is a suitable index set.) The mesh Tht can be nonconforming with hanging nodes. By ij we denote a common edge between two neighbouring elements Ki and Kj . We set st (i) = {j ∈ It ; Kj is a neighbour of Ki }. The boundary ∂ht is formed by a finite number of faces of elements Ki adjacent to ∂ht . We denote all these boundary faces by Sj , where j ∈ Ibt ⊂ Z− = {−1, −2, . . .}. Now we set γt (i) = {j ∈ Ibt ; Sj is a face of Ki ∈ Tht } and ij = Sj for Ki ∈ Tht such that Sj ⊂ ∂Ki , j ∈ Ibt . For Ki not containing any boundary face Sj we set γt (i) = ∅. Obviously, st (i) ∩ γt (i) = ∅ for all i ∈ It . If we write St (i) = st (i) ∪ γt (i), we have ∂Ki = ij , ∂Ki ∩ ∂ht = ij . (16) j ∈ St (i)
j ∈ γt (i)
The symbol nij = ((nij )1 , (nij )2 ) will denote the unit outer normal to ∂Ki on the face ij . The approximate solution will be sought at each time instant t as an element of the finite-dimensional space S ht = S r,−1 (ht , Tht ) = {ϕh ; ϕh |K ∈ P r (K)
∀K ∈ Tht }4 ,
(17)
where r ≥ 0 is an integer and P r (K) denotes the space of all polynomials on K of degree ≤ r. Function ϕ ∈ S ht are in general discontinuous on interfaces ij . By ϕ|ij and ϕ|ji we denote the values of ϕ on ij considered from the interior and the exterior of Ki , respectively. The symbols 1 (ϕ|ij + ϕ|ji ), 2
ϕij =
[ϕ]ij = ϕ|ij − ϕ|ji
(18)
denote the avarege and jump of a function ϕ on ij = ij . In order to derive the discrete problem, we assume that w is a sufficiently regular solution of system (14), multiply (14) by a test function ϕ ∈ S ht , integrate over any element Ki , i ∈ It , apply Green’s theorem and sum over all i ∈ It . We get the relation Ki ∈ Tht
Ki
2 2 DA w(t) ∂ϕ gs (w(t)) · dx − gs (w(t))(nij )s · ϕ dS · ϕ dx = Dt ∂xs Ki ∈ Tht Ki s=1 Ki ∈ Tht j ∈ St (i) ij s=1 − div z(w · ϕ) dx. (19) Ki ∈ Tht
Ki
Now the exact solution w(t) is approximated by an element wh (t) ∈ S ht and the fluxes through the faces ij are approximated with the aid of a numerical flux H = H(u, w, n). It means that
2
ij s=1
gs (w(t)) nij
s
· ϕ dS ≈ ij
H g (wh (t)|ij , wh (t)|ji , nij ) · ϕ dS.
(20)
Here H g is an analogy to the Vijayasundaram numerical flux consistent with the fluxes gs , s = 1, 2. Taking into account (15), (5) and (6), we have Dgs (w) Df s (w) = − zs I = As − zs I, Dw Dw
(21)
and can write ˜ P(w, n) =
2 Dgs (w) s=1
Dw
ns =
2
(As ns − zs ns I) = P(w, n) − (z · n)I.
(22)
s=1
This and (8) imply that ˜ −1 , P˜ = TT
˜ = diag(λ˜ 1 , . . . , λ˜ 4 ),
λ˜ i = λi − z · n,
i = 1, . . . , 4.
(23)
1616
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
Now we define the “positive” and “negative” parts of the matrix P˜ as ˜ ± = diag(λ˜ ± , . . . , λ˜ ± ), 1 4
P˜ ± = T± T−1 ,
(24)
where λ+ = max(λ, 0), λ− = min(λ, 0). Using the above concepts, we introduce the modified Vijayasundaram numerical flux (cf. [15] or [9]) as wL + w R wL + w R H g (wL , wR , n) = P˜ + (25) , n wL + P˜ − , n wR . 2 2 If we define the forms (ψh , ϕh )h = ψh · ϕh dx,
(26)
ht
b˜ h (wh , ϕh ) = −
Ki ∈ Tht
dh (wh , ϕh ) = −
2
Ki s=1
Ki ∈ Tht
Ki
gs (wh ) ·
∂ϕh dx + ∂xs
Ki ∈ Tht j ∈ St (i) ij
H g (wh |ij , wh |ji , nij ) · ϕh dS,
divz(wh · ϕh ) dx,
(27)
(28)
we arrive at the semi-implicit discrete problem. We define an approximate solution as a function wh = wh (t) satisfying the conditions (a) (b)
wh (t) ∈ S ht , ∀t ∈ [0, T ],
DA wh (t) + b˜ h (wh (t), ϕh ) − dh (wh (t), ϕh ) = 0 , ϕh Dt
∀ϕh ∈ S ht ,
∀t ∈ (0, T ),
(29)
h
(c)
wh (0) = h w0 ,
where h w0 is the L2 (0 )-projection of w0 from the initial condition (3) on the space S h0 . 4.2. Time discretization The discrete problem (29) is equivalent to a large system of ordinary differential equations. For solving this system one can apply several numerical schemes like Runge-Kutta methods. However, they are conditionally stable and the time step is strongly limited by the CFL-stability condition. For low Mach number flows these methods fail. Therefore, we apply here a semi-implicit scheme, which is a generalization of the technique proposed in [10]. We introduce the partition 0 = t0 < t1 < . . . of the time interval [0, T ] and set τj = tj+1 − tj . The function wh (·, tj ) j will be approximated by wj , defined in tj . Let us assume that the approximate solution wh has already been computed for j = 0, . . . , k. We are interested in the computation of the approximate solution wk+1 at time instant tk+1 . If we set h ˆ h (x) = wh (Atj (A−1 w tk+1 )(x)), j
j
x ∈ htk+1 ,
(30)
then, on the basis of (11), we can approximate the ALE derivative using the first order backward difference:
ˆ kh −w wk+1 DA wh (t) h , ϕh ≈ , ϕh . Dt τk h t=tk+1
h
(31)
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1617
Further, on the basis of (7) and the definition of the modified Vijayasundaram numerical flux we define a partial linearization bh of the form b˜ h : 2 ∂ϕ (x) ˆ kh , wk+1 ˆ k (x)) − zs (x)I)wk+1 (x)) · h , ϕ bh (w ) = − ((As (w dx h h ∂xs Ki Ki ∈ Thtk+1
+
s=1
Ki ∈ Thtk+1 j ∈ Stk+1 (i) ij
˜ − ˆ k , nij )wk+1 | ) · ϕh dS. ˆ kh ij , nij )wk+1 (P˜ + (w ji h ij h |ij + P (w h (32)
The term dh linear with respect to wh will be treated implicitly. These considerations lead us to the following semi-implicit scheme: For k = 0.1, . . . find wk+1 such that h (a) (b)
∈ S htk+1 , wk+1 h
k+1 ˆ kh wh − w k+1 ˆ kh , wk+1 , ϕh + bh (w h , ϕh ) − dh (wh , ϕh ) = 0 τk
∀ϕh ∈ S htk+1 ,
(33)
h
(c)
w0h = h w0 .
This is a first order accurate scheme in time. In the solution of nonstationary flows, it is necessary to apply a scheme, which is sufficiently accurate in space as well as in time. One possibility is to apply the following two step second order time discretization based on the backward difference formula (BDF): DA wh (x, t) 2τk + τk−1 k+1 τk + τk−1 k τk ˆ h (x) + ˆ k−1 (x), x ∈ htk+1 , ≈ w (x) − w w Dt τk (τk + τk−1 ) h τk τk−1 τk−1 (τk + τk−1 ) h t=tk+1 (34) for the approximation of the ALE derivative of the solution at tk+1 and the second order approximation ˜ k+1 w = h
τk + τk−1 k τk k−1 ˆh − ˆ w w , τk−1 τk−1 h
(35)
ˆ kh in the form bh . As a result we get the of wh (tk+1 ) obtained with the aid of extrapolation, which replaces the state w ∈ S htk+1 such that following two-step second-order scheme: For each k ≥ 1 find wk+1 h 2τk + τk−1 τk + τk−1 k τk k+1 ˆ k−1 , ϕh )h ˜ k+1 ˆ h , ϕ h )h − (wk+1 (w (w h , ϕh )h + bh (w h , wh , ϕh ) = τk (τk + τk−1 ) τk τk−1 τk−1 (τk + τk−1 ) h × ∀ϕh ∈ S htk+1 ,
k = 0, 1, . . . ,
(36)
where w0h = h w0 , w1h obtained by the Runge-Kutta method. In order to guarantee the stability of the scheme, we use the CFL condition 1 τk maxKi ∈ Thtk (maxj ∈ S(i) |ij |λmax ) ≤ CFL, (37) ˆ kh |ij ,nij ) P˜ (w |Ki | where |Ki | denotes the area of Ki , |ij | is the length of the edge ij , CFL is a given constant and λmax ˆ k| P˜ (w
h ij ,nij )
is the
˜ w ˆ kh |ij , nij ), here the maximum is taken over ij . This condition is a heuristic maximal eigenvalue of the matrix P( extension of a similar condition applied in the finite volume method (see, e.g. [9]). Numerical experiments show that the presented method is practically unconditionally stable and the CFL number can be chosen large—cf. [10]. 5. Boundary conditions If ij ⊂ ∂ht , i.e. j ∈ γt (i), it is necessary to specify the boundary state w|ji appearing in the numerical flux H g in the definition of the inviscid form bh . Similarly as in [10], the choice of the boundary conditions is governed by the ˜ sign of the eigenvalues λ˜ i , i = 1, . . . , 4, of the matrix P(w, n).
1618
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
On the inlet and outlet, which are assumed fixed, we proceed in the same way as in [10], Section 4, where we prescribe the state w|ji in such a way that the locally linearized Euler equations are well posed. Using the rotational invariance, we transform the Euler equations to the coordinates x˜ 1 , parallel with the normal direction n = (n1 , n2 ) = nij to the boundary, and x˜ 2 , tangential to the boundary, neglect the derivative with respect to x˜ 2 and linearize the system around the state qij = Q(n)w|ij , where ⎛ ⎞ 1, 0, 0, 0 ⎜ 0, n , n , 0 ⎟ 1 2 ⎜ ⎟ Q(n) = ⎜ (38) ⎟ ⎝ 0, −n2 , n1 , 0 ⎠ 0,
0,
0,
1
is the rotational matrix. Then we obtain the linear system ∂q ∂q = 0, + A1 (qij ) ∂t ∂˜x1
(39)
for the transformed vector-valued function q = Q(n)w, considered in the set (−∞, 0) × (0, ∞) and equipped with the initial and boundary conditions q(˜x1 , 0) = qij , x˜ 1 < 0,
and q(0, t) = qji , t > 0.
(40)
The goal is to choose qji in such a way that this initial-boundary value problem is well posed, i.e. has a unique solution. The method of characteristics leads to the following process: Let us put q∗ji = Q(n)w∗ji , where w∗ji is a given boundary state at the inlet or outlet. We calculate the eigenvectors rs corresponding to the eigenvalues λs , s = 1, . . . , 4, of the matrix A1 (qij ), arrange them as columns in the matrix T and calculate T−1 . Now we set α = T−1 qij ,
β = T−1 q∗ji .
and define the state qji by the relations 4 αs , λs ≥ 0, qji := γs rs , γs = βs , λs < 0. s=1
(41)
(42)
Finally, the sought boundary state w|ji is defined as w|ji = wji = Q−1 (n)qji .
(43)
On the impermeable moving wall we prescribe the normal component of the velocity v · n = z · n,
(44)
˜ n) vanish, where n is the unit outer normal to Wt and z is the wall velocity. This means that two eigenvalues of P(w, one is positive and one is negative. Then, in analogy to [9], Section 3.3.6, we should prescribe one quantity, namely v · n, and extrapolate three quantities—tangential velocity, density and pressure. However, here we define the numerical flux on Wt as the physical flux through the boundary with the assumption (44) taken into account. Thus, on Wt we write 2
gs (w)ns = (v · n − z · n)w + p (0, n1 , n2 , v · n)T = p (0, n1 , n2 , z · n)T =: H g .
(45)
s=1
As was shown in [1] or [9], Section 4.6.8, the piecewise linear approximation of a curved boundary can lead to a spurious nonphysical entropy production near vertices. In order to avoid this drawback, isoparametric curved elements should be used. In our code we follow [9], Section 4.6.8, where the reference triangle is transformed by a bilinear mapping onto the approximation of a curved triangle adjacent to the boundary ∂.
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1619
Fig. 1. Triangulations at time instants t = 0.0, 6.1392 and 10.0592.
6. Limiting procedure at discontinuities For high speed flows with shock waves and contact discontinuities it is necessary to avoid the Gibbs phenomenon manifested by spurious overshoots and undershoots in computed quantities near discontinuities. In order to avoid the Gibbs phenomenon, we apply the limiting procedure from [10] based on the discontinuity indicator 2 [ρˆ hk ] dS gk (i) = , Ki ∈ Thtk+1 , (46) 3/4 ∂Ki hKi |Ki | ˆ kh . Then we define the discrete introduced in [7]. The density ρˆ hk represents the first component of the state vector w discontinuity indicator Gk (i) = 0,
if gk (i) < 1,
Gk (i) = 1,
if gk (i) ≥ 1,
and the artificial viscosity forms k ˆ kh , wk+1 βh (w , ϕ ) = ν h G (i) 1 Ki h h i ∈ Itk+1
and ˆ kh , wk+1 Jh (w h , ϕh ) = ν2
i ∈ I j ∈ stk+1 (i)
Ki
Ki ∈ Thtk+1 ,
∇wk+1 · ∇ϕh dx h
1 k (G (i) + Gk (j)) 2
(47)
(48)
ij
[wk+1 h ] · [ϕh ] dS,
(49)
with ν1 , ν2 = O(1). Then the resulting scheme obtained by limiting of (33)(b) has the following form: we seek wk+1 ∈ S htk+1 such that h
ˆ kh −w wk+1 h ˆ kh , wk+1 ˆ kh , wk+1 ˆ kh , wk+1 , ϕh + bh (w h , ϕh ) + βh (w h , ϕh ) + Jh (w h , ϕh ) = 0, ∀ ϕh ∈ S htk+1 , τk h
k = 0, 1, . . . , (Similarly we obtain a stabilized version of scheme (36).)
1620
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
Fig. 2. Velocity (left) and pressure (right) isolines at time instants t = 0.3992, 2.2192, 3.1992, 4.1792, 6.1392 and 7.1192.
In practical computations, the integrals appearing in the definition of the approximated solution are evaluated with the aid of quadrature formulae. Linear algebraic systems are solved either with the aid of the direct solver UMFPACK [4] or by the GMRES method with block diagonal preconditioning. 7. Numerical experiments We consider inviscid compressible flow in the rectangular channel with the initial shape 0 = [−2, 2] × [0, 1], where the lower wall of the channel is moving. The ALE mapping is equal to identity in the sets [−2, −1] × [0, 1] and
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1621
Fig. 3. Velocity (left) and pressure (right) isolines at time instants t = 8.0992, 10.0592, 12.0192, 13.9792, 14.9592 and 15.9392.
[1, 2] × [0, 1]. In the rest of the channel we construct the ALE mapping so that the lower wall is represented at time t by the graph of the function 0.45 sin(0.4t) (cos(πX1 ) + 1),
X1 ∈ (−1, 1).
(50)
This movement is interpolated to the rest of the domain resulting in the ALE mapping At . The computation was carried out with the dimensionless form of the Euler equations, using the dimensionless conservation variables, on a triangulation with 1160 elements. Fig. 1 shows the triangulation at time t = 0.0, 6.1392 and 10.0592. The triangular mesh was constructed by the technique from [5]. The inlet Mach number Minlet = 0.067,
1622
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
the dimensionless inlet density ρinlet = 1.0. At the outlet, the dimensionless pressure poutlet = 159.12. In the artificial viscosity forms (48) and (49) the values ν1 = ν2 = 0.2 were used. Figs. 2 and 3 show velocity and pressure isolines at different time instants t = 0.3992, 2.2192, 3.1992, 4.1792, 6.1392, 7.1192, 8.0992, 10.0592, 12.0192, 13.9792, 14.9592 and 15.9392 during one period. In the solution we can observe a vortex formation, when the lower wall starts to descend. This vortex is convected through the domain. Moreover, we see that a contact discontinuity is developed, when the channel becomes narrow. The contact discontinuity is characterized by the discontinuity in the velocity, whereas the pressure remains continuous. In Fig. 2, one can see some “artefacts” in pressure isolines near the moving wall at time t = 6.1392, which might be caused by the realization of the boundary conditions on the curved boundary. This indicates that the treatment of boundary conditions should be further analyzed testing also other approaches as, e.g. from [11]. 8. Conclusion We have presented an efficient higher-order numerical scheme for the solution of the compressible Euler equations in time-dependent domains. It is based on several important ingredients: • • • • • •
the ALE method applied to the compressible Euler equations, the application of the discontinuous Galerkin method for the space discretization, semi-implicit time discretization, special treatment of boundary conditions, artificial viscosity applied in the vicinity of discontinuities, the use of superparametric elements near curved parts of the boundary.
The presented method behaves as unconditionally stable and appears to be robust with respect to the magnitude of the Mach number. Future work will be concentrated on the following topics: • • • •
further analysis of various treatments of boundary conditions, realization of a remeshing in case of closing the channel, coupling of the developed method with the solution of equations describing flow-induced vibrations of structures, extension of the method to the solution of compressible viscous flow.
Acknowledgments This work is a part of the research project No. MSM 0021620839 of the Ministry of Education of the Czech Republic. The research of M. Feistauer was partly supported by the Grant No. 201/08/0012 of the Czech Grant Agency. The research of J. Prokopová was supported by the Grant No. 46807 of the Grant Agency of the Charles University. The research of V. Kuˇcera was partly supported by the project LC06052 of the Ministry of Education of the Czech Republic (Jindˇrich Neˇcas Center for Mathematical Modelling). The authors acknowledge the support of these institutions. References [1] F. Bassi, S. Rebay, High-order accurate discontinuous finite element solution of the 2D Euler equations, J. Comput. Phys. 138 (1997) 251–285. [2] C.E. Baumann, J.T. Oden, A discontinuous hp finite element method for the Euler and Navier–Stokes equations, Int. J. Numer. Methods Fluids 31 (1999) 79–95. [3] B. Cockburn, G.E. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods, Lecture Notes in Computational Science and Engineering 11, Springer, Berlin, 2000. [4] T.A. Davis, I.S. Duff, A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Trans. Math. Softw. 25 (1999) 1–19. [5] V. Dolejˇsí, Anisotropic mesh adaptation for finite volume and finite element methods on triangular meshes, Comput. Vis. Sci. 1 (3) (1998) 165–178. [6] V. Dolejˇsí, M. Feistauer, A semi-implicit discontinuous Galerkin finite element method for the numerical solution of inviscid compressible flow, J. Comput. Phys. 198 (2004) 727–746. [7] V. Dolejˇsí, M. Feistauer, C. Schwab, On some aspects of the discontinuous Galerkin finite element method for conservation laws, Math. Comput. Simul. 61 (2003) 333–346.
M. Feistauer et al. / Mathematics and Computers in Simulation 80 (2010) 1612–1623
1623
[8] M. Feistauer, V. Dolejˇsí, V. Kuˇcera, On the discontinuous Galerkin method for the simulation of compressible flow with wide range of Mach numbers, Comput. Vis. Sci. 10 (2007) 17–27. [9] M. Feistauer, J. Felcman, I. Straˇskaraba, Mathematical and Computational Methods for Compressible Flow, Clarendon Press, Oxford, 2003. [10] M. Feistauer, V. Kuˇcera, On a robust discontinuous Galerkin technique for the solution of compressible flow, J. Comput. Phys. 224 (2007) 208–221. [11] L. Krivodonova, M. Berger, High-order accurate implementation of solid wall boundary conditions in curved geometries, J. Comput. Phys. 211 (2006) 492–512. [12] T. Nomura, T.J.R. Hughes, An arbitrary Lagrangian–Eulerian finite element method for interaction of fluid and a rigid body, Comput. Methods Appl. Mech. Eng. 95 (1992) 115–138. [13] P. Sváˇcek, M. Feistauer, J. Horáˇcek, Numerical simulation of flow induced airfoil vibrations with large amplitudes, J. Fluids Struct. 23 (2007) 391–411. [14] J.J.W. van der Vegt, H. van der Ven, Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flow, J. Comput. Phys. 182 (2002) 546–585. [15] G. Vijayasundaram, Transonic flow simulation using upstream centered scheme of Godunov type in finite elements, J. Comput. Phys. 63 (1986) 416–433.