- Email: [email protected]
Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow
Accepted Manuscript
Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow ˇ Monika Balazsov a, Miloslav Feistauer, Petr Svaˇ ´ ´ Jan Cesenek, ´ cek, Jarom´ır Horaˇ ´ cek PII: DOI: Reference:
S0045-7930(18)30391-8 10.1016/j.compfluid.2018.07.002 CAF 3943
To appear in:
Computers and Fluids
Received date: Revised date: Accepted date:
22 September 2017 28 June 2018 2 July 2018
ˇ Please cite this article as: Monika Balazsov a, Miloslav Feistauer, Petr Svaˇ ´ ´ Jan Cesenek, ´ cek, Jarom´ır Horaˇ ´ cek, Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow, Computers and Fluids (2018), doi: 10.1016/j.compfluid.2018.07.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • viscous incompressible and viscous compressible flow for low Mach numbers
• flow in a simplified channel of the human vocal tract • flow around an airfoil
AC
CE
PT
ED
M
AN US
• flow-induced airfoil vibrations
CR IP T
• robustness of the discontinuous Galerkin method with respect to low Mach numbers
1
ACCEPTED MANUSCRIPT
Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow✩
CR IP T
2 ˇ Monika Bal´azsov´a1 , Jan Cesenek , Miloslav Feistauer1 , Petr Sv´acˇ ek3 , Jarom´ır Hor´acˇ ek4
Abstract
The novelty of the paper is the application of the discontinuous Galerkin method (DGM) to the solution of viscous compressible flow with low Mach number (M ≈ 10−2 − 10−3 ). The compressible Navier-Stokes
AN US
equations are written in the conservative form using conservative variables. They are discretized by the DGM in space. In time we use either the backward difference formula (BDF) or the discontinuous Galerkin method. The resulting technique is called the space-time discontinuous Galerkin method (STDGM). The goal of the presented analysis is to prove the ability of the developed DGM to solve viscous compressible low Mach number flows described by the compressible Navier-Stokes equations without any
M
modification and to obtain results comparable with simulations of incompressible viscous flow described by the incompressible Navier-Stokes equations, solved by the finite element method (FEM). We present
ED
numerical results for a flow in a simplified channel of the human vocal tract, for a flow around the airfoil NACA0012 and moreover for flow-induced airfoil vibrations. The results prove the robustness of the discontinuous Galerkin method with respect to the low Mach number.
PT
Keywords: incompressible Navier-Stokes equations, compressible Navier-Stokes equations, low Mach number flow, finite element method, discontinuous Galerkin method, fluid-structure interaction
CE
2000 MSC: 65N50, 65N15, 65D05
AC
1. Introduction
In the numerical solution of compressible flow, it is necessary to overcome several issues: resolve accu-
rately shock waves, contact discontinuities and boundary layers. All these phenomena are connected with ✩ The work of M. Bal´ azsov´a was supported by the Grant No. SVV-2017-260455 of the Charles University. The research of M. Feistauer was supported by the Grant No. 17-01747S of the Czech Science Foundation and finally the research of P. Sv´acˇ ek and J. Hor´acˇ ek was supported by the Grant No. 16-01246S of the Czech Science Foundation. ˇ Email addresses: [email protected] (Monika Bal´azsov´a), [email protected] (Jan Cesenek), [email protected] (Miloslav Feistauer), [email protected] (Petr Sv´acˇ ek), [email protected] (Jarom´ır Hor´acˇ ek) 1 Charles University Prague, Faculty of Mathematics and Physics, Sokolovsk´ a 83, 186 75 Praha, Czech Republic 2 Aeronautical Research and Test Institute, Beranov´ ych 130, 199 05 Prague 9, Czech Republic 3 Department of Technical Mathematics, Faculty of Mechanical Engineering, Czech Technical University, Karlovo n´ amˇest´ı 13, 121 35 Prague 2, Czech Republic 4 Institute of Termomechanics, Academy of Sciences of the Czech Republic, Dolejˇskova 5, 182 00 Prague 8, Czech Republic
Preprint submitted to Computers & Fluids
June 28, 2018
ACCEPTED MANUSCRIPT
the simulation of high speed flow with high Mach numbers. However, it appears that the solution of low Mach number flow is also rather difficult. This is caused by the stiff behaviour of numerical schemes and acoustic phenomena appearing in low Mach number flows at incompressible limit. This led to the development of special finite volume techniques allowing the simulation of compressible flow at incompressible
CR IP T
limit, which are based on modifications of Euler equations describing inviscid compressible flow. From a wide literature we mention, for example, works [17], [25], [20], Chapter 5, or [32], Chapter 14). However, these techniques could not be applied to the solution of high speed flow. Therefore, further attempts were concentrated on the extension of these methods to the solution of flows at all speeds. A success in this direction was achieved by several authors. Let us mention, for example, the works by Wesseling et al. (e. g. [30],
AN US
[31], [33]), Parker and Munz ([23]), Meister ([18], [19]), Nerinckx et al. ([21]) and Darwish et al. ([7]). Main ingredients of these techniques are finite volume schemes applied on staggered grids, combined with multigrid, the use of the pressure-correction, multiple pressure variables and flux preconditioning. The paper [27] uses the immersed boundary method. There is a question, if it is possible to develop a numerical method allowing the solution of compressible flow with a wide range of Mach numbers without modifi-
M
cation of the governing equations. This was successful in papers [13] and [12] dealing with the solution of the inviscid compressible Euler equations with a wide range of Mach numbers. This technique is based
ED
on the discontinuous Galerkin method (DGM), which combines advantages of the finite volume as well as finite element methods. It employs piecewise polynomial approximations without any requirement on the continuity on interfaces between neighbouring elements.
PT
The DGM was first applied to the numerical simulation of the compressible Euler equations, by Bassi and Rebay in [1], where the space discontinuous Galerkin discretization was combined with explicit Runge-
CE
Kutta time discretization. In [2] Baumann and Oden apply an hp version for the solution of the space DG discretization with explicit time stepping to compressible flow. Van der Vegt and van der Ven use the space-
AC
time discontinuous Galerkin method for the solution of the Euler equations in [29], where the discrete problem is solved with the aid of a multigrid accelerated pseudo-time-integration. The DGM applied to compressible flow was combined with adaptive techniques (see, e.g. [16]). An extension of results from [2] is the subject of [9]. In [11] the DGM was combined with the PN P M method. It is possible to say that
during the last decades, a number of papers devoted to the solution of compressible flow by the DGM was published. An older survey can be found, for example, in [6], whereas recent references are mentioned in [10]. The subject of this paper is the numerical simulation of viscous low Mach number compressible flow described by the compressible Navier-Stokes equations and compared with results obtained with the use of 2
ACCEPTED MANUSCRIPT
incompressible Navier-Stokes equations. The novelty of the presented research lies in the demonstration of the applicability of the developed DG method for the numerical solution of compressible low Mach number flow, practically at incompressible limit. The robustness of the DGM with respect to the Mach number is proved by the comparison of results with the solution of incompressible Navier-Stokes equations and in
CR IP T
several cases with Sandia experimental data.
In Section 2 we formulate the continuous problem of incompressible as well as compressible viscous flow. Section 3 is devoted to the discretization of the problem. For the discretization of the incompressible flow we use standard conforming stabilized finite element method (FEM). The discretization of compressible flow is carried out by the discontinuous Galerkin method in space. In both cases, for the time dis-
AN US
cretization we use the backward difference formula (BDF). In the case of the compressible Navier-Stokes equations, also discontinuous Galerkin in time is used. Then we get the so-called space-time discontinuous Galerkin method (STDGM). Section 4 presents the description of the simulation of flow induced airfoil vibrations. Here the above problems are reformulated in time-dependent domains with the aid of the arbitrary Lagrangian-Eulerian (ALE) method. In Section 5 we present several numerical experiments showing
M
the comparison of incompressible flow and compressible viscous flow with low Mach numbers. In the first case we compare incompressible flow and compressible viscous flow in a fixed 2D channel, which repre-
ED
sents a simplified model of the glottal region of the human vocal tract, i.e. airways in the human subglottal and supraglottal spaces with a sudden constriction at the glottis. Our choice of the domain was motivated by number of works on fluid-structure interaction problems for vocal folds, e.g. [14], [4] and [24]. In the
PT
second case we investigate the incompressible flow and the compressible viscous flow around the fixed airfoil NACA0012 for different angles of attack. Moreover, we are concerned with the simulation of airfoil
CE
vibrations induced by incompressible flow and low Mach number compressible flow.
AC
2. Governing equations The aim of this paper is to show the ability of the DGM to approximate compressible flow for low
Mach numbers, particularly in the case when the flow is almost incompressible. To this end, numerical results for the compressible flow for low Mach numbers are compared with the numerical results for the incompressible flow. In what follows, the incompressible and compressible flow models are formulated.
3
ACCEPTED MANUSCRIPT
2.1. Incompressible flow The incompressible viscous flow in a computational polygonal domain Ω ⊂ R2 and a time interval [0, T ], T > 0 is governed by the Navier-Stokes equations in the dimensionless form
∇ · v = 0,
(2.1)
CR IP T
∂v 1 + (v · ∇)v + ∇p − △v = 0, ∂t Re
where v = (v1 , v2 ) = v(x, t) denotes the velocity vector, p = p(x, t) denotes the pressure, Re is the Reynolds number, defined as Lre f vre f ρre f , µ
(2.2)
AN US
Re =
where Lre f is the reference length, vre f is the reference velocity, ρre f is the reference density (for incompressible flow it is the constant fluid density) and µ is the dynamic viscosity of the fluid. System (2.1) is equipped with the initial condition
M
v(x, 0) = v0 (x)
x ∈ Ω,
(2.3)
ED
and with the boundary conditions
a) v = vin = (vin1 , vin2 )T
on ΓI ,
b) v = 0
PT
1 ∂v 1 =0 c) − (p − pout )n − (v · n)− v + 2 Re ∂n
on ΓW ,
(2.4)
on ΓO ,
CE
where pout is a prescribed pressure distribution on the outlet and n denotes the unit outer normal vector to the boundary ∂Ω = ΓI ∪ ΓO ∪ ΓW , where ΓI denotes the inlet, ΓO the outlet and ΓW the fixed impermeable
AC
walls. The boundary condition (2.4c) is a modification of the so-called ’do-nothing’ boundary condition, cf. [3]. We set λ− = min{0, λ} and λ+ = max{0, λ} for λ ∈ R. 2.2. Compressible flow The compressible viscous flow in the domain Ω ⊂ R2 and the time interval [0, T ] is described by the system of the Navier-Stokes equations in the dimensionless form ∂w X ∂ f s (w) + ∂t ∂x s s=1 2
=
2 X ∂R s (w, ∇w) s=1
4
∂x s
,
(2.5)
ACCEPTED MANUSCRIPT
where w is the so called state vector defined as w = (w1 , w2 , w3 , w4 )T = (ρ, ρv1 , ρv2 , E)T ∈ R4 and w0 (x). Functions
CR IP T
f s (w) = (ρv s , ρv1 v s + δ1s p, ρv2 v s + δ2s p, (E + p)v s )T , !T γ ∂θ V V V V R s (w, ∇w) = 0, τ s1 , τ s2 , τ s1 v1 + τ s2 v2 + Re Pr ∂x s
represent the inviscid and viscous terms, respectively. We use the following notation: ρ - density, δi, j Kronecker symbol, E - total energy, θ - absolute temperature, τV is the viscous part of the stress tensor with components
! 1 ∂vi ∂v j 2 + − δi j ∇ · v , = Re ∂x j ∂xi 3
i, j = 1, 2
AN US
τVij
(2.6)
and τ denotes the stress tensor with components
τi j = −pδi j + τVij ,
i, j = 1, 2.
(2.7)
by (2.2), Pr =
cpµ k
M
Further, we use the following notation: γ > 1 - Poisson adiabatic constant, Re - Reynolds number defined - Prandtl number, where c p is the specific heat at constant pressure and k is the heat
ED
conductivity.
The above system is completed by the thermodynamical relations
PT
! 1 2 p = (γ − 1) E − ρ|v| , 2
θ=
E 1 2 − |v| . ρ 2
(2.8)
CE
The resulting system is equipped with the initial condition
x ∈ Ω0 ,
(2.9)
AC
w(x, 0) = w0 (x),
and boundary conditions
a) ρ|ΓI = ρin , b) v|ΓI = vin , 2 X ∂θ = 0 on ΓI , c) τVij ni v j + k ∂n i, j=1 d) v = 0, f)
2 X
e)
∂θ =0 ∂n
τVij ni = 0,
on ΓW ,
j = 1, 2,
i=1
5
(2.10)
g)
∂θ = 0 on ΓO . ∂n
ACCEPTED MANUSCRIPT
Since we consider low Mach number flow, we also consider the following condition
p = pout
on ΓO ,
(2.11)
CR IP T
see also [10], Sections 9.1.2 and 8.3. 2.3. Flow in a time-dependent domain Ωt
In this case the impermeable wall ΓWt of the domain occupied by the fluid may depend on time t ∈ [0, T ] and the boundary conditions (2.4), b) and (2.10), d) are replaced by the condition
AN US
v|ΓWt = zD = (zD1 , zD2 ) − velocity of the moving wall ΓWt .
(2.12)
In order to treat the time dependence of the domain, we use the so-called arbitrary Lagrangian-Eulerian (ALE) technique, proposed in [22]. We define a reference domain Ω0 and introduce a regular one-to-one
M
ALE mapping of Ω0 onto Ωt (cf. [22], [28] and [29])
ED
At : Ω0 −→ Ωt , i.e. X ∈ Ω0 7−→ x = x(X, t) = At (X) ∈ Ωt . Here we use the notation X for points in Ω0 and x for points in Ωt .
CE
PT
Further, we define the domain velocity:
˜z(X, t) =
∂ At (X), ∂t
z(x, t) = ˜z(At−1 (x), t),
t ∈ [0, T ], X ∈ Ω0 , t ∈ [0, T ], x ∈ Ωt
(2.13) (2.14)
AC
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
(2.15)
where f˜(X, t) = f (At (X), t), X ∈ Ω0 , x = At (X). As a direct consequence of the chain rule we get the relation ∂f DA f = + div (z f ) − f div z. Dt ∂t
6
(2.16)
ACCEPTED MANUSCRIPT
This leads to the ALE formulation of the incompressible Navier-Stokes equations in the form DA v 1 + ((v − z) · ∇)v + ∇p − △v = 0. Dt Re
(2.17)
CR IP T
The compressible Navier-Stokes equations have the ALE form X ∂R s (w, ∇w) DA w X ∂ g s (w) + + w divz = , Dt ∂x s ∂x s s=1 s=1 2
2
where s = 1, 2,
(2.19)
AN US
g s (w) := f s (w) − z s w,
(2.18)
are the ALE modified inviscid fluxes.
3. Discretization
The space semidiscretization of the above initial boundary-value problems is carried out on a partition
M
Th of the domain Ω, formed by a finite number of closed triangles with disjoint interiors. We assume that the intersection of two elements K, K ′ ∈ Th is either empty or a common edge or a common vertex of these
ED
elements and the end points of ΓI , ΓO , ΓW are vertices of some elements K ∈ Th . In order to discretize the problem in time, we consider an equidistant partition 0 = t0 < t1 < · · · < tN =
PT
T, tk = k∆t of the time interval [0, T ] with a time step ∆t > 0. 3.1. Incompressible flow
AC
CE
The time derivative in (2.1) is then approximated by the second order backward difference formula 3vn+1 − 4vn + vn−1 ∂v , t=t ≈ ∂t n+1 2∆t
(3.20)
where vk denotes the approximation of v(tk ) for k = 0, . . . , N, and we set v−1 := v0 . Now, let us focus on the space discretization at a given time level tn+1 using the finite element method. The standard conforming finite element spaces defined with the aid of the above mentioned admissible triangulation Th of Ω are used. Particularly, the finite element spaces for the flow velocity Yh and for pressure Qh are defined as i2 h Yh = ϕ = (ϕ1 , ϕ2 ) ∈ C(Ω) : ϕi ∈ P2 (K), ∀K ∈ Th , i = 1, 2 , K
7
ACCEPTED MANUSCRIPT
and
Qh = ϕ ∈ C(Ω) : ϕ ∈ P1 (K), ∀K ∈ Th , K
where by the symbol Pr (K) the space of all polynomials on the triangle K of degree less or equal r is
CR IP T
denoted. Further, the space Wh of test functions is defined as Wh = {ϕ ∈ Yh : ϕ = 0 on ΓI ∪ ΓW } .
The finite element approximation is introduced in the standard way using the multiplication of equations (2.1) by the test functions ϕ ∈ Wh and q ∈ Qh , integration over Ω and using the Green theorem together
AN US
with the boundary condition (2.4c). This leads to the definition of the forms Z
1 1 3 v · ϕ + ((v · ∇)v) · ϕ − ((v · ∇)ϕ) · v − p(∇ · ϕ) dx 2∆t 2 2 Ω Z Z 1 1 (v · n)+ v · ϕ dS + ∇v · ∇ϕ + (∇ · v)q dx, + 2 Re Ω ΓO Z 1 (4vn − vn−1 ) · ϕ dx, f (V) = 2∆t Ω =
(3.21)
(3.22)
M
a(U, V)
for any U = (v, p) ∈ Yh × Qh and V = (ϕ, q) ∈ Wh × Qh . In order to overcome the possible instability of
ED
the convection dominated flow in the case of Galerkin approximations, the residual based stabilization is employed, cf. [15]. It uses the SUPG/PSPG stabilization terms defined by X
PT
L(U, V)
=
CE
F (U; V) =
δK
K∈Th
X
K∈Th
δK
Z
Z
K
K
! 1 3v − △v + (v · ∇)v + ∇p · (v · ∇) ϕ + ∇q dx, 2∆t Re
4vn − vn−1 · (v · ∇)ϕ + ∇q dx, 2∆t
(3.23) (3.24)
AC
in combination with the stabilization term for div-div stabilization given as
P(U, V) =
X
τK
K∈Th
Z
K
(∇ · v)(∇ · ϕ) dx.
(3.25)
The choice of the parameters δK = h2K and τK = 1 is carried out according to [15] or [28], where hK denotes the local element length. Using the above defined non-linear forms, the non-linear stabilized discrete problem at a time instant t = tn+1 reads: Definition An approximate solution of problem (2.1), (2.4) is defined as U = (v, p) ∈ Yh ×Qh , p := pn+1 , 8
ACCEPTED MANUSCRIPT
v := vn+1 , such that v satisfies approximately the conditions (2.4 a-b) and
a(U, V) + L(U, V) + P(U, V) = f (V) + F (U; V),
CR IP T
holds for all V = (ϕ, q) ∈ Wh × Qh .
(3.26)
The solution of the arising system of non-linear equations is realized by the Oseen linearization. The linearized problem is solved with the aid of a direct solver. See, e.g., [28]. 3.2. Compressible flow
For the space discretization we use the discontinuous Galerkin method (DGM) applied on the triangu-
AN US
lation Th , introduced in the previous section. Furthermore, we consider a partition 0 = t0 < t1 . . . < T as above with a time step ∆t. Let us denote by Fh the set of all faces of elements from Th . It consists of the set of all interior and boundary faces: Fh = FhI ∪ FhB . Each face Γ ∈ Fh is associated with a unit normal nΓ , which is an outer normal to ∂Ω for Γ ∈ FhB . The approximate solution is sought in the space of piecewise
M
polynomial functions
Shr = (S hr )4 ,
S hr = {v; v|K ∈ Pr (K) ∀K ∈ Th }.
(3.27)
ED
where
The symbols [u]Γ and huiΓ will denote the jump and the mean value of u ∈ Shr on the face Γ ∈ FhI and [u]Γ = huiΓ = u|Γ for Γ ∈ FhB .
PT
In order to derive the space semidiscrete problem, we assume that there exists an exact solution w ∈ C 1 ([0, T ]; Shr ) of the Navier-Stokes equation (2.5). We multiply (2.5) by a test function ϕ h ∈ Shr , integrate
AC
CE
over an element K ∈ Th , apply Green’s theorem and sum over all elements. Then we can formally write X Z ∂w · ϕ dx + Inv + Vis = 0, K ∂t K∈T
(3.28)
h
where
Inv =
XZ
K∈Th
Vis = −
∂K s=1
XZ
K∈Th
2 X
f s (w)n s · ϕ dS −
2 X
∂K s=1
2 XZ X
K∈Th
R s (w, ∇w)n s · ϕ dS +
K s=1
f s (w) ·
2 XZ X
K∈Th
K s=1
ϕ ∂ϕ dx, ∂x s
R s (w, ∇w) ·
(3.29) ϕ ∂ϕ dx, ∂x s
represent the inviscid and viscous terms and n = (n1 , n2 ) is the outer unit normal to ∂K. To discretize the inviscid term, we use the so-called numerical flux H, which has following properties: 9
ACCEPTED MANUSCRIPT
H(w1 , w2 , n) is locally Lipschitz-continuous with respect to w1 , w2 , P H(w, w, n) is consistent: H(w, w, n) = ds=1 f (w)n s ,
(H1) (H2)
H(w1 , w2 , n) is conservative: H(w1 , w2 , n) = −H(w2 , w1 , −n).
(H3)
Z X 2
f s (w)n s · ϕ dS ≈
Γ s=1
Z
Γ
CR IP T
Integrals over faces Γ ∈ Fh are approximated with the aid of the numerical flux H: (R) H(w(L) Γ , wΓ , nΓ ) ϕ dS .
(3.30)
Then the inviscid terms Inv are discretized by convection form bh (w, ϕ ). Hence
AN US
Inv ≈ bh (w, ϕ ), where XZ
Γ
Γ∈FhI
−
ϕ] dS + nΓ ) [ϕ
XZ
Γ∈FhB
2 XZ X
f s (w) ·
Γ
(L) H(w(L) Γ , wΓ , nΓ ) ϕ dS
(3.31)
ϕ ∂ϕ dx. ∂x s
M
bh (w, ϕ ) =
(R) H(w(L) Γ , wΓ ,
K s=1
K∈Th
ED
For details we can refer to the paper [13].
Now we can turn our attention to the approximation of the viscous term Vis. It is based on the relation
PT
R s (w, ∇w) =
2 X
K s,k (w)
k=1
∂w , ∂xk
s = 1, 2,
(3.32)
CE
where K s,k are 4 × 4 matrices dependent on w. For more details see [10].
AC
ah (w, ϕ ) =
2 XZ X
K∈Th
−
K s,k=1
! + 2 *X 2 XZ X ϕ ∂w ∂w ∂ϕ ϕ] dS (3.33) K s,k (w) K s,k (w) · dx − n s · [ϕ ∂xk ∂x s ∂xk Γ s=1 k=1 I
2 XZ X
Γ∈FhB
−θ
Γ s,k=1
Γ∈Fh
K s,k (w)
2 XZ X
Γ∈FhB
Γ s,k=1
+ 2 *X 2 XZ X ∂w ∂ϕ KTs,k (w) n s · ϕ dS − θ n s · [w] dS ∂xk ∂xk Γ s=1 k=1 I Γ∈Fh
KTs,k (w)
∂ϕ n s · (w − wB ) dS , ∂xk
where θ = 1, 0 or −1 (SIPG, IIPG, resp. NIPG variants) and wB is a boundary state obtained by Dirichlet data and extrapolation. In the scheme we include interior and boundary penalty terms, vanishing for a sufficiently smooth exact
10
ACCEPTED MANUSCRIPT
solution satisfying the boundary conditions, defined as XZ
Γ∈FhI
Γ
ϕ]dS + σ[w] · [ϕ
XZ
Γ∈FhB
Γ
σ(w − wB ) · ϕ dS ,
where the penalty weight σ is chosen as
σ|Γ =
CW . Γ ∈ Fh , diam(Γ) Re
(3.34)
CR IP T
Jh (w, ϕ ) =
with a suitable constant CW > 0, which guarantees the stability of the method. Its choice depends on the
AN US
variant of the used DGM (NIPG, IIPG or SIPG). Now we approximate the viscous term in the following form:
Vis ≈ ah (w, ϕ) + Jh (w, ϕ)
M
Finally we set
ϕ ∈ Shr . ∀ϕ
w, ϕ ∈ Shr ,
(3.35)
ED
ch (w, ϕ ) = bh (w, ϕ ) + ah (w, ϕ ) + Jh (w, ϕ ),
and arrive at the definiton of an approximate solution. Definition We say that the functions wh is the space semidiscrete solution of our problem, if the fol-
AC
CE
PT
lowing conditions are satisfied:
wh ∈ C 1 ([0, T ]; Shr ),
(3.36)
wh (0) = Πh w0 ,
(3.37)
d (wh , ϕh )Ω + ch (wh , ϕh ) = 0 dt
ϕh ∈ Shr , ∀t ∈ (0, T ), ∀ϕ
(3.38)
where Πh w0 is an Shr -approximation of w0 from the initial condition. In the time discretization we use the second-order backward difference formula (BDF) as for the in-
compressible flow. Definition We say that the finite sequence of functions wkh = wh (tk ),
11
k = 0, . . . , N is the approximate
ACCEPTED MANUSCRIPT
solution of (5)–(10) computed by the BDF-DGM, if the following conditions are satisfied:
0 We set w−1 h = wh .
(3.39) ϕh ∈ Shr , ∀ϕ
n = 1, 2, . . . .
(3.40)
CR IP T
w0h = Πh w0 , wn+1 ∈ Shr , n = 0, 1, . . . , h n+1 n n−1 3wh − 4wh + wh , ϕ h + ch (wn+1 h , ϕh) = 0 2∆t Ω
Another technique for the time discretization is the time discontinuous Galerkin method. Then we obtain the space-time discontinuous Galerkin method (STDGM). We again consider a partition 0 = t0 < t1 < . . . < tN = T of the time interval [0, T ] and denote Im = (tm−1 , tm ), I m = [tm−1 , tm ]. We define the
r,q S h,∆t
AN US
r,q r,q 4 space Sh,∆t = (S h,∆t ) , where
q X = φ ; φ| = ζ φ , where φ ∈ S , ζ ∈ P (I ) i i i hr i q m I m i=0
with integers p, q ≥ 1, Pq (Im ) denoting the space of all polynomials in t on Im of degree ≤ q and the space
M
r,q S hr defined in (3.27). For ϕ ∈ Sh,∆t we introduce the following notation:
ϕ±m = ϕ(tm± ) = lim ϕ(t),
ED
t→tm ±
{ϕ}m = ϕ+m − ϕ−m .
(3.41)
Definition The space-time discontinuous Galerkin approximate solution of (5)–(10) is defined as a
PT
function wh∆t satisfying
wh∆t (0+), ϕh
Ω
= w0 , ϕh
Ω
∀ϕh ∈ Shr
(3.42)
AC
CE
and the following conditions:
r,q 1) wh∆t ∈ Sh,∆t , ! ! Z ∂wh∆t (t), ϕh∆t + ah (wh∆t , ϕh∆t , t) dt 2) ∂t I Ω Zm + bh (wh∆t , ϕh∆t , t) + Jh (wh∆t , ϕh∆t , t) dt + ({wh∆t }m−1 , ϕh∆t (tm−1 +))Ω Im Z r,q = ℓh (whD , ϕh∆t , t) dt ∀ϕh∆t ∈ Sh,∆t , m = 1, . . . , N.
(3.43)
Im
The nonlinear discrete problems (3.39)–(3.40) and (3.42)–(3.43) are solved at each time level by a Oseen-like iterations. See, e.g., [13], [12].
12
ACCEPTED MANUSCRIPT
M(t) Uµ a a
L(t)
kaa
EA T
H
CR IP T
kHH
Figure 1: Elastically supported airfoil with two degrees of freedom.
4. Flow-induced airfoil vibrations
AN US
We consider an elastically supported airfoil ΓWt with two degrees of freedom: the vertical displacement H (positively oriented downwards) and the angle α of rotation around an elastic axis EA (positively oriented clockwise). See Figure 1.
The motion of the airfoil for large vibration amplitudes is described by the system of ordinary differen-
M
tial equations for unknowns H and α:
(4.44)
S α H¨ cos α + Iα α¨ + kαα α = M(t),
(4.45)
PT
ED
mH¨ + S α α¨ cos α − S α α˙ 2 sin α + kHH H = −L(t),
where the dot and two dots denote the first-order and the second-order time derivative, respectively. This
CE
system is derived from the Lagrange equations of the second art, see [28]. We use the following notation: L(t) - aerodynamic lift force, M(t) - aerodynamic torsional moment, m - airfoil mass, S α - static moment of the airfoil around the elastic axis EA, Iα - inertia moment of the airfoil around the elastic axis EA, kHH -
AC
bending stiffness, kαα - torsional stiffness. ˙ System (4.44) – (4.45) is equipped with the initial conditions prescribing H(0), α(0), H(0), α(0). ˙ The
aerodynamic lift force L and the torsional moment M are defined by
L=
−lρre f v2re f c
Z
2 X
τ2 j n j dS ,
M=
ΓWt j=1
lρre f v2re f c2
Z
ΓWt i, j=1
where
r1ort = −(x2 − xEA2 ),
r2ort = x1 − xEA1 , 13
2 X
τi j n j riort dS ,
(4.46)
ACCEPTED MANUSCRIPT
l is the airfoil depth and c is the length of the airfoil cord. By τi j we denote the components of the aerodynamic stress tensor, defined by (2.6) and (2.7) (taking into account that ∇ · v = 0 in the case of the incompressible flow), n = (n1 , n2 ) is the unit outer norm to ∂Ωt on ΓWt (pointing into the airfoil) and xEA = (xEA1 , xEA2 ) is the position of the elastic axis at time t.
CR IP T
System (4.44) – (4.45) is transformed to a first-order system and solved by the Runge-Kutta method. The construction of the ALE mapping, the domain velocity and the discretization in the framework of the ALE method has been carried out according to the papers [28] for incompressible case and [5] for compressible flow.
AN US
5. Numerical experiments 5.1. Flow in a simplified channel of human vocal tract
In numerical experiments we consider a channel modeling a simplified subglottal and supraglottal spaces constricted by two bumps of the length 15.4 mm in the glottis region, which represent a simplified model of vocal folds. The numerical results are obtained for a triangulation consisting of 7790 elements
M
(see Figure 2). We prescribe the following dimensional data: the inlet density ρ¯ in = 1.225 kg m−3 , the outlet atmospheric pressure p¯ out = 97611 Pa, the inlet velocity v¯ in = 4 m s−1 (Re = 4356, inlet Mach number
ED
M = 0.012) or v¯ in = 0.4 m s−1 (Re = 436, M = 0.0012), the fluid viscosity µ = 18 · 10−6 kg m−1 s−1 , the heat conduction coefficient k = 2.428 · 10−2 kg m s−2 K−1 , the specific heat at constant pressure
PT
c p = 1009.992 m2 s−2 K−1 and the Prandtl number Pr = 0.748. The Reynolds number was computed by (2.2), with Lre f = 0.016 m, ρre f = ρ¯ in and vre f = v¯ in .
CE
We use the time step ∆t = 0.00002 s and set T = 0.2 s. For computing the compressible viscous flow we use the NIPG version of the discontinuous Galerkin method with constant CW = 500 inside of the
AC
computational domain and CW = 5000 on the boundary. We compare the distribution of the pressure drop and velocity along the symmetry axis of the channel
for compressible and incompressible flow at time instants t = 0.02, 0.08 and 0.2 s. Results for v¯ in = 4 m s−1
are presented in Figures 3-6 and results for v¯ in = 0.4 m s−1 are shown in Figures 7-10. Here the notation x = x1 is used.
14
CR IP T
ACCEPTED MANUSCRIPT
Figure 2: Triangulation of the computational domain with boundary parts and sizes of the individual parts in mm.
p [Pa]
p [Pa]
1500
1500 1000 500 0 -500 -1000 -1500 -2000 -2500 -3000
500 0 -500 -1000 -1500 -3
-2
-1
0
1
2
3
4
5
x / Lref
AN US
1000
6
-3
-2
-1
0
1 2 x / Lref
3
4
5
6
ED
M
Figure 3: Distribution of the pressure drop for vin = 4 m s−1 along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
|u| [m/s] 70
60
PT
60
|u| [m/s]
50 40
CE
30
50 40 30 20
20
10
10 0
AC
-3
-2
-1
0
0 1
2
3
4
5
6
-3
x / Lref
-2
-1
0
1 2 x / Lref
3
4
5
6
Figure 4: Distribution of the magnitude of the velocity for vin = 4 m s−1 along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
15
ACCEPTED MANUSCRIPT
AN US
CR IP T
1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 61
AC
CE
PT
ED
M
Figure 5: Magnitude of the velocity for vin = 4 m s−1 in the case of incompressible flow at time instants t = 0.02, 0.08 and 0.2 s.
Figure 6: Magnitude of the velocity for vin = 4 m s−1 in the case of compressible flow at time instants t = 0.02, 0.08 and 0.2 s.
16
p [Pa] 10 8 6 4 2 0 -2 -4 -6 -8 -10
p [Pa] 15 10 5 0 -5 -10 -15 -20 -25 -3
-2
-1
0
1 2 x / Lref
3
4
5
6
-3
-2
CR IP T
ACCEPTED MANUSCRIPT
-1
0
1 2 x / Lref
3
4
5
6
AN US
Figure 7: Distribution of the pressure drop for vin = 0.4 m s−1 , along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
|u| [m/s]
|u| [m/s]
4.5 4 3.5 3 2.5 2 1.5 1 0.5 0
7 6 5 4
M
3
-2
-1
0
1 2 x / Lref
3
4
5
ED
-3
6
2 1 0 -3
-2
-1
0
1 2 x / Lref
3
4
5
6
AC
CE
PT
Figure 8: Distribution of the magnitude of the velocity for vin = 0.4 m s−1 , along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
17
ACCEPTED MANUSCRIPT
AN US
CR IP T
0.1 0.4 0.7 1 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7 4 4.3
AC
CE
PT
ED
M
Figure 9: Magnitude of the velocity for vin = 0.4 m s−1 in the case of incompressible flow at time instants t = 0.02, 0.08 and 0.2 s.
Figure 10: Magnitude of the velocity for vin = 0.4 m s−1 in the case of compressible flow at time instants t = 0.02, 0.08 and 0.2 s.
From the above figures we see that the pressure drop between the inlet and outlet computed for incompressible flow is in a good correspondence with the results for compressible flow. The same is true for the velocity magnitude in the narrowest part of the channel. Both the velocity and pressure distributions are identical from the inlet up to the narrowest part of the channel. It is important that the vortex structures are similar for the incompressible and compressible flow. The differences may be caused by different jet declination behind the channel constriction and by the Coanda effect, which is moreover influenced by the 18
ACCEPTED MANUSCRIPT
interaction of the jet with the large-scale supraglottal eddies. 5.2. Flow around the NACA0012 profile We consider laminar flow around the NACA0012 profile for six different flow regimes characterized by
and the Reynolds number Re:
α = 1, Re = 160000,
(C2)
α = 3, Re = 160000,
(C3)
α = 5, Re = 160000,
(C4) (C5) (C6)
AN US
(C1)
CR IP T
the far-field velocity v¯ in = 8 m s−1 (the corresponding inlet Mach number is M = 0.024), angle of attack α
α = 1, Re = 5000,
α = 2, Re = 5000,
α = 3, Re = 5000.
ED
M
We evaluate the aerodynamic coefficients drag cD and lift cL . They are defined as
(cD , cL ) = 2
Z
τn dS ,
(5.47)
ΓW
PT
where ΓW is the surface of the profile, τ is the dimensionless stress tensor and n = (n1 , n2 ) is the unit outer normal to ∂Ω on ΓW .
AC
11).
CE
For each flow regime we carried out computations on the grid consisting of 17158 elements (see Figure
19
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 11: The computational grid around the NACA0012 profile.
M
We prescribe the following data: ρ¯ in = 1.225 kg m−3 , p¯ out = 97611 Pa, k = 2.428 · 10−2 kg m s−2 K−1 , c p = 1009.992 m2 s−2 K−1 and two values of the fluid viscosity µ = 1.8375 · 10−5 kg m−1 s−1 (corresponding
ED
to Re = 160000, Pr = 0.764) and µ = 5.88 · 10−4 kg m−1 s−1 (corresponding to Re = 5000, Pr = 24.459). Moreover ∆t, T and values of CW are chosen similarly as in Section 5.1. For computing the compressible viscous flow we again use the NIPG version of the discontinuous Galerkin method.
PT
For flow regimes (C1)–(C3) we observe small fluctuations of the drag cD and lift cL in dependence on time, as can be seen from Figures 12–17 both for incompressible and compressible cases. After a transient
CE
regime (0 < t < 0.5s) numerical results oscillate around mean values. Table 1 shows the values of the corresponding drag and lift coefficients and in case Re = 160000 also
AC
experimental data from [26]. For the case Re = 5000 we have not found any experimental data. For flow regimes (C1)–(C3) we write the values in the form ”mean value ± standard deviation”. For flow regimes
(C4)–(C6) we obtained a stationary flow and present the steady-state values of the drag and lift coefficients. In Table 1 it is possible to see an interesting fact that the experimental values of drag and lift are in the
intervals between the mean values numerically computed with the use of incompressible and compressible models. It is also interesting that in all cases ccompr. < cincompr. and ccompr. > cincompr. . D D L L
20
ACCEPTED MANUSCRIPT
ccompr. D
cincompr. D
exp. data
ccompr. L
cincompr. L
exp. data
(C1) (C2) (C3) (C4) (C5) (C6)
0.0096 ± [0.001] 0.0104 ± [0.002] 0.0136 ± [0.004] 0.0507 0.0517 0.0531
0.0116 ± [0.001] 0.0132 ± [0.002] 0.0169 ± [0.004] 0.0533 0.0541 0.0566
0.0104 0.0114 0.0140 -
0.155 ± [0.041] 0.382 ± [0.027] 0.562 ± [0.036] 0.0374 0.0727 0.0988
0.086 ± [0.046] 0.30 ± [0.043] 0.498 ± [0.047] 0.0366 0.0696 0.1062
0.11 0.33 0.55 -
0.02
0.02
0.018
0.018
0.016
0.016
0.014
0.014 CD
0.012 0.01
CR IP T
Flow
AN US
CD
Table 1: Values of the computed aerodynamic coefficients cD and cL for incompressible and compressible viscous flow and experimental data.
0.012 0.01
0.008
0.008
0.006
0.006
0.5
1
1.5
2
2.5
(a) Incompressible flow
0.5
1
1.5
2
2.5
t[s]
(b) Compressible flow
M
t[s]
ED
Figure 12: Values of the computed aerodynamic coefficients cD for the incompressible flow (left) and compressible flow (right) in flow regime (C1) in dependence on time t[s].
0.4
0.2
CE
0.15
0.3 0.25 CL
0.3 0.25 CL
0.4 0.35
PT
0.35
0.2 0.15
0.1
0.1
0.05
0.05
0
AC
0.5
1
0 1.5
2
2.5
0.5
t[s]
1
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
Figure 13: Values of the computed aerodynamic coefficients cL for the incompressible flow (left) and compressible flow (right) in flow regime (C1) in dependence on time t[s].
21
0.025
0.025
0.02
0.02
0.015
0.01
0.015
0.01
0.005
0.005 0.5
1
1.5
2
2.5
0.5
t[s]
CR IP T
CD
CD
ACCEPTED MANUSCRIPT
1
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
AN US
Figure 14: Values of the computed aerodynamic coefficients cD for the incompressible flow (left) and compressible flow (right) in flow regime (C2) in dependence on time t[s].
0.5
0.5
0.45
0.45
0.4
0.4
CL
0.35
0.3
M
CL
0.35
0.25 0.2 0.15
0.3
0.25 0.2 0.15
0.1
0.1 1
1.5
2
ED
0.5
2.5
0.5
1
t[s]
1.5
2
2.5
t[s]
(b) Compressible flow
PT
(a) Incompressible flow
CE
Figure 15: Values of the computed aerodynamic coefficients cL for the incompressible flow (left) and compressible flow (right) in flow regime (C2) in dependence on time t[s].
0.04
AC
0.04
0.035
0.03
0.03
0.025
0.025 CD
CD
0.035
0.02
0.015
0.02 0.015
0.01
0.01
0.005
0.005
0
0 0.5
1
1.5
2
2.5
0.5
t[s]
1
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
Figure 16: Values of the computed aerodynamic coefficients cD for the incompressible flow (left) and compressible flow (right) in flow regime (C3) in dependence on time t[s].
22
ACCEPTED MANUSCRIPT
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3 0.5
1
1.5
2
2.5
0.5
t[s]
1
CR IP T
CL
0.8
CL
0.8
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
AN US
Figure 17: Values of the computed aerodynamic coefficients cL for the incompressible flow (left) and compressible flow (right) in flow regime (C3) in dependence on time t[s].
5.3. Numerical experiments of the flow-induced airfoil vibrations
Following the example of numerical simulations of flow-induced vibrations published in [28], we prescribe the input data as follows: m = 0.086622 kg, S α = −0.000779673 kg m, Iα = 0.000487291 kg m−2 ,
M
kHH = 105.109 N m−1 , kαα = 3.695582 N m rad−1 , the fluid viscosity µ = 1.8375 · 10−5 kg m−1 s−1 , the density ρ¯ in = 1.225 kg m−3 , the pressure p¯ out = 101250 Pa, k = 2.428 · 10−2 kg m s−2 K−1 , c p =
ED
1009.992 m2 s−2 K−1 , Pr = 0.764. The elastic axis EA is located at 40% of the profile cord.The airfoil cord c = 0.3 m, airfoil depth l = 0.05 m, and the following initial conditions for the structural equations
PT
˙ H(0) = −20 mm, α(0) = 6, H(0) = α(0) ˙ = 0. We prescribed inlet velocities v¯ in = 10 m s−1 (corresponding to Re = 200000 and Mach number M = 0.0294) and v¯ in = 20 m s−1 (corresponding to Re = 400000 and
CE
Mach number M = 0.0588). The computation started at a time instant t = −δ < 0 with a fixed airfoil. Then at time t = 0 the airfoil was released and the fluid-structure interaction process followed.
AC
Incompressible flow was computed by the method described in section 3.1 with the use of its ALE version. For compressible flow we used the SIPG STDG method with quadratic polynomials in space (r = 2) and linear polynomials in time (q = 1) (the results are denoted by ST-DGr2q1). The parameter CW was chosen similarly as in section 5.1. Figures 18–21 show the comparison of the displacement H and the angle α computed for incompressible and compressible flow. The time responses of the airfoil are damped due to the unsteady aerodynamic forces given by (4.46) because the considered inlet flow velocities are far below the critical flow velocity 42 m s−1 for flutter, see [28]. The responses computed by DGFEM for the compressible flow are in good agreement with the results computed by FEM for the incompressible flow; the maximal displacements for t > 0 are very similar for 23
ACCEPTED MANUSCRIPT
both the translation H and the rotation α, and the differences in frequencies of oscillations are less than 1%
8
15
6
10
4
5
CR IP T
20
2
0
α
H [mm]
for the inlet flow velocity vin = 10 m s−1 and less then 4% for vin = 20 m s−1 .
0
-5 -10
-2
-15
-4
-20 0
0.5
1
1.5
-6
2
0
0.5
t[s]
1
1.5
2
AN US
t[s]
Figure 18: Values of the computed displacement H and angle α for the compressible (black) and the incompressible flow (green) for v¯ in = 10 m s−1 in dependence on time t[s].
8
15
6
M
0 -5 -10 -15 -20 0
0.2
0.4
ED
H [mm]
5
0.6
0.8
4 2
α
10
0 -2 -4 -6
1
0
t[s]
0.2
0.4
0.6
0.8
1
PT
t[s]
CE
Figure 19: Values of the computed displacement H and angle α for the compressible (black) and the incompressible flow (green) for v¯ in = 20 m s−1 in dependence on time t[s].
AC
6. Conclusion
The paper is devoted to the analysis of the discontinuous Galerkin solution of compressible low Mach
number flows. The goal is to prove the robustness of the developed DGM with respect to the Mach number at incompressible limit. To this end, several test problems were solved by the DGM and results were compared with their numerical simulation with the aid of the stabilized finite element method applied to incompressible Navier-Stokes equations. The performed numerical tests show that the correspondence of compressible and incompressible results is good not only for the stationary cases but also for the case of time dependent flow. The time dependent character of the solution is caused either by moving vortices or by interaction of fluid with an elastically supported vibrating airfoil. In all these cases the DGM provides results comparable to the solution of the incompressible model. 24
ACCEPTED MANUSCRIPT
As a conclusion, it is possible to say that the described numerical method for the solution of the compressible Navier-Stokes equations allows the simulation of very low Mach number flows and yields the results comparable with results obtained with the aid of the incompressible model, both in internal aerodynamics, for the airflow in a constricted channel (M = 0.0012), and in outer aerodynamics for the fixed
CR IP T
airfoil (M = 0.024) as well as for the flow induced airfoil vibrations (M = 0.0294). The results prove the robustness of the discontinuous Galerkin method with respect to the low Mach number. Acknowledgment
The work of M. Bal´azsov´a was supported by the Grant No. SVV-2017-260455 of the Charles University. The
AN US
research of M. Feistauer was supported by the Grant No. 17-01747S of the Czech Science Foundation and finally the research of P. Sv´acˇ ek and J. Hor´acˇ ek was supported by the Grant No. 16-01246S of the Czech Science Foundation.
References
[1] Bassi, F., Rebay, S.: High-order accurate discontinuous finite element solution of the 2D Euler equations. J. Comput. Phys. 138 (1997), 251–285.
Methods Fluids 31 (1999), 79–95.
M
[2] Baumann, C. E., Oden, J. T.: A discontinuous hp finite element method for the Euler and Navier-Stokes equations. Int. J. Numer.
ED
[3] Bruneau, Ch.-H., Fabrie, P.: Effective downstream boundary conditions for incompressible Navier–Stokes equations. Int. J. Nu-
mer. Methods Fluids 19(8) (1994), 693–705.
ˇ [4] Cesenek, J., Feistauer, M., Hor´acˇ ek, J., Kuˇcera, V., Prokopov´a, J.: Simulation of compressible viscous flow in time-dependent
PT
domains. Appl. Math. Comput. 13 (2013), 7139–7150. ˇ [5] Cesenek, J., Feistauer, M., Kos´ık, A.: DGFEM for the analysis of airfoil vibrations induced by compressible flow. ZAMM 93(13) (2013), 387–402.
CE
[6] Cockburn, B., Karniadakis, G.E., Shu C.W., Eds.: Discontinuous Galerkin Methods, Lecture Notes in Computational Science and Engineering 11, Springer, Berlin (2000). [7] Darwish, M., Moukalled, F., Sekar, B.: A robust multi-grid pressure-based algorithm for multi-fluid flow at all speeds, Int. J.
AC
Numer. Methods Fluids 41 (2003) 1221–1251.
[8] Davis, T. A., Duff, I. S.: A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Trans. Math. Softw. 25 (1999) 1–19.
[9] Dolejˇs´ı, V.: Semi-implicit interior penalty discontinuous Galerkin methods for viscous compressible flows. Commun. Comput.
Phys. 4 (2008), 231–274. [10] Dolejˇs´ı, V., Feistauer, M.: Discontinuous Galerkin Method - Analysis and Applications to Compressible Flow. Springer, (2015). [11] Dumbser, M.: Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier-Stokes equations.
Comput. Fluids 39(1) (2010), 60–76. [12] Feistauer, M., Dolejˇs´ı, V., Kuˇcera, V.: On the discontinuous Galerkin method for the simulation of compressible flow with wide range of Mach numbers. Comput. Vis. Sci. 10 (2007), 17–27. [13] Feistauer, M., Kuˇcera, V.: On a robust discontinuous Galerkin technique for the solution of compressible Flow. J. Comput. Phys. 224 (2007), 208–221.
25
ACCEPTED MANUSCRIPT
[14] Feistauer, M., Hasnedlov´a-Prokopov´a, J., Hor´acˇ ek, J., Kos´ık, A., Kuˇcera, V.: DGFEM for dynamical systems describing interaction of compressible fluid and structures. J. Comput. Appl. Math. 254 (2013), 17–30. [15] Gelhard, T., Lube, G., Olshanskii, M. A., Starcke, J.-H.: Stabilized finite element schemes with LBB-stable elements for incompressible flows. J. Comput. Appl. Math. 177 (2005), 243–267. [16] Hartmann, R., Houston, P.: Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations, J.
CR IP T
Comput. Phys. 183 (2002), 508–532.
[17] Klein, R.: Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics 1: one-dimensional flow,
J. Comput. Phys. 121 (1995) 213–237.
[18] Meister, A.: Comparison of different Krylov subspace methods embedded in an implicit finite volume scheme for the computation of viscous and inviscid flow fields on unstructured grids, J. Comput. Phys. 140 (1998) 311–345.
[19] Meister, A.: Viscous flow fields at all speeds: Analysis and numerical simulation, J. Appl. Math. and Physics 54 (2003) 1010–
AN US
1049.
[20] Meister, A., Struckmeier, J.: Hyperbolic Partial Differential Equations, Theory, Numerics and Applications, Vieweg, 2002. [21] Nerinckx, K., Vierendeels, J., Dick, E.: Mach-uniformity through the coupled pressure and temperature correction algorithm. J.
Comput. Physics 206 (2005), 597-623.
[22] Nomura, T., Hughes, T.J.R.: An arbitrary Lagrangian-Eulerian finite element method for interaction of fluid and a rigid body,
Comput. Methods Appl. Mech. Engrg. 95 (1992), 115–138.
[23] Park, J.H., Munz, C.-D.: Multiple pressure variables methods for fluid flow at all Mach numbers, Int. J. Numer. Methods Fluids
M
49 (2005) 905–931.
[24] Punˇcoch´aˇrov´a-Poˇr´ızkov´a, P., F¨urst, J., Hor´acˇ ek, J., Kozel, K.: Numerical solutions of unsteady flows with low inlet Mach
ED
numbers. Math. Comput. Simul. 80(8) (2010), 1795–1805.
[25] Roller, S., Munz, C.-D., Geratz, K. J., Klein, R.: The multiple pressure variables method for weakly compressible fluids, Z. Angew. Math. Mech. 77 (1997) 481–484.
PT
[26] Sheldahl, R. E., Klimas, P. C.: Aerodynamic Characteristics of Seven Airfoil Sections Through 180 Degrees Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbines, SAND80-2114, March 1981, Sandia National Laboratories, Albuquerque, New Mexico.
CE
[27] Seo, J. H., Mittal, R.: A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries, J. Comput. Phys. 77 (1997) 481–484. [28] Sv´acˇ ek, P., Feistauer, M., Hor´acˇ ek, J.: Numerical simulation of flow induced airfoil vibrations with large amplitudes. J. of Fluids
AC
Struct. 23 (2007), 391–411. [29] van der Vegt, J. J. W., van der Ven, H.: Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flow. J. Comput. Phys. 182 (2002), 546–585.
[30] van der Heul, D.R., Vuik, C., Wesseling, P.: A conservative pressure-correction method for flow at all speeds, Comput. Fluids 32 (2003) 1113–1132. [31] Wenneker, I., Segal, A., Wesseling, P.: A Mach-uniform unstructured steggered grid method, Int. J. Numer. Methods Fluids 40 (2002) 1209–1235. [32] Wesseling, P.: Principles of Computational Fluid Dynamics, Springer, Berlin, 2001. [33] Wesseling, P., van der Heul, D.R.: Uniformly effective numerical methods for hyperbolic systems, Computing 66 (2001) 249– 267.
26
Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow ˇ Monika Balazsov a, Miloslav Feistauer, Petr Svaˇ ´ ´ Jan Cesenek, ´ cek, Jarom´ır Horaˇ ´ cek PII: DOI: Reference:
S0045-7930(18)30391-8 10.1016/j.compfluid.2018.07.002 CAF 3943
To appear in:
Computers and Fluids
Received date: Revised date: Accepted date:
22 September 2017 28 June 2018 2 July 2018
ˇ Please cite this article as: Monika Balazsov a, Miloslav Feistauer, Petr Svaˇ ´ ´ Jan Cesenek, ´ cek, Jarom´ır Horaˇ ´ cek, Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow, Computers and Fluids (2018), doi: 10.1016/j.compfluid.2018.07.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • viscous incompressible and viscous compressible flow for low Mach numbers
• flow in a simplified channel of the human vocal tract • flow around an airfoil
AC
CE
PT
ED
M
AN US
• flow-induced airfoil vibrations
CR IP T
• robustness of the discontinuous Galerkin method with respect to low Mach numbers
1
ACCEPTED MANUSCRIPT
Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow✩
CR IP T
2 ˇ Monika Bal´azsov´a1 , Jan Cesenek , Miloslav Feistauer1 , Petr Sv´acˇ ek3 , Jarom´ır Hor´acˇ ek4
Abstract
The novelty of the paper is the application of the discontinuous Galerkin method (DGM) to the solution of viscous compressible flow with low Mach number (M ≈ 10−2 − 10−3 ). The compressible Navier-Stokes
AN US
equations are written in the conservative form using conservative variables. They are discretized by the DGM in space. In time we use either the backward difference formula (BDF) or the discontinuous Galerkin method. The resulting technique is called the space-time discontinuous Galerkin method (STDGM). The goal of the presented analysis is to prove the ability of the developed DGM to solve viscous compressible low Mach number flows described by the compressible Navier-Stokes equations without any
M
modification and to obtain results comparable with simulations of incompressible viscous flow described by the incompressible Navier-Stokes equations, solved by the finite element method (FEM). We present
ED
numerical results for a flow in a simplified channel of the human vocal tract, for a flow around the airfoil NACA0012 and moreover for flow-induced airfoil vibrations. The results prove the robustness of the discontinuous Galerkin method with respect to the low Mach number.
PT
Keywords: incompressible Navier-Stokes equations, compressible Navier-Stokes equations, low Mach number flow, finite element method, discontinuous Galerkin method, fluid-structure interaction
CE
2000 MSC: 65N50, 65N15, 65D05
AC
1. Introduction
In the numerical solution of compressible flow, it is necessary to overcome several issues: resolve accu-
rately shock waves, contact discontinuities and boundary layers. All these phenomena are connected with ✩ The work of M. Bal´ azsov´a was supported by the Grant No. SVV-2017-260455 of the Charles University. The research of M. Feistauer was supported by the Grant No. 17-01747S of the Czech Science Foundation and finally the research of P. Sv´acˇ ek and J. Hor´acˇ ek was supported by the Grant No. 16-01246S of the Czech Science Foundation. ˇ Email addresses: [email protected] (Monika Bal´azsov´a), [email protected] (Jan Cesenek), [email protected] (Miloslav Feistauer), [email protected] (Petr Sv´acˇ ek), [email protected] (Jarom´ır Hor´acˇ ek) 1 Charles University Prague, Faculty of Mathematics and Physics, Sokolovsk´ a 83, 186 75 Praha, Czech Republic 2 Aeronautical Research and Test Institute, Beranov´ ych 130, 199 05 Prague 9, Czech Republic 3 Department of Technical Mathematics, Faculty of Mechanical Engineering, Czech Technical University, Karlovo n´ amˇest´ı 13, 121 35 Prague 2, Czech Republic 4 Institute of Termomechanics, Academy of Sciences of the Czech Republic, Dolejˇskova 5, 182 00 Prague 8, Czech Republic
Preprint submitted to Computers & Fluids
June 28, 2018
ACCEPTED MANUSCRIPT
the simulation of high speed flow with high Mach numbers. However, it appears that the solution of low Mach number flow is also rather difficult. This is caused by the stiff behaviour of numerical schemes and acoustic phenomena appearing in low Mach number flows at incompressible limit. This led to the development of special finite volume techniques allowing the simulation of compressible flow at incompressible
CR IP T
limit, which are based on modifications of Euler equations describing inviscid compressible flow. From a wide literature we mention, for example, works [17], [25], [20], Chapter 5, or [32], Chapter 14). However, these techniques could not be applied to the solution of high speed flow. Therefore, further attempts were concentrated on the extension of these methods to the solution of flows at all speeds. A success in this direction was achieved by several authors. Let us mention, for example, the works by Wesseling et al. (e. g. [30],
AN US
[31], [33]), Parker and Munz ([23]), Meister ([18], [19]), Nerinckx et al. ([21]) and Darwish et al. ([7]). Main ingredients of these techniques are finite volume schemes applied on staggered grids, combined with multigrid, the use of the pressure-correction, multiple pressure variables and flux preconditioning. The paper [27] uses the immersed boundary method. There is a question, if it is possible to develop a numerical method allowing the solution of compressible flow with a wide range of Mach numbers without modifi-
M
cation of the governing equations. This was successful in papers [13] and [12] dealing with the solution of the inviscid compressible Euler equations with a wide range of Mach numbers. This technique is based
ED
on the discontinuous Galerkin method (DGM), which combines advantages of the finite volume as well as finite element methods. It employs piecewise polynomial approximations without any requirement on the continuity on interfaces between neighbouring elements.
PT
The DGM was first applied to the numerical simulation of the compressible Euler equations, by Bassi and Rebay in [1], where the space discontinuous Galerkin discretization was combined with explicit Runge-
CE
Kutta time discretization. In [2] Baumann and Oden apply an hp version for the solution of the space DG discretization with explicit time stepping to compressible flow. Van der Vegt and van der Ven use the space-
AC
time discontinuous Galerkin method for the solution of the Euler equations in [29], where the discrete problem is solved with the aid of a multigrid accelerated pseudo-time-integration. The DGM applied to compressible flow was combined with adaptive techniques (see, e.g. [16]). An extension of results from [2] is the subject of [9]. In [11] the DGM was combined with the PN P M method. It is possible to say that
during the last decades, a number of papers devoted to the solution of compressible flow by the DGM was published. An older survey can be found, for example, in [6], whereas recent references are mentioned in [10]. The subject of this paper is the numerical simulation of viscous low Mach number compressible flow described by the compressible Navier-Stokes equations and compared with results obtained with the use of 2
ACCEPTED MANUSCRIPT
incompressible Navier-Stokes equations. The novelty of the presented research lies in the demonstration of the applicability of the developed DG method for the numerical solution of compressible low Mach number flow, practically at incompressible limit. The robustness of the DGM with respect to the Mach number is proved by the comparison of results with the solution of incompressible Navier-Stokes equations and in
CR IP T
several cases with Sandia experimental data.
In Section 2 we formulate the continuous problem of incompressible as well as compressible viscous flow. Section 3 is devoted to the discretization of the problem. For the discretization of the incompressible flow we use standard conforming stabilized finite element method (FEM). The discretization of compressible flow is carried out by the discontinuous Galerkin method in space. In both cases, for the time dis-
AN US
cretization we use the backward difference formula (BDF). In the case of the compressible Navier-Stokes equations, also discontinuous Galerkin in time is used. Then we get the so-called space-time discontinuous Galerkin method (STDGM). Section 4 presents the description of the simulation of flow induced airfoil vibrations. Here the above problems are reformulated in time-dependent domains with the aid of the arbitrary Lagrangian-Eulerian (ALE) method. In Section 5 we present several numerical experiments showing
M
the comparison of incompressible flow and compressible viscous flow with low Mach numbers. In the first case we compare incompressible flow and compressible viscous flow in a fixed 2D channel, which repre-
ED
sents a simplified model of the glottal region of the human vocal tract, i.e. airways in the human subglottal and supraglottal spaces with a sudden constriction at the glottis. Our choice of the domain was motivated by number of works on fluid-structure interaction problems for vocal folds, e.g. [14], [4] and [24]. In the
PT
second case we investigate the incompressible flow and the compressible viscous flow around the fixed airfoil NACA0012 for different angles of attack. Moreover, we are concerned with the simulation of airfoil
CE
vibrations induced by incompressible flow and low Mach number compressible flow.
AC
2. Governing equations The aim of this paper is to show the ability of the DGM to approximate compressible flow for low
Mach numbers, particularly in the case when the flow is almost incompressible. To this end, numerical results for the compressible flow for low Mach numbers are compared with the numerical results for the incompressible flow. In what follows, the incompressible and compressible flow models are formulated.
3
ACCEPTED MANUSCRIPT
2.1. Incompressible flow The incompressible viscous flow in a computational polygonal domain Ω ⊂ R2 and a time interval [0, T ], T > 0 is governed by the Navier-Stokes equations in the dimensionless form
∇ · v = 0,
(2.1)
CR IP T
∂v 1 + (v · ∇)v + ∇p − △v = 0, ∂t Re
where v = (v1 , v2 ) = v(x, t) denotes the velocity vector, p = p(x, t) denotes the pressure, Re is the Reynolds number, defined as Lre f vre f ρre f , µ
(2.2)
AN US
Re =
where Lre f is the reference length, vre f is the reference velocity, ρre f is the reference density (for incompressible flow it is the constant fluid density) and µ is the dynamic viscosity of the fluid. System (2.1) is equipped with the initial condition
M
v(x, 0) = v0 (x)
x ∈ Ω,
(2.3)
ED
and with the boundary conditions
a) v = vin = (vin1 , vin2 )T
on ΓI ,
b) v = 0
PT
1 ∂v 1 =0 c) − (p − pout )n − (v · n)− v + 2 Re ∂n
on ΓW ,
(2.4)
on ΓO ,
CE
where pout is a prescribed pressure distribution on the outlet and n denotes the unit outer normal vector to the boundary ∂Ω = ΓI ∪ ΓO ∪ ΓW , where ΓI denotes the inlet, ΓO the outlet and ΓW the fixed impermeable
AC
walls. The boundary condition (2.4c) is a modification of the so-called ’do-nothing’ boundary condition, cf. [3]. We set λ− = min{0, λ} and λ+ = max{0, λ} for λ ∈ R. 2.2. Compressible flow The compressible viscous flow in the domain Ω ⊂ R2 and the time interval [0, T ] is described by the system of the Navier-Stokes equations in the dimensionless form ∂w X ∂ f s (w) + ∂t ∂x s s=1 2
=
2 X ∂R s (w, ∇w) s=1
4
∂x s
,
(2.5)
ACCEPTED MANUSCRIPT
where w is the so called state vector defined as w = (w1 , w2 , w3 , w4 )T = (ρ, ρv1 , ρv2 , E)T ∈ R4 and w0 (x). Functions
CR IP T
f s (w) = (ρv s , ρv1 v s + δ1s p, ρv2 v s + δ2s p, (E + p)v s )T , !T γ ∂θ V V V V R s (w, ∇w) = 0, τ s1 , τ s2 , τ s1 v1 + τ s2 v2 + Re Pr ∂x s
represent the inviscid and viscous terms, respectively. We use the following notation: ρ - density, δi, j Kronecker symbol, E - total energy, θ - absolute temperature, τV is the viscous part of the stress tensor with components
! 1 ∂vi ∂v j 2 + − δi j ∇ · v , = Re ∂x j ∂xi 3
i, j = 1, 2
AN US
τVij
(2.6)
and τ denotes the stress tensor with components
τi j = −pδi j + τVij ,
i, j = 1, 2.
(2.7)
by (2.2), Pr =
cpµ k
M
Further, we use the following notation: γ > 1 - Poisson adiabatic constant, Re - Reynolds number defined - Prandtl number, where c p is the specific heat at constant pressure and k is the heat
ED
conductivity.
The above system is completed by the thermodynamical relations
PT
! 1 2 p = (γ − 1) E − ρ|v| , 2
θ=
E 1 2 − |v| . ρ 2
(2.8)
CE
The resulting system is equipped with the initial condition
x ∈ Ω0 ,
(2.9)
AC
w(x, 0) = w0 (x),
and boundary conditions
a) ρ|ΓI = ρin , b) v|ΓI = vin , 2 X ∂θ = 0 on ΓI , c) τVij ni v j + k ∂n i, j=1 d) v = 0, f)
2 X
e)
∂θ =0 ∂n
τVij ni = 0,
on ΓW ,
j = 1, 2,
i=1
5
(2.10)
g)
∂θ = 0 on ΓO . ∂n
ACCEPTED MANUSCRIPT
Since we consider low Mach number flow, we also consider the following condition
p = pout
on ΓO ,
(2.11)
CR IP T
see also [10], Sections 9.1.2 and 8.3. 2.3. Flow in a time-dependent domain Ωt
In this case the impermeable wall ΓWt of the domain occupied by the fluid may depend on time t ∈ [0, T ] and the boundary conditions (2.4), b) and (2.10), d) are replaced by the condition
AN US
v|ΓWt = zD = (zD1 , zD2 ) − velocity of the moving wall ΓWt .
(2.12)
In order to treat the time dependence of the domain, we use the so-called arbitrary Lagrangian-Eulerian (ALE) technique, proposed in [22]. We define a reference domain Ω0 and introduce a regular one-to-one
M
ALE mapping of Ω0 onto Ωt (cf. [22], [28] and [29])
ED
At : Ω0 −→ Ωt , i.e. X ∈ Ω0 7−→ x = x(X, t) = At (X) ∈ Ωt . Here we use the notation X for points in Ω0 and x for points in Ωt .
CE
PT
Further, we define the domain velocity:
˜z(X, t) =
∂ At (X), ∂t
z(x, t) = ˜z(At−1 (x), t),
t ∈ [0, T ], X ∈ Ω0 , t ∈ [0, T ], x ∈ Ωt
(2.13) (2.14)
AC
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
(2.15)
where f˜(X, t) = f (At (X), t), X ∈ Ω0 , x = At (X). As a direct consequence of the chain rule we get the relation ∂f DA f = + div (z f ) − f div z. Dt ∂t
6
(2.16)
ACCEPTED MANUSCRIPT
This leads to the ALE formulation of the incompressible Navier-Stokes equations in the form DA v 1 + ((v − z) · ∇)v + ∇p − △v = 0. Dt Re
(2.17)
CR IP T
The compressible Navier-Stokes equations have the ALE form X ∂R s (w, ∇w) DA w X ∂ g s (w) + + w divz = , Dt ∂x s ∂x s s=1 s=1 2
2
where s = 1, 2,
(2.19)
AN US
g s (w) := f s (w) − z s w,
(2.18)
are the ALE modified inviscid fluxes.
3. Discretization
The space semidiscretization of the above initial boundary-value problems is carried out on a partition
M
Th of the domain Ω, formed by a finite number of closed triangles with disjoint interiors. We assume that the intersection of two elements K, K ′ ∈ Th is either empty or a common edge or a common vertex of these
ED
elements and the end points of ΓI , ΓO , ΓW are vertices of some elements K ∈ Th . In order to discretize the problem in time, we consider an equidistant partition 0 = t0 < t1 < · · · < tN =
PT
T, tk = k∆t of the time interval [0, T ] with a time step ∆t > 0. 3.1. Incompressible flow
AC
CE
The time derivative in (2.1) is then approximated by the second order backward difference formula 3vn+1 − 4vn + vn−1 ∂v , t=t ≈ ∂t n+1 2∆t
(3.20)
where vk denotes the approximation of v(tk ) for k = 0, . . . , N, and we set v−1 := v0 . Now, let us focus on the space discretization at a given time level tn+1 using the finite element method. The standard conforming finite element spaces defined with the aid of the above mentioned admissible triangulation Th of Ω are used. Particularly, the finite element spaces for the flow velocity Yh and for pressure Qh are defined as i2 h Yh = ϕ = (ϕ1 , ϕ2 ) ∈ C(Ω) : ϕi ∈ P2 (K), ∀K ∈ Th , i = 1, 2 , K
7
ACCEPTED MANUSCRIPT
and
Qh = ϕ ∈ C(Ω) : ϕ ∈ P1 (K), ∀K ∈ Th , K
where by the symbol Pr (K) the space of all polynomials on the triangle K of degree less or equal r is
CR IP T
denoted. Further, the space Wh of test functions is defined as Wh = {ϕ ∈ Yh : ϕ = 0 on ΓI ∪ ΓW } .
The finite element approximation is introduced in the standard way using the multiplication of equations (2.1) by the test functions ϕ ∈ Wh and q ∈ Qh , integration over Ω and using the Green theorem together
AN US
with the boundary condition (2.4c). This leads to the definition of the forms Z
1 1 3 v · ϕ + ((v · ∇)v) · ϕ − ((v · ∇)ϕ) · v − p(∇ · ϕ) dx 2∆t 2 2 Ω Z Z 1 1 (v · n)+ v · ϕ dS + ∇v · ∇ϕ + (∇ · v)q dx, + 2 Re Ω ΓO Z 1 (4vn − vn−1 ) · ϕ dx, f (V) = 2∆t Ω =
(3.21)
(3.22)
M
a(U, V)
for any U = (v, p) ∈ Yh × Qh and V = (ϕ, q) ∈ Wh × Qh . In order to overcome the possible instability of
ED
the convection dominated flow in the case of Galerkin approximations, the residual based stabilization is employed, cf. [15]. It uses the SUPG/PSPG stabilization terms defined by X
PT
L(U, V)
=
CE
F (U; V) =
δK
K∈Th
X
K∈Th
δK
Z
Z
K
K
! 1 3v − △v + (v · ∇)v + ∇p · (v · ∇) ϕ + ∇q dx, 2∆t Re
4vn − vn−1 · (v · ∇)ϕ + ∇q dx, 2∆t
(3.23) (3.24)
AC
in combination with the stabilization term for div-div stabilization given as
P(U, V) =
X
τK
K∈Th
Z
K
(∇ · v)(∇ · ϕ) dx.
(3.25)
The choice of the parameters δK = h2K and τK = 1 is carried out according to [15] or [28], where hK denotes the local element length. Using the above defined non-linear forms, the non-linear stabilized discrete problem at a time instant t = tn+1 reads: Definition An approximate solution of problem (2.1), (2.4) is defined as U = (v, p) ∈ Yh ×Qh , p := pn+1 , 8
ACCEPTED MANUSCRIPT
v := vn+1 , such that v satisfies approximately the conditions (2.4 a-b) and
a(U, V) + L(U, V) + P(U, V) = f (V) + F (U; V),
CR IP T
holds for all V = (ϕ, q) ∈ Wh × Qh .
(3.26)
The solution of the arising system of non-linear equations is realized by the Oseen linearization. The linearized problem is solved with the aid of a direct solver. See, e.g., [28]. 3.2. Compressible flow
For the space discretization we use the discontinuous Galerkin method (DGM) applied on the triangu-
AN US
lation Th , introduced in the previous section. Furthermore, we consider a partition 0 = t0 < t1 . . . < T as above with a time step ∆t. Let us denote by Fh the set of all faces of elements from Th . It consists of the set of all interior and boundary faces: Fh = FhI ∪ FhB . Each face Γ ∈ Fh is associated with a unit normal nΓ , which is an outer normal to ∂Ω for Γ ∈ FhB . The approximate solution is sought in the space of piecewise
M
polynomial functions
Shr = (S hr )4 ,
S hr = {v; v|K ∈ Pr (K) ∀K ∈ Th }.
(3.27)
ED
where
The symbols [u]Γ and huiΓ will denote the jump and the mean value of u ∈ Shr on the face Γ ∈ FhI and [u]Γ = huiΓ = u|Γ for Γ ∈ FhB .
PT
In order to derive the space semidiscrete problem, we assume that there exists an exact solution w ∈ C 1 ([0, T ]; Shr ) of the Navier-Stokes equation (2.5). We multiply (2.5) by a test function ϕ h ∈ Shr , integrate
AC
CE
over an element K ∈ Th , apply Green’s theorem and sum over all elements. Then we can formally write X Z ∂w · ϕ dx + Inv + Vis = 0, K ∂t K∈T
(3.28)
h
where
Inv =
XZ
K∈Th
Vis = −
∂K s=1
XZ
K∈Th
2 X
f s (w)n s · ϕ dS −
2 X
∂K s=1
2 XZ X
K∈Th
R s (w, ∇w)n s · ϕ dS +
K s=1
f s (w) ·
2 XZ X
K∈Th
K s=1
ϕ ∂ϕ dx, ∂x s
R s (w, ∇w) ·
(3.29) ϕ ∂ϕ dx, ∂x s
represent the inviscid and viscous terms and n = (n1 , n2 ) is the outer unit normal to ∂K. To discretize the inviscid term, we use the so-called numerical flux H, which has following properties: 9
ACCEPTED MANUSCRIPT
H(w1 , w2 , n) is locally Lipschitz-continuous with respect to w1 , w2 , P H(w, w, n) is consistent: H(w, w, n) = ds=1 f (w)n s ,
(H1) (H2)
H(w1 , w2 , n) is conservative: H(w1 , w2 , n) = −H(w2 , w1 , −n).
(H3)
Z X 2
f s (w)n s · ϕ dS ≈
Γ s=1
Z
Γ
CR IP T
Integrals over faces Γ ∈ Fh are approximated with the aid of the numerical flux H: (R) H(w(L) Γ , wΓ , nΓ ) ϕ dS .
(3.30)
Then the inviscid terms Inv are discretized by convection form bh (w, ϕ ). Hence
AN US
Inv ≈ bh (w, ϕ ), where XZ
Γ
Γ∈FhI
−
ϕ] dS + nΓ ) [ϕ
XZ
Γ∈FhB
2 XZ X
f s (w) ·
Γ
(L) H(w(L) Γ , wΓ , nΓ ) ϕ dS
(3.31)
ϕ ∂ϕ dx. ∂x s
M
bh (w, ϕ ) =
(R) H(w(L) Γ , wΓ ,
K s=1
K∈Th
ED
For details we can refer to the paper [13].
Now we can turn our attention to the approximation of the viscous term Vis. It is based on the relation
PT
R s (w, ∇w) =
2 X
K s,k (w)
k=1
∂w , ∂xk
s = 1, 2,
(3.32)
CE
where K s,k are 4 × 4 matrices dependent on w. For more details see [10].
AC
ah (w, ϕ ) =
2 XZ X
K∈Th
−
K s,k=1
! + 2 *X 2 XZ X ϕ ∂w ∂w ∂ϕ ϕ] dS (3.33) K s,k (w) K s,k (w) · dx − n s · [ϕ ∂xk ∂x s ∂xk Γ s=1 k=1 I
2 XZ X
Γ∈FhB
−θ
Γ s,k=1
Γ∈Fh
K s,k (w)
2 XZ X
Γ∈FhB
Γ s,k=1
+ 2 *X 2 XZ X ∂w ∂ϕ KTs,k (w) n s · ϕ dS − θ n s · [w] dS ∂xk ∂xk Γ s=1 k=1 I Γ∈Fh
KTs,k (w)
∂ϕ n s · (w − wB ) dS , ∂xk
where θ = 1, 0 or −1 (SIPG, IIPG, resp. NIPG variants) and wB is a boundary state obtained by Dirichlet data and extrapolation. In the scheme we include interior and boundary penalty terms, vanishing for a sufficiently smooth exact
10
ACCEPTED MANUSCRIPT
solution satisfying the boundary conditions, defined as XZ
Γ∈FhI
Γ
ϕ]dS + σ[w] · [ϕ
XZ
Γ∈FhB
Γ
σ(w − wB ) · ϕ dS ,
where the penalty weight σ is chosen as
σ|Γ =
CW . Γ ∈ Fh , diam(Γ) Re
(3.34)
CR IP T
Jh (w, ϕ ) =
with a suitable constant CW > 0, which guarantees the stability of the method. Its choice depends on the
AN US
variant of the used DGM (NIPG, IIPG or SIPG). Now we approximate the viscous term in the following form:
Vis ≈ ah (w, ϕ) + Jh (w, ϕ)
M
Finally we set
ϕ ∈ Shr . ∀ϕ
w, ϕ ∈ Shr ,
(3.35)
ED
ch (w, ϕ ) = bh (w, ϕ ) + ah (w, ϕ ) + Jh (w, ϕ ),
and arrive at the definiton of an approximate solution. Definition We say that the functions wh is the space semidiscrete solution of our problem, if the fol-
AC
CE
PT
lowing conditions are satisfied:
wh ∈ C 1 ([0, T ]; Shr ),
(3.36)
wh (0) = Πh w0 ,
(3.37)
d (wh , ϕh )Ω + ch (wh , ϕh ) = 0 dt
ϕh ∈ Shr , ∀t ∈ (0, T ), ∀ϕ
(3.38)
where Πh w0 is an Shr -approximation of w0 from the initial condition. In the time discretization we use the second-order backward difference formula (BDF) as for the in-
compressible flow. Definition We say that the finite sequence of functions wkh = wh (tk ),
11
k = 0, . . . , N is the approximate
ACCEPTED MANUSCRIPT
solution of (5)–(10) computed by the BDF-DGM, if the following conditions are satisfied:
0 We set w−1 h = wh .
(3.39) ϕh ∈ Shr , ∀ϕ
n = 1, 2, . . . .
(3.40)
CR IP T
w0h = Πh w0 , wn+1 ∈ Shr , n = 0, 1, . . . , h n+1 n n−1 3wh − 4wh + wh , ϕ h + ch (wn+1 h , ϕh) = 0 2∆t Ω
Another technique for the time discretization is the time discontinuous Galerkin method. Then we obtain the space-time discontinuous Galerkin method (STDGM). We again consider a partition 0 = t0 < t1 < . . . < tN = T of the time interval [0, T ] and denote Im = (tm−1 , tm ), I m = [tm−1 , tm ]. We define the
r,q S h,∆t
AN US
r,q r,q 4 space Sh,∆t = (S h,∆t ) , where
q X = φ ; φ| = ζ φ , where φ ∈ S , ζ ∈ P (I ) i i i hr i q m I m i=0
with integers p, q ≥ 1, Pq (Im ) denoting the space of all polynomials in t on Im of degree ≤ q and the space
M
r,q S hr defined in (3.27). For ϕ ∈ Sh,∆t we introduce the following notation:
ϕ±m = ϕ(tm± ) = lim ϕ(t),
ED
t→tm ±
{ϕ}m = ϕ+m − ϕ−m .
(3.41)
Definition The space-time discontinuous Galerkin approximate solution of (5)–(10) is defined as a
PT
function wh∆t satisfying
wh∆t (0+), ϕh
Ω
= w0 , ϕh
Ω
∀ϕh ∈ Shr
(3.42)
AC
CE
and the following conditions:
r,q 1) wh∆t ∈ Sh,∆t , ! ! Z ∂wh∆t (t), ϕh∆t + ah (wh∆t , ϕh∆t , t) dt 2) ∂t I Ω Zm + bh (wh∆t , ϕh∆t , t) + Jh (wh∆t , ϕh∆t , t) dt + ({wh∆t }m−1 , ϕh∆t (tm−1 +))Ω Im Z r,q = ℓh (whD , ϕh∆t , t) dt ∀ϕh∆t ∈ Sh,∆t , m = 1, . . . , N.
(3.43)
Im
The nonlinear discrete problems (3.39)–(3.40) and (3.42)–(3.43) are solved at each time level by a Oseen-like iterations. See, e.g., [13], [12].
12
ACCEPTED MANUSCRIPT
M(t) Uµ a a
L(t)
kaa
EA T
H
CR IP T
kHH
Figure 1: Elastically supported airfoil with two degrees of freedom.
4. Flow-induced airfoil vibrations
AN US
We consider an elastically supported airfoil ΓWt with two degrees of freedom: the vertical displacement H (positively oriented downwards) and the angle α of rotation around an elastic axis EA (positively oriented clockwise). See Figure 1.
The motion of the airfoil for large vibration amplitudes is described by the system of ordinary differen-
M
tial equations for unknowns H and α:
(4.44)
S α H¨ cos α + Iα α¨ + kαα α = M(t),
(4.45)
PT
ED
mH¨ + S α α¨ cos α − S α α˙ 2 sin α + kHH H = −L(t),
where the dot and two dots denote the first-order and the second-order time derivative, respectively. This
CE
system is derived from the Lagrange equations of the second art, see [28]. We use the following notation: L(t) - aerodynamic lift force, M(t) - aerodynamic torsional moment, m - airfoil mass, S α - static moment of the airfoil around the elastic axis EA, Iα - inertia moment of the airfoil around the elastic axis EA, kHH -
AC
bending stiffness, kαα - torsional stiffness. ˙ System (4.44) – (4.45) is equipped with the initial conditions prescribing H(0), α(0), H(0), α(0). ˙ The
aerodynamic lift force L and the torsional moment M are defined by
L=
−lρre f v2re f c
Z
2 X
τ2 j n j dS ,
M=
ΓWt j=1
lρre f v2re f c2
Z
ΓWt i, j=1
where
r1ort = −(x2 − xEA2 ),
r2ort = x1 − xEA1 , 13
2 X
τi j n j riort dS ,
(4.46)
ACCEPTED MANUSCRIPT
l is the airfoil depth and c is the length of the airfoil cord. By τi j we denote the components of the aerodynamic stress tensor, defined by (2.6) and (2.7) (taking into account that ∇ · v = 0 in the case of the incompressible flow), n = (n1 , n2 ) is the unit outer norm to ∂Ωt on ΓWt (pointing into the airfoil) and xEA = (xEA1 , xEA2 ) is the position of the elastic axis at time t.
CR IP T
System (4.44) – (4.45) is transformed to a first-order system and solved by the Runge-Kutta method. The construction of the ALE mapping, the domain velocity and the discretization in the framework of the ALE method has been carried out according to the papers [28] for incompressible case and [5] for compressible flow.
AN US
5. Numerical experiments 5.1. Flow in a simplified channel of human vocal tract
In numerical experiments we consider a channel modeling a simplified subglottal and supraglottal spaces constricted by two bumps of the length 15.4 mm in the glottis region, which represent a simplified model of vocal folds. The numerical results are obtained for a triangulation consisting of 7790 elements
M
(see Figure 2). We prescribe the following dimensional data: the inlet density ρ¯ in = 1.225 kg m−3 , the outlet atmospheric pressure p¯ out = 97611 Pa, the inlet velocity v¯ in = 4 m s−1 (Re = 4356, inlet Mach number
ED
M = 0.012) or v¯ in = 0.4 m s−1 (Re = 436, M = 0.0012), the fluid viscosity µ = 18 · 10−6 kg m−1 s−1 , the heat conduction coefficient k = 2.428 · 10−2 kg m s−2 K−1 , the specific heat at constant pressure
PT
c p = 1009.992 m2 s−2 K−1 and the Prandtl number Pr = 0.748. The Reynolds number was computed by (2.2), with Lre f = 0.016 m, ρre f = ρ¯ in and vre f = v¯ in .
CE
We use the time step ∆t = 0.00002 s and set T = 0.2 s. For computing the compressible viscous flow we use the NIPG version of the discontinuous Galerkin method with constant CW = 500 inside of the
AC
computational domain and CW = 5000 on the boundary. We compare the distribution of the pressure drop and velocity along the symmetry axis of the channel
for compressible and incompressible flow at time instants t = 0.02, 0.08 and 0.2 s. Results for v¯ in = 4 m s−1
are presented in Figures 3-6 and results for v¯ in = 0.4 m s−1 are shown in Figures 7-10. Here the notation x = x1 is used.
14
CR IP T
ACCEPTED MANUSCRIPT
Figure 2: Triangulation of the computational domain with boundary parts and sizes of the individual parts in mm.
p [Pa]
p [Pa]
1500
1500 1000 500 0 -500 -1000 -1500 -2000 -2500 -3000
500 0 -500 -1000 -1500 -3
-2
-1
0
1
2
3
4
5
x / Lref
AN US
1000
6
-3
-2
-1
0
1 2 x / Lref
3
4
5
6
ED
M
Figure 3: Distribution of the pressure drop for vin = 4 m s−1 along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
|u| [m/s] 70
60
PT
60
|u| [m/s]
50 40
CE
30
50 40 30 20
20
10
10 0
AC
-3
-2
-1
0
0 1
2
3
4
5
6
-3
x / Lref
-2
-1
0
1 2 x / Lref
3
4
5
6
Figure 4: Distribution of the magnitude of the velocity for vin = 4 m s−1 along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
15
ACCEPTED MANUSCRIPT
AN US
CR IP T
1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 61
AC
CE
PT
ED
M
Figure 5: Magnitude of the velocity for vin = 4 m s−1 in the case of incompressible flow at time instants t = 0.02, 0.08 and 0.2 s.
Figure 6: Magnitude of the velocity for vin = 4 m s−1 in the case of compressible flow at time instants t = 0.02, 0.08 and 0.2 s.
16
p [Pa] 10 8 6 4 2 0 -2 -4 -6 -8 -10
p [Pa] 15 10 5 0 -5 -10 -15 -20 -25 -3
-2
-1
0
1 2 x / Lref
3
4
5
6
-3
-2
CR IP T
ACCEPTED MANUSCRIPT
-1
0
1 2 x / Lref
3
4
5
6
AN US
Figure 7: Distribution of the pressure drop for vin = 0.4 m s−1 , along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
|u| [m/s]
|u| [m/s]
4.5 4 3.5 3 2.5 2 1.5 1 0.5 0
7 6 5 4
M
3
-2
-1
0
1 2 x / Lref
3
4
5
ED
-3
6
2 1 0 -3
-2
-1
0
1 2 x / Lref
3
4
5
6
AC
CE
PT
Figure 8: Distribution of the magnitude of the velocity for vin = 0.4 m s−1 , along the symmetry axis of the channel in dependence of the dimensionless length x/Lre f (see Figure 2), for incompressible flow (left) and compressible flow (right). The dashed blue diamond marked line stands for time instant t = 0.02 s, the green circle marked line for t = 0.08 s and finally the red triangle marked line for t = 0.2 s.
17
ACCEPTED MANUSCRIPT
AN US
CR IP T
0.1 0.4 0.7 1 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7 4 4.3
AC
CE
PT
ED
M
Figure 9: Magnitude of the velocity for vin = 0.4 m s−1 in the case of incompressible flow at time instants t = 0.02, 0.08 and 0.2 s.
Figure 10: Magnitude of the velocity for vin = 0.4 m s−1 in the case of compressible flow at time instants t = 0.02, 0.08 and 0.2 s.
From the above figures we see that the pressure drop between the inlet and outlet computed for incompressible flow is in a good correspondence with the results for compressible flow. The same is true for the velocity magnitude in the narrowest part of the channel. Both the velocity and pressure distributions are identical from the inlet up to the narrowest part of the channel. It is important that the vortex structures are similar for the incompressible and compressible flow. The differences may be caused by different jet declination behind the channel constriction and by the Coanda effect, which is moreover influenced by the 18
ACCEPTED MANUSCRIPT
interaction of the jet with the large-scale supraglottal eddies. 5.2. Flow around the NACA0012 profile We consider laminar flow around the NACA0012 profile for six different flow regimes characterized by
and the Reynolds number Re:
α = 1, Re = 160000,
(C2)
α = 3, Re = 160000,
(C3)
α = 5, Re = 160000,
(C4) (C5) (C6)
AN US
(C1)
CR IP T
the far-field velocity v¯ in = 8 m s−1 (the corresponding inlet Mach number is M = 0.024), angle of attack α
α = 1, Re = 5000,
α = 2, Re = 5000,
α = 3, Re = 5000.
ED
M
We evaluate the aerodynamic coefficients drag cD and lift cL . They are defined as
(cD , cL ) = 2
Z
τn dS ,
(5.47)
ΓW
PT
where ΓW is the surface of the profile, τ is the dimensionless stress tensor and n = (n1 , n2 ) is the unit outer normal to ∂Ω on ΓW .
AC
11).
CE
For each flow regime we carried out computations on the grid consisting of 17158 elements (see Figure
19
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 11: The computational grid around the NACA0012 profile.
M
We prescribe the following data: ρ¯ in = 1.225 kg m−3 , p¯ out = 97611 Pa, k = 2.428 · 10−2 kg m s−2 K−1 , c p = 1009.992 m2 s−2 K−1 and two values of the fluid viscosity µ = 1.8375 · 10−5 kg m−1 s−1 (corresponding
ED
to Re = 160000, Pr = 0.764) and µ = 5.88 · 10−4 kg m−1 s−1 (corresponding to Re = 5000, Pr = 24.459). Moreover ∆t, T and values of CW are chosen similarly as in Section 5.1. For computing the compressible viscous flow we again use the NIPG version of the discontinuous Galerkin method.
PT
For flow regimes (C1)–(C3) we observe small fluctuations of the drag cD and lift cL in dependence on time, as can be seen from Figures 12–17 both for incompressible and compressible cases. After a transient
CE
regime (0 < t < 0.5s) numerical results oscillate around mean values. Table 1 shows the values of the corresponding drag and lift coefficients and in case Re = 160000 also
AC
experimental data from [26]. For the case Re = 5000 we have not found any experimental data. For flow regimes (C1)–(C3) we write the values in the form ”mean value ± standard deviation”. For flow regimes
(C4)–(C6) we obtained a stationary flow and present the steady-state values of the drag and lift coefficients. In Table 1 it is possible to see an interesting fact that the experimental values of drag and lift are in the
intervals between the mean values numerically computed with the use of incompressible and compressible models. It is also interesting that in all cases ccompr. < cincompr. and ccompr. > cincompr. . D D L L
20
ACCEPTED MANUSCRIPT
ccompr. D
cincompr. D
exp. data
ccompr. L
cincompr. L
exp. data
(C1) (C2) (C3) (C4) (C5) (C6)
0.0096 ± [0.001] 0.0104 ± [0.002] 0.0136 ± [0.004] 0.0507 0.0517 0.0531
0.0116 ± [0.001] 0.0132 ± [0.002] 0.0169 ± [0.004] 0.0533 0.0541 0.0566
0.0104 0.0114 0.0140 -
0.155 ± [0.041] 0.382 ± [0.027] 0.562 ± [0.036] 0.0374 0.0727 0.0988
0.086 ± [0.046] 0.30 ± [0.043] 0.498 ± [0.047] 0.0366 0.0696 0.1062
0.11 0.33 0.55 -
0.02
0.02
0.018
0.018
0.016
0.016
0.014
0.014 CD
0.012 0.01
CR IP T
Flow
AN US
CD
Table 1: Values of the computed aerodynamic coefficients cD and cL for incompressible and compressible viscous flow and experimental data.
0.012 0.01
0.008
0.008
0.006
0.006
0.5
1
1.5
2
2.5
(a) Incompressible flow
0.5
1
1.5
2
2.5
t[s]
(b) Compressible flow
M
t[s]
ED
Figure 12: Values of the computed aerodynamic coefficients cD for the incompressible flow (left) and compressible flow (right) in flow regime (C1) in dependence on time t[s].
0.4
0.2
CE
0.15
0.3 0.25 CL
0.3 0.25 CL
0.4 0.35
PT
0.35
0.2 0.15
0.1
0.1
0.05
0.05
0
AC
0.5
1
0 1.5
2
2.5
0.5
t[s]
1
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
Figure 13: Values of the computed aerodynamic coefficients cL for the incompressible flow (left) and compressible flow (right) in flow regime (C1) in dependence on time t[s].
21
0.025
0.025
0.02
0.02
0.015
0.01
0.015
0.01
0.005
0.005 0.5
1
1.5
2
2.5
0.5
t[s]
CR IP T
CD
CD
ACCEPTED MANUSCRIPT
1
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
AN US
Figure 14: Values of the computed aerodynamic coefficients cD for the incompressible flow (left) and compressible flow (right) in flow regime (C2) in dependence on time t[s].
0.5
0.5
0.45
0.45
0.4
0.4
CL
0.35
0.3
M
CL
0.35
0.25 0.2 0.15
0.3
0.25 0.2 0.15
0.1
0.1 1
1.5
2
ED
0.5
2.5
0.5
1
t[s]
1.5
2
2.5
t[s]
(b) Compressible flow
PT
(a) Incompressible flow
CE
Figure 15: Values of the computed aerodynamic coefficients cL for the incompressible flow (left) and compressible flow (right) in flow regime (C2) in dependence on time t[s].
0.04
AC
0.04
0.035
0.03
0.03
0.025
0.025 CD
CD
0.035
0.02
0.015
0.02 0.015
0.01
0.01
0.005
0.005
0
0 0.5
1
1.5
2
2.5
0.5
t[s]
1
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
Figure 16: Values of the computed aerodynamic coefficients cD for the incompressible flow (left) and compressible flow (right) in flow regime (C3) in dependence on time t[s].
22
ACCEPTED MANUSCRIPT
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3 0.5
1
1.5
2
2.5
0.5
t[s]
1
CR IP T
CL
0.8
CL
0.8
1.5
2
2.5
t[s]
(a) Incompressible flow
(b) Compressible flow
AN US
Figure 17: Values of the computed aerodynamic coefficients cL for the incompressible flow (left) and compressible flow (right) in flow regime (C3) in dependence on time t[s].
5.3. Numerical experiments of the flow-induced airfoil vibrations
Following the example of numerical simulations of flow-induced vibrations published in [28], we prescribe the input data as follows: m = 0.086622 kg, S α = −0.000779673 kg m, Iα = 0.000487291 kg m−2 ,
M
kHH = 105.109 N m−1 , kαα = 3.695582 N m rad−1 , the fluid viscosity µ = 1.8375 · 10−5 kg m−1 s−1 , the density ρ¯ in = 1.225 kg m−3 , the pressure p¯ out = 101250 Pa, k = 2.428 · 10−2 kg m s−2 K−1 , c p =
ED
1009.992 m2 s−2 K−1 , Pr = 0.764. The elastic axis EA is located at 40% of the profile cord.The airfoil cord c = 0.3 m, airfoil depth l = 0.05 m, and the following initial conditions for the structural equations
PT
˙ H(0) = −20 mm, α(0) = 6, H(0) = α(0) ˙ = 0. We prescribed inlet velocities v¯ in = 10 m s−1 (corresponding to Re = 200000 and Mach number M = 0.0294) and v¯ in = 20 m s−1 (corresponding to Re = 400000 and
CE
Mach number M = 0.0588). The computation started at a time instant t = −δ < 0 with a fixed airfoil. Then at time t = 0 the airfoil was released and the fluid-structure interaction process followed.
AC
Incompressible flow was computed by the method described in section 3.1 with the use of its ALE version. For compressible flow we used the SIPG STDG method with quadratic polynomials in space (r = 2) and linear polynomials in time (q = 1) (the results are denoted by ST-DGr2q1). The parameter CW was chosen similarly as in section 5.1. Figures 18–21 show the comparison of the displacement H and the angle α computed for incompressible and compressible flow. The time responses of the airfoil are damped due to the unsteady aerodynamic forces given by (4.46) because the considered inlet flow velocities are far below the critical flow velocity 42 m s−1 for flutter, see [28]. The responses computed by DGFEM for the compressible flow are in good agreement with the results computed by FEM for the incompressible flow; the maximal displacements for t > 0 are very similar for 23
ACCEPTED MANUSCRIPT
both the translation H and the rotation α, and the differences in frequencies of oscillations are less than 1%
8
15
6
10
4
5
CR IP T
20
2
0
α
H [mm]
for the inlet flow velocity vin = 10 m s−1 and less then 4% for vin = 20 m s−1 .
0
-5 -10
-2
-15
-4
-20 0
0.5
1
1.5
-6
2
0
0.5
t[s]
1
1.5
2
AN US
t[s]
Figure 18: Values of the computed displacement H and angle α for the compressible (black) and the incompressible flow (green) for v¯ in = 10 m s−1 in dependence on time t[s].
8
15
6
M
0 -5 -10 -15 -20 0
0.2
0.4
ED
H [mm]
5
0.6
0.8
4 2
α
10
0 -2 -4 -6
1
0
t[s]
0.2
0.4
0.6
0.8
1
PT
t[s]
CE
Figure 19: Values of the computed displacement H and angle α for the compressible (black) and the incompressible flow (green) for v¯ in = 20 m s−1 in dependence on time t[s].
AC
6. Conclusion
The paper is devoted to the analysis of the discontinuous Galerkin solution of compressible low Mach
number flows. The goal is to prove the robustness of the developed DGM with respect to the Mach number at incompressible limit. To this end, several test problems were solved by the DGM and results were compared with their numerical simulation with the aid of the stabilized finite element method applied to incompressible Navier-Stokes equations. The performed numerical tests show that the correspondence of compressible and incompressible results is good not only for the stationary cases but also for the case of time dependent flow. The time dependent character of the solution is caused either by moving vortices or by interaction of fluid with an elastically supported vibrating airfoil. In all these cases the DGM provides results comparable to the solution of the incompressible model. 24
ACCEPTED MANUSCRIPT
As a conclusion, it is possible to say that the described numerical method for the solution of the compressible Navier-Stokes equations allows the simulation of very low Mach number flows and yields the results comparable with results obtained with the aid of the incompressible model, both in internal aerodynamics, for the airflow in a constricted channel (M = 0.0012), and in outer aerodynamics for the fixed
CR IP T
airfoil (M = 0.024) as well as for the flow induced airfoil vibrations (M = 0.0294). The results prove the robustness of the discontinuous Galerkin method with respect to the low Mach number. Acknowledgment
The work of M. Bal´azsov´a was supported by the Grant No. SVV-2017-260455 of the Charles University. The
AN US
research of M. Feistauer was supported by the Grant No. 17-01747S of the Czech Science Foundation and finally the research of P. Sv´acˇ ek and J. Hor´acˇ ek was supported by the Grant No. 16-01246S of the Czech Science Foundation.
References
[1] Bassi, F., Rebay, S.: High-order accurate discontinuous finite element solution of the 2D Euler equations. J. Comput. Phys. 138 (1997), 251–285.
Methods Fluids 31 (1999), 79–95.
M
[2] Baumann, C. E., Oden, J. T.: A discontinuous hp finite element method for the Euler and Navier-Stokes equations. Int. J. Numer.
ED
[3] Bruneau, Ch.-H., Fabrie, P.: Effective downstream boundary conditions for incompressible Navier–Stokes equations. Int. J. Nu-
mer. Methods Fluids 19(8) (1994), 693–705.
ˇ [4] Cesenek, J., Feistauer, M., Hor´acˇ ek, J., Kuˇcera, V., Prokopov´a, J.: Simulation of compressible viscous flow in time-dependent
PT
domains. Appl. Math. Comput. 13 (2013), 7139–7150. ˇ [5] Cesenek, J., Feistauer, M., Kos´ık, A.: DGFEM for the analysis of airfoil vibrations induced by compressible flow. ZAMM 93(13) (2013), 387–402.
CE
[6] Cockburn, B., Karniadakis, G.E., Shu C.W., Eds.: Discontinuous Galerkin Methods, Lecture Notes in Computational Science and Engineering 11, Springer, Berlin (2000). [7] Darwish, M., Moukalled, F., Sekar, B.: A robust multi-grid pressure-based algorithm for multi-fluid flow at all speeds, Int. J.
AC
Numer. Methods Fluids 41 (2003) 1221–1251.
[8] Davis, T. A., Duff, I. S.: A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Trans. Math. Softw. 25 (1999) 1–19.
[9] Dolejˇs´ı, V.: Semi-implicit interior penalty discontinuous Galerkin methods for viscous compressible flows. Commun. Comput.
Phys. 4 (2008), 231–274. [10] Dolejˇs´ı, V., Feistauer, M.: Discontinuous Galerkin Method - Analysis and Applications to Compressible Flow. Springer, (2015). [11] Dumbser, M.: Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier-Stokes equations.
Comput. Fluids 39(1) (2010), 60–76. [12] Feistauer, M., Dolejˇs´ı, V., Kuˇcera, V.: On the discontinuous Galerkin method for the simulation of compressible flow with wide range of Mach numbers. Comput. Vis. Sci. 10 (2007), 17–27. [13] Feistauer, M., Kuˇcera, V.: On a robust discontinuous Galerkin technique for the solution of compressible Flow. J. Comput. Phys. 224 (2007), 208–221.
25
ACCEPTED MANUSCRIPT
[14] Feistauer, M., Hasnedlov´a-Prokopov´a, J., Hor´acˇ ek, J., Kos´ık, A., Kuˇcera, V.: DGFEM for dynamical systems describing interaction of compressible fluid and structures. J. Comput. Appl. Math. 254 (2013), 17–30. [15] Gelhard, T., Lube, G., Olshanskii, M. A., Starcke, J.-H.: Stabilized finite element schemes with LBB-stable elements for incompressible flows. J. Comput. Appl. Math. 177 (2005), 243–267. [16] Hartmann, R., Houston, P.: Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations, J.
CR IP T
Comput. Phys. 183 (2002), 508–532.
[17] Klein, R.: Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics 1: one-dimensional flow,
J. Comput. Phys. 121 (1995) 213–237.
[18] Meister, A.: Comparison of different Krylov subspace methods embedded in an implicit finite volume scheme for the computation of viscous and inviscid flow fields on unstructured grids, J. Comput. Phys. 140 (1998) 311–345.
[19] Meister, A.: Viscous flow fields at all speeds: Analysis and numerical simulation, J. Appl. Math. and Physics 54 (2003) 1010–
AN US
1049.
[20] Meister, A., Struckmeier, J.: Hyperbolic Partial Differential Equations, Theory, Numerics and Applications, Vieweg, 2002. [21] Nerinckx, K., Vierendeels, J., Dick, E.: Mach-uniformity through the coupled pressure and temperature correction algorithm. J.
Comput. Physics 206 (2005), 597-623.
[22] Nomura, T., Hughes, T.J.R.: An arbitrary Lagrangian-Eulerian finite element method for interaction of fluid and a rigid body,
Comput. Methods Appl. Mech. Engrg. 95 (1992), 115–138.
[23] Park, J.H., Munz, C.-D.: Multiple pressure variables methods for fluid flow at all Mach numbers, Int. J. Numer. Methods Fluids
M
49 (2005) 905–931.
[24] Punˇcoch´aˇrov´a-Poˇr´ızkov´a, P., F¨urst, J., Hor´acˇ ek, J., Kozel, K.: Numerical solutions of unsteady flows with low inlet Mach
ED
numbers. Math. Comput. Simul. 80(8) (2010), 1795–1805.
[25] Roller, S., Munz, C.-D., Geratz, K. J., Klein, R.: The multiple pressure variables method for weakly compressible fluids, Z. Angew. Math. Mech. 77 (1997) 481–484.
PT
[26] Sheldahl, R. E., Klimas, P. C.: Aerodynamic Characteristics of Seven Airfoil Sections Through 180 Degrees Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbines, SAND80-2114, March 1981, Sandia National Laboratories, Albuquerque, New Mexico.
CE
[27] Seo, J. H., Mittal, R.: A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries, J. Comput. Phys. 77 (1997) 481–484. [28] Sv´acˇ ek, P., Feistauer, M., Hor´acˇ ek, J.: Numerical simulation of flow induced airfoil vibrations with large amplitudes. J. of Fluids
AC
Struct. 23 (2007), 391–411. [29] van der Vegt, J. J. W., van der Ven, H.: Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flow. J. Comput. Phys. 182 (2002), 546–585.
[30] van der Heul, D.R., Vuik, C., Wesseling, P.: A conservative pressure-correction method for flow at all speeds, Comput. Fluids 32 (2003) 1113–1132. [31] Wenneker, I., Segal, A., Wesseling, P.: A Mach-uniform unstructured steggered grid method, Int. J. Numer. Methods Fluids 40 (2002) 1209–1235. [32] Wesseling, P.: Principles of Computational Fluid Dynamics, Springer, Berlin, 2001. [33] Wesseling, P., van der Heul, D.R.: Uniformly effective numerical methods for hyperbolic systems, Computing 66 (2001) 249– 267.
26