On numerical simulation of ﬂow problems in three dimension: energy conservation in ﬂuid-structure interactions

. This paper is interested in three-dimensional ﬂow problems, where the attention is paid to two topics: ﬁrst, the numerical approximation of turbulent ﬂows, second, the relation between the incompressible ﬂow and the motion of immersed solid body. The considered system shall be studied with respect to the approximation point of view, and the variational formulation shall be discussed.


Introduction
Mathematical modelling is important in many engineering problems of fluid and solid mechanics, the use of proper mathematical model is usually crucial to get correct results. Recently, also the numerical simulations are being used for solution of large number of technical or scientific problems and, due to the increase of the computer power, the simulations of complex mathematical models are possible, cf. [1], [2], [3], [4]. Particularly, the turbulent flow under different flow conditions is modelled and numerically approximated in various applications, see e.g. [5], [6], [7] and similarly fluid-structure interactions are investigated [8], [9], [10], [11], [12], [13]. On the other hand, the effect of numerical errors and/or mathematical model is usually not analyzed.
Here, the attention is paid to the correctness of the variational formulation of the three dimensional fluid flow interacting with a solid motion, where a basic energy estimate is derived for a simplified problem. The outflow boundary condition based on [14] are used in order to treat the possible difficulties caused by the outflow boundary (see also [15], [16])). Next, the numerical approximation of the flow in a three dimensional channel is considered, see also [17]. The governing equations are then approximated by the in-house implementation of finite-element and finite volume methods, the description of the finite element method is given and numerical results are shown.

Mathematical Description
Mathematical model. For the flow model, we consider the time dependent computational domain Ω t ⊂ R 3 with the Lipschitz continuous boundary ∂Ω t . The fluid motion in the domain Ω t is described by the Navier-Stokes system a e-mail: Petr.Svacek@fs.cvut.cz are the components of the symmetric part of the gradient of u denoted by S = S(u) , p is the kinematic pressure (i.e., the pressure divided by the constant fluid density ρ), ν is the kinematic viscosity of the fluid (i.e. the viscosity divided by the density ρ).
Let us note, that in the case of the Reynolds averaged Navier-Stokes equations the viscosity coefficient is replaced by ν eff = ν + ν T , where ν T is a turbulent viscosity (obtained by an additional model). In this case, only the mean part of the velocity vector and mean part of the pressure are modelled, cf. [18]. In what follow, we shall focus only on Navier-Stokes equations.
The system (1) is equipped with boundary conditions prescribed on the mutually disjoint parts of the boundary where p re f denotes a reference pressure, α − = min(0, α) denotes the negative part of the number α ∈ R and w D is the velocity of the boundary Γ Wt . For the sake of simplicity, let us set p re f = 0. Further, the system (1) is equipped with an initial condition u(x, 0) = u 0 (x), x ∈ Ω 0 . The boundary condition (2c) is the modification of the well known donothing boundary condition, see [14].
ALE method. In order to practically treat the motion of the domain Ω t , the Arbitrary Lagrangian-Eulerian (ALE) method is used, cf. [19]. The ALE mapping A = A(ξ, t) = A t (ξ) defined for all t ∈ (0, T ) and ξ ∈ Ω 0 is assumed to be sufficiently smooth with the continuous and bounded Jacobian and the domain velocity The Navier-Stokes equations (1) can be equivalently written either in the non-conservative formulation or in the ALE conservative way Structure motion. In this paper, the solid body which can be displaced by h = h(t) ∈ R 3 . For the sake of brevity, only the rotation around the elastic axis (EA) placed at the center of gravity (CoG) parallel to z-axis is allowed by the angle α = α(t). The equation of motion reads where m is the mass of the airfoil, I α is the moment of inertia around the (EA). The diagonal matrix K denote the stiffness coefficients. On the right-hand side the aerodynamical force F(t) and the aerodynamical moment M 3 (t) are involved given by where ε mik is the Levi-Civita symbol, r = (r 1 , r 2 , r 3 ) and and x CG is the position of the CoG at the time instant t and n is the unit outward normal to ∂Ω t . Let us mention, that this is only a model problem, used to derive the variational weak form of the problem and show, that a basic energy estimate can be well estabilished in this case. However, almost the same approach can be used for more relevant technical applications as motion of a solid particle driven by a fluid flow.

Weak formulation
For the weak formulation the function spaces Q t = L 2 (Ω t ) and W t = H 1 (Ω t ) are employed for the pressure and the velocity, respectively, defined at any t ∈ (0, T ). The space for test functions is defined by Flow problem. The weak formulation of the Navier-Stokes equations reads: Find U = (u, p) ∈ W t × Q t such that u approximately satisfy boundary conditions (2) and and where by (·, ·) G the dot-product in L 2 (G) is denoted. Similarly, the weak formulation of ALE non-conservative equations (3) read: Find U = (u, p) ∈ W t × Q t such that u approximately satisfy boundary conditions (2a,b) and It should be noted, that both the formulations (8) and (9) are formally equivalent.
Aerodynamical forces. In order to formulate the aerodynamical forces F weakly, we shall use a function ϕ ∈ H 1 (Ω t ) such that ϕ(x, t) = 1 for x ∈ Γ Wt , and its compact support supp ϕ ⊂ Ω t ∪ Γ Wt . Multiplying the first equation applying Green's theorem to viscous and pressure terms, using the notation w = u − w D , we get Thus with the aid of (6), and having Ψ k,h i = δ ik on Γ Wt , we get the weak form of the components of the aerodynamical force F : where (σ i, j n j ) · Ψ α i = ε 3ik σ i, j n j r k on Γ Wt , and thus Structure displacement and grid velocity. A point ξ ∈ Γ W 0 from the airfoil surface is at time instant t transformed to the point x(t) ∈ Γ Wt given by where ξ CoG denotes the location of the CoG in the reference domain (i.e. at the time t = 0). Now, with the aid of the previously defined functions Ψ i,h and Ψ α , the domain velocity w D on the surface of the airfoil Γ Wt satisfies the relation for any x ∈ Γ Wt and t ∈ [0, T ]. Without loss of generality we can assume that the function Ψ h,i (x, t), Ψ α (x, t) are chosen such that equation (13) holds for any x ∈ Ω t . Further, we multiply equation (10) byḣ , multiply equation (11) bẏ α, sum them, and geṫ Substituting for L and M from equations (5) to the lefthand side of equation (14) leads tȯ hL +αM = d dt Further, the right hand side of the equation (14) has a similar form to the (non-conservative) ALE weak formulation (9), if we take the test function z = w D =ḣΨ h (x, t) + αΨ α (x, t). In order to get an energy estimate for solution of (9), we continue only with a simplified problem.
Energy estimate for the simplified model. Assuming now, that u D ≡ 0 and Γ O = ∅, we write the solution u as u = (u − w D ) + w D , where the first part z = (u − w D ) ∈ X t , and where for the last term we already found the identity (14). Now, if u D ≡ 0, then an apriori estimate for the simplified coupled problem can be derived in the form which is the energy estimate consisting of the sum of the fluid kinetic energy, the dissipation term, the structural potential energy and the structural kinetic energy, respectively.

Numerical approximation
Time discretization For the sake of simplicity, we consider the equidistant partition t k = kΔt of the time interval I with a time step Δt > 0, and denote the approximations u k ≈ u(·, t k ) and p k ≈ p(·, t k ). Moreover we approximate the domain velocity w D at time level t k by w k D . Furthermore, in what follows only the non-conservative ALE formulation shall be used, although in ALE frame the conservative ALE formulation can have an influence on the quality of the numerical solution. Let us focus on the description of the discretization at an arbitrary time step t = t n+1 , which is kept fixed in this section (for the sake of simplicity we shall omit the subscripts t and t n+1 ). We shall consider all the function spaces X, W , Q defined for this time instant t = t n+1 on the domain Ω := Ω t n+1 . Then the time derivative in the weak formulation (8) is approximated at the time t = t n+1 by the second order backward difference formula, i.e.
where U = (u, p), V = (z, q), U * = (u * , p). Weak ALE formulation of the time discretized problem reads: Find U = (u n+1 , p n+1 ) ∈ W × Q such that for all test functions V = (z, q) where q ∈ Q and z ∈ X holds a(U, U, V) = L(V). Finite element spaces Further, the solution u and p is sought on a couple of the finite element spaces W Δ ⊂ H 1 (Ω n+1 ) and Q Δ ⊂ L 2 (Ω n+1 ) for the approximation of the velocity components and the pressure. The finite element spaces W Δ , X Δ , Q Δ are defined over an admissible triangulation T Δ of the domain Ω, cf. [20], formed by a finite number of closed elements K ∈ T Δ with the following properties: The intersection of two different elements K, K ∈ T Δ is either empty or a common edge or a common vertex or a common (triangular or quadrilateral) face of these elements, A3 K ∈ T Δ is either tetrahedron, hexahedron, pyramid or prism, see figures 2-3.
Further, on the reference tetrahedral and hexahedral ele-mentsK we define P K as piecewise linear functions P K = P 1 (K) and piecewise trilinear functions P K = Q 1 (K), respectively. For the prism reference element the space P K is chosen as linear functions in the reference variables ξ 1 , ξ 2 and bi-linear in ξ 1 , ξ 3 and ξ 2 , ξ 3 . The definition of P K on the pyramidal element is slightly more complicated, see [21], but the definition allows the linearity of the shape functions on each face, which makes the finite element space consistent.
The finite element spaces are defined by By X Δ ⊂ W Δ the subspace of the test functions is denoted. The introduced finite element spaces do not satisfy the Babuška-Brezzi (BB) inf-sup condition [see, e.g., [22], [23] or [24]], and the problem requires to use a stabilization, here based on the pressure stabilizing/Petrov-Galerkin method (PSPG), see, e.g., [25].
Stabilized problem In order to stabilize the method the fully stabilized scheme is used, which consists of streamlineupwind/Petrov-Galerkin(SUPG) and PSPG stabilization combined with the div-div stabilization, cf. [25]. The stabilized discrete problem reads: Find U = (u n+1 Δ , p n+1 Δ ) ∈ W Δ × Q Δ such that u n+1 satisfies approximately the Dirichlet boundary conditions (2,a,b) and a(U; U, V) + L(U; U, V) + P(U, V) = L(V) + F (V), (18) holds for all V = (z, q) ∈ X Δ × Q Δ , where the terms L and F are the SUPG/PSPG terms defined by where the function w n+1 = u * − w n+1 D stands for the transport velocity, and the terms R a K and R f K are parts of the local residual defined by The div-div stabilizing terms P(U, V) read Here, the choice of the parameters δ K and τ K is carried out according to [25] or [26] on the basis of the local element length h K , i.e.
where the local Reynolds number Re loc is defined as EPJ Web of Conferences 02089-p.4

Numerical results and conclusion
Numerical results The finite element method was implemented by an in-house software, and applied for solution of flow in bent channels shown in figure 4. The comparison of the results obtained by the presented finiteelement method to the finite volume method is shown in figures 5 and 6 for the flow in bent channels with 90 or 180 degrees. The Reynolds number for the considered case was chosen as Re = 600. For the bent channel with 90 degrees, the separation zone appears close to the bent part of the boundary. Both the methods predicted similar structure of the separation zones. For the finite volume code, structured multiblock hexahedral mesh was used with approximately 60000 cells, for the finite element the unstructured hexahedral mesh was used with either 60000 or 120000 elements (without significant differences). For both methods, the mesh was refined nearby the boundary to capture the boundary layer. The results shows very good agreement of both methods, see figure 5. Similarly, the flow in the bent channel with 180 degrees was approximated, where slightly higher number of elements were used (approximately 70000 both for finite volume and finite elment cases). Both methods leads to very similar results, see figure 6.

Conclusion
The presented paper is devoted to two topics, the first one is the variational formulation of a problem of a solid body immersed in flowing fluid, where the problem was thouroughly analyzed and an a priori energy estimate was shown under simplifying assumptions. Let us mention, that this estimate can be very important for numerical approximation espesically in the case of fluidstructure interactions [27] or motion of particles immersed in fluid driven by its flow. Let us mention, that the estimates shown here are the base for the estimates required for the discrete solution, for the discrete problem in 2D the problem was analyzed in [27], [28]. This situation is however much more complicated for 3D. Furthermore, the finite element method was used for approximation of 3D flow in two channels and the results were verified by comparison with numerical results of finite volume method. The solution by finite element method is based on the solution of the coupled problem both for velocity and pressure unknowns. On the other hand, for the finite volume method is based on the artificial compressibility method and the viscous terms are reconstructed from dual volumes. As the underlying numerical methods differs significantly, the presented results can serve as a verification of the results.