Comparison of advanced turbulence modeling approaches for fluid-structure interaction

Comparison of advanced turbulence modeling approaches for fluid-structure interaction VI International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2015 B. Schrefler, E. O˜nate and M. Papadrakakis(Eds) COMPARISON OF ADVANCED TURBULENCE MODELING APPROACHES FOR FLUID-STRUCTURE INTERACTION A. Ali∗†, T. Reimann†, D.C. Sternel†and M. Sch¨afer† †Institute of Numerical Methods in Mechanical Engineering Technische Universit¨at Darmstadt Dolivostraße 15, 64293 Darmstadt, Germany e-mail: [email protected], web page: http://fnb.tu-darmstadt.de Key words: Fluid-Structure Interaction, Turbulence, Delayed Detached-Eddy Simulation, Large Eddy Simulation Abstract. In this study we present the results of the benchmark of a turbulent FluidStructure Interaction test case. An implicit partitioned approach is employed to couple the fluid and structure subproblems. We employ three different techniques to model the turbulence in fluid motion. A 2-d unsteady Reynolds Averaged Navier-Stokes approach, with an elliptic relaxation based turbulence model (ζ−f), successfully captures the oscillation mode. Further investigations are performed with a Delayed Detached Eddy Simulation and a Large Eddy Simulation model. The ζ−fmodel is used as a baseline unsteady Reynolds Averaged Navier-Stokes model for the Delayed Detached Eddy Simulation. A comparison of the structural deflections from the simulations show a reasonable agreement with the experiment. In light of the presented results, the suitability of the modeling approaches is discussed. 1 INTRODUCTION Fluid-Structure Interaction (FSI) phenomena are important to study the design of many engineering applications. Experiments for most of the real-world FSI applications are not feasible – measurement techniques can not look inside all the important parameters, and are too expensive. With increasing computational power, simulation of such multi physics scenarios are becoming feasible and can give new knowledge. The capabilities of these numerical methods to study FSI needs to be validated against benchmark test cases. A great number of FSI applications have turbulent fluid motion, thus making it important to study the FSI phenomena in turbulent flows. The numerical and experimental studies on FSI, conducted in the last decade mostly focused on laminar flows. In an effort 1 512 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer to provide experimental data for validation of numerical tools, Gomes and Lienhart [6] proposed the first validation test case for FSI with incompressible turbulent fluid motion. The structural model of this test case, contains a thin flexible sheet attached to a revolvable cylinder with a rectangular mass attached to the other end of the sheet. The structure is placed inside a vertical tunnel with the flow Reynolds Number (Re) of 15,000, based on the cylinder diameter. The structure exhibits a periodic oscillation close to its natural frequency. This instability mechanism is characterized as Instability Induced Excitation (IIE) [14] with the structure swiveling in the first mode. Recently De Nayer and Kalmbach [4], as well as Kalmbach and Breuer [11] have proposed two different test cases with turbulent fluid motion. The proposed benchmark test cases offer a simpler structure geometry and provide a different excitation mechanism and motion mode, which is not present in [6]. The Rebased on the cylinder diameter is in sub-critical regime [15] (103<R e<2×105). This flow configuration is considered challenging for the turbulence models, since the boundary layer is laminar and transition to turbulence occurs in the separated shear layers and the wake. The present investigation aims to access the capabilities of different turbulence modeling techniques for coupled FSI problems. The turbulent test case presented in [6] is simulated employing three different turbulence modeling techniques. The ζ−fmodel proposed in [8] is utilized to perform the study in 2-d Unsteady Reynolds Averaged NavierStokes (URANS) flow simulation. The Smagorinsky model [18] with dynamic procedure suggested by Germano [5] is used to carry out a Large Eddy Simulation (LES). A hybrid URANS/LES based on the Detached Eddy Simulation model (DDES) [19] with ζ−fas baseline URANS model, is also tested. The DDES formulation of ζ−fmodel has been discussed and verified in [24]. The structural deflections, swiveling frequency and the end mass phase delay from simulations are compared with the experimental data. 2 GOVERNING EQUATIONS For the fluid subdomain Ωf, the fluid is assumed to be Newtonian with incompressible fluid motion. The basic conservation equations governing transport of mass and momentum are given as ∂vi ∂xi =0,(1) ρf Dvi Dt =−∂p ∂xi +µf ∂2vi ∂x2 j +ρffi,(2) where viis the velocity vector, pis the static pressure, µfis the dynamic fluid viscosity, ρfis the fluid density and firepresent the external force vector. For the structure subdomain Ωs, we define a material point Xin the reference configuration. The function χrepresents the transformation from Xto xas xi=χ(Xj,t),(3) 2 513 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer where Xjare the components of position vector Xand xirepresents the components of current spatial position vector x. The displacements are then defined as ui=xi−Xi.(4) The basic equation of momentum balance for solid domain Ωsis written as ρs ∂2χ(Xj,t) ∂t2=∂SjiFij ∂Xj +ρsfi,(5) where Sji is the second Piola-Kirchhoff stress tensor, ρsis the density of solid material and fsrepresents the external forces on the solid. Fij =∂xi/∂Xjis the deformation gradient. In the present study, the material is modeled utilizing a simple hyper-elastic material model, the Saint Venant-Kirchhoff law (for details see [16, 23]). For the second Piola-Kirchhoff stress tensor the model states Sij =λsEkkδij +2µsEij,(6) where the Green-Lagrange strain tensor is represented as Eij =1 2(Fkifkj −δij),(7) with λsand µsas Lam´e constants. The two subproblems are coupled at the boundary with suitable interface and boundary conditions. The standard boundary conditions apply on the fluid boundaries Γfand the structure boundaries Γs. The following conditions on velocities and stresses are applied at the fluid-structure interface vf iΓf∩Γs =˙ub iand σs ijΓf∩Γs=σf ijΓf∩Γs ,(8) where ˙ub iis the velocity of the interface and σs ij and σf ij represent the Cauchy stress tensor of the solid and the fluid domain, respectively. 3 MODELING APPROACH This section gives a brief description of turbulence modeling approaches in this study. 3.1 ζ−fModel The ζ−fmodel proposed in [8] is a linear eddy-viscosity model. The model is capable of predicting the anisotropic behavior of turbulence near walls by evaluating the eddyviscosity νtbased on the wall normal velocity scale ratio ζas νt=Cµζkτ, where ζ=v2 2/k and kis the kinetic energy of turbulence. The constitutive model equations are given as 3 514 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer Dk Dt =P−+∂ ∂xiν+νt σk∂k ∂xi,(9) D Dt =C1P−C2 τ+∂ ∂xiν+νt σ∂ ∂xi,(10) Dζ Dt =f−ζ kP+∂ ∂xiν+νt σζ∂ζ ∂xi,(11) L2∇f−f=1 τc1+C2 P ζ−2 3,(12) where is the dissipation of turbulent kinetic energy, Pis the production of turbulent kinetic energy and fis the elliptic relaxation term that models the pressure-velocity correlations. Land τare the length and time scales of turbulence, whereas other unknown terms in the given set of equations are model constants. For a detailed model description see [8]. 3.2 Dynamic Smagorinsky Model The LES simulation in this study is performed by applying the Smagorinsky model [18] to estimate the Sub-Grid Scale (SGS) turbulent viscosity νSGS as νSGS =Cs∆2|S|,(13) where |S|= (2SijSij)1/2is the magnitude of strain-rate tensor Sij =∂vi/∂xj+∂vj/∂xi, Csis the model constant and ∆ = (∆1∆2∆3)1/3is the filter width, with ∆irepresenting filter width in each spatial direction. An overbar on Sij represents a filtered quantity. The model constant Cs=Cs(x, t) is calculated dynamically as proposed by Germano et al. [5]. The resulting equation system to estimate Csis solved using least squares method as suggested by Lilly [12]. The Csvalues are clipped as Cs(x, t) = max{Cs(x, t),0}to avoid negative values of Cs. 3.3 ζ−fDDES The Detached Eddy Simulation (DES) concept first proposed by Spalart [20], is to combine URANS and LES to have a model with better prediction of turbulence than a URANS and computationally less expensive than an LES. In the DES approach, the URANS model is modified to achieve a SGS model in regions where grid is fine enough for an LES. The switching between two modes is based on the length scale lturb as lturb = min(lRANS,C DES∆),(14) where ∆ = max(∆1,∆2,∆3) and CDES is a model constant. lturb is introduced in the URANS model by modifying the dissipation term in the transport equation (9) for the 4 515 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer turbulent kinetic energy as =k3/2/lturb.lturb either becomes the original URANS length scale (lRANS <C DES∆) or the SGS length scale (lRANS >C DES∆). The Delayed Detached Eddy Simulation (DDES) model was proposed in [19] as an improvement for some deficiencies of the original DES model. The modification and verification of ζ−fmodel to perform a DDES are presented in [24]. For DDES, the DES length scale is modified to incorporate a shielding function to preserve the URANS mode in boundary layers, because of the grid clustering near the boundaries CDES∆<l RANS. A quantity rdis defined as rd=νt+ν √vi,jvi,jκ2d2,(15) where νis the molecular viscosity, vi,j are the velocity gradients, κis the K´arm´an constant and dis the wall distance. The shielding function fdis defined as a function of rdas fd=1−tanh([8rd]3).(16) The function fdis designed to be 0, to prevent activation of LES mode in boundary layer regions. The new length scale for DDES is then defined as lDDES =lRANS −fdmax(0,d−CDES∆ψ),(17) where ψis the term added to eliminate the influence of low-Returbulence models in the SGS mode, which is given as ψ=C1 C2Cµζ3/2 .(18) 4 EXPERIMENTAL SETUP The structural model for this test case consists of a flexible stainless steel sheet of 0.4 mm thickness, with density ρflexible sheet = 7855 kg/m3and Young’s modulus Eflexible sheet = 2×1011 N/m2. The flexible sheet is attached to a revolvable circular cylinder of aluminum. A rectangular stainless steel mass is attached to the other end of the sheet. The physical Figure 1: Structural model and dimensions in mm. 5 516 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer Figure 2: CFD grid around structure. dimensions and shape of the structure are illustrated in Figure 1. Both, the front cylinder and the rectangular mass, can be considered rigid. The structure is placed inside a vertical tunnel with a test section cross-sectional area of 180 mm ×240 mm and a length of 338 mm. The fluid is water at 25◦C, with kinematic fluid viscosity νfluid =0.97 ×10−6 m2/s and density ρfluid = 998 kg/m3. The bulk fluid velocity at the inlet is 0.68 m/s. The structure exposed to the incoming flow velocity oscillates around a mean position, and the flexible sheet deflects in the first mode. The detailed experimental setup, the measurement techniques, and the physical properties of the structure are described in [6]. 5 COMPUTATIONAL APPROACH The structural subproblem is solved employing the finite-element solver FEAP [21]. For the fluid domain, the finite-volume solver FASTEST [13] with block-structured body-fitted grids is used. The parallelization in FASTEST is achieved with domain decomposition and communication via MPI. The data transfer between the two codes and the interpolation on non-matching grid interfaces is performed via the interface coupling code MpCCI [9]. For details concerning the coupling algorithm see [17]. Both, the structural and fluid solver employ fully implicit second-order temporal discretizations. For the spatial discretization, the structural solver involves hexahedra elements with enhanced-strain formulation, whereas the fluid solver utilizes a second-order MUSCL [25] for the 2-d URANS simulation, second-order Central Differencing Scheme (CDS) for the LES and a blending between CDS and the GAMMA scheme [10] for the DDES. The blending between two schemes is done via a function introduced in [22], which ensures the calculation of fluxes with the CDS in LES regions and with the GAMMA scheme in RANS regions of 6 517 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer the DDES. Table 1: Number of CV and time step size Abbreviation Turbulence No. of Time step Averaged no. model CVs ∆t[s] of periods Sim-1 ζ−f0.37 ×1062.5×10−49 Sim-2 ζ−fDDES 12.0×1061.5×10−413 Sim-3 Dyn. Smag. 40.0×1061.5×10−47 Table 1 summarize the number of Control Volumes (CV), the time step sizes ∆t, and the number of motion cycles performed for averaging of the structural deflections. Figure 2 presents a view of the grid used for 2-d URANS (Sim-1), in x-y plane around the structure, where every 4th grid-line is shown. The grids are designed to have y+<1, for the first cell adjacent to solid walls. The time step sizes are lower bound by the artificial added mass effect [2]. The convergence of the coupled problem was observed to deteriorate, when reducing the time step size. The CFL number based on the time step sizes varied between 1.4 and 2.0. One reason for large variations in CFL number is the grid movement with structural deflections, where maximum CFL numbers are observed when the structural deflection or the velocity approaches a maximum. 6 RESULTS AND DISCUSSION The periodic motion of the structure is adequately predicted by the simulations, whereas quantitative comparison among the simulations and the experiment is based on the averaged structural deflections. The averaging is performed in time-phase as suggested in [6], after cycle-to-cycle variations of the end mass displacements reach a minimum. The number of motion periods averaged for each simulation are listed in Table 1. Figure 3 compares the absolute velocity contours of the experiment and Sim-1 at different phase angles. The arrangement of flow instabilities from the simulation is comparable with the experiment, despite the 2-d approach in Sim-1. Table 2 draws a quantitative comparison for the oscillation frequency of the structure fFSI, the end mass phase delay φshift and yextrema of the end mass normalized by the cylinder diameter. The quantities in Table 2 are time-phase averaged. The time-phase averaged cylinder rotation angle, and the end mass excursions are plotted in Figure 4a and Figure 4b, respectively. A slight asymmetry in the simulation data can be observed from Table 2 and Figure 4. This asymmetry is more noticeable in Sim-1 and Sim-3, where the motion cycles performed for averaging are less than that of Sim-2. The end mass displacement and the cylinder rotation from Sim-1 are underestimated, with yextrema of the end mass 20% lower than that of experimental values. The overdamping 7 518 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer 0o45o90o135o Figure 3: Phase-resolved contours of absolute velocity at different phase angles for Sim-1, comparison of experiment (top) and simulation (bottom). Table 2: Oscillation frequency, the end mass phase shift and yextrema of the end mass displacement fFSI[hz] φshift[deg] (uy)∗ max (uy)∗ min GL10 [6] 4.45 95 1.12 -1.11 Sim-1 4.58 109 0.85 -0.90 Sim-2 4.37 84 1.24 -1.25 Sim-3 4.37 83 1.25 -1.26 of the 3-d flow configuration in a 2-d flow simulation can explain the under-prediction of the the structural deflections. A study conducted by Breuer [1] for 2-d simulation of a circular cylinder at Re= 3900, reports the overdamping of turbulent fluid motion. Sim-2 and Sim-3 exhibit a close agreement in predicted values of the structural deflections. yextrema of the end mass displacement are about 13% higher than the experimental extrema of the end mass. A reason for this pronounced increase in the structural deflections could be the negligence of structural damping. A study of the laminar version of this test case [7] produced good agreement when simulating the Movement Induced Instability (MII), whereas large differences are observed for IIE where the structure oscillates in the first mode, as it is the case with this turbulent benchmark test case. The new test cases proposed in [4, 3] (also introduced in Section 1) study the effects of material damping in two different modes of the structural oscillation. The study depicts a higher importance of the material damping in the first mode of the structural oscillation, where the damping model significantly effects the structural deflections in numerical simulation. Nevertheless the material damping for the rubber (used in [4, 3]) would be higher than steel, and the premise that material damping is the cause of the over-prediction in Sim-2 and Sim-3, might not apply. Other possible reason for the differences between simulations 8 519 A. Ali, T. Reimann, D.C. Sternel and M. Sch¨afer 0 90 180 270 360 Phase angle [o] -40 -20 0 20 40 Rotation angle [ o ] GL10 Sim-1 Sim-2 Sim-3 (a) 0.02 0.04 0.06 0.08 x [m] -0.02 0 0.02 y [m] GL10 Sim-1 Sim-2 Sim-3 (b) Figure 4: (a)Cylinder rotation angle plotted against time phase angle. (b)Trailing edge coordinates. and the experiment could be the ignored side walls, taken as symmetry boundaries in the simulations to reduce the computational cost. The effect of side wall boundary layers are ignored on the assumption that the test section width in experiment is too high (about 8 times the diameter of cylinder) for the side wall boundary layers to have a significant effect. 7 CONCLUSIONS We have presented the results for a turbulent FSI benchmark test case. Three different turbulence modeling techniques have been studied. The 2-d URANS depicts a reasonable agreement with the experiment, regardless of the 3-d flow configuration. The excessive fluid damping in 2-d URANS is considered to be the cause of underrated structural deflections. The LES and the DDES simulation reproduce a close agreement between each other and an acceptable agreement with the experiment. The probable causes of overestimation of deflection in two simulations are also discussed. Further investigations with a variation of the structural material model are planned, as well as the simulation of the test cases [4, 3] with ζ−fmodel. 8 ACKNOWLEDGMENTS The authors gratefully acknowledge the Research Training Group 1344 of the German Research Foundation for the funding of this project. 9 520