Numerical investigation of flow around hairy flaps cylinder using FSI Capabilities

Flow around a cylinder has been extensively studied due to its practical importance in engineering, much attention has been devoted to drag reduction and vortex shedding suppression using active and passive control devices. A strong motivation for studying such a practical problem come from the fact that large amplitude of lift fluctuation and alternate vortex shedding are generally design concerns for engineers. Several numerical and experimental studies have been devoted for studying flow around a cylinder with a flexible plate attached to its centerline. This investigation has been known to be one of the most successful ways of controlling vortex shedding. Other active device for vortex shedding cycle reduction in order to minimize sound pressure is a cylinder equipped with hairy flaps. Experimental studies of air around a cylinder equipped with hairy flaps show that hairy flaps can reduce the wake deficit by modifying the shedding cycle behind the cylinder. For the model and simulation of such a coupling problem, fluid structure capabilities need to be performed. Fluid structure coupling problems can be solved using different solvers; a monolithic solver and a partitioned process. Monolithic process is a fully implicit method preserving energy at the fluid structure interface. However, its implementation is more complex when specific methods are required for both fluid and structure solvers. When efficient fluid and structure software packages are available, a partitioned procedure can be used in order to couple the two codes. The present work is devoted to simulation of fluid structure interaction problems and flow around thin flexible hairy flaps, using a partitioned procedure. One of the main problems encountered in the simulation is the automatic remeshing when the flaps come into contact and the fluid mesh between the flaps undergoes high mesh distortion. In this paper, numerical simulation has been performed and Strouhal number for two different Reynolds number Re=1.46 104 and 1.89 104 have been investigated. For both Reynolds number experimental data is available. For comparison with flow around a plain cylinder, numerical simulation was also performed and Strouhal number compared to experimental value. ___________________________________ *Corresponding Author: mhamed.souli@univ-lille1.fr 190 Numerical investigation of flow around hairy flaps cylinder using FSI Capabilities


INTRODUCTION
Fluid Structure Interaction (FSI) plays a crucial role in design for engineering and manufacturing.For a new product design several numerical and experimental tests need to be performed before getting to production process.These can be encountered in many engineering applications, where mechanical structures can be subjected to complex flow causing large deformation and damages.In Civil engineering, the design of fluid tanks storage require requires FSI simulation for the knowledge of sloshing frequencies and hydrodynamic pressure distribution on the structure.Empirical data presented in current tank seismic design codes as Eurocode, are based on simplified assumptions for the geometry and material tank properties and cannot be used for complex geometry.Fuel tanks storage may undergo different types of loading, including seismic loading, where the behaviour of storage tanks includes material nonlinearities, which are caused by material yielding.The ALE (Arbitrary Lagrangian Eulerian) formulation based on Finite Element analysis, presented in Ozdemir et al. [1], has been used for this application, taking into account material properties of the structure as well as the complex geometry of the tank.The formulation uses a moving mesh with a mesh velocity defined base on the structure boundary motion.These simulations can be very useful for engineers and designers to define appropriate material properties and shell thickness of the structure to be resistant under seismic loading.In previous research papers, Longatte et al. [2] developed and used FSI capabilities for nuclear industry, where mechanical structures like heat exchanger are subjected to complex flows causing possible vibrations and damages.
For a flow around a plain cylinder, several numerical and experimental investigations can be found in the literature.Due to practical importance in engineering of flow around a cylinder, much attention has been devoted to drag reduction and vortex shedding suppression using active and passive control devices.For drag reduction and vortex shedding suppression on flow past a cylinder, different researchers used different approaches such as splitter plate placed at upstream [3].
Numerical and experimental studies have been devoted studying flow around a cylinder with a flexible plate attached to its centerline.This investigation has been known to be one of the most successful ways of controlling vortex shedding.In [5] Kamps et al. used hairy flaps attached to the cylinder as a specific active device for vortex shedding cycle reduction in order to minimize sound pressure.Experimental studies of air around a cylinder equipped with hairy flaps, as described in Figure 1, show that hairy flaps can reduce the wake deficit by modifying the shedding cycle behind the cylinder.For flow around a cylinder equipped with hairy flaps, few wind tunnel experimental investigations were conducted for studying the effects of flaps motion on the reduction of wake and vortex shedding cycle.To the best of our knowledge and by searching in the literature through published papers, there is no reported work on numerical models, and no full numerical simulation has been carried out for this specific application, using full finite elements analysis of deformable thin flaps.In [4] a homogenization approach for numerical simulation was used to model the fluid flow around a circular cylinder partially coated with hair.
From experimental investigation, carried out by Kamps et al [5], results show that, above a certain Reynolds number, the hairy flaps lead to a jump in the vortex shedding frequency.This phenomenon is similarly observed in the water flow experiments as a jump in the nondimensional Strouhal number that is related to the change of the shedding cycle.Numerical simulation of this specific problem, including complexities due to presence of very thin flexible hairy flaps, is addressed in this paper for two different Reynolds numbers.
The complexity of the numerical simulation is due to the fact that the flaps get too close or into contact with each other.An efficient contact algorithm needs to be developed to compute contact forces that need to be added to hydrodynamic forces from the fluid flow.The presence of highly deformable flaps attached to the cylinder lead to a complex fluid structure interaction problem.For high flaps deformation, the solver can automatically re-mesh the fluid domain based on the nodes located at the new fluid structure interface with respect to Delaunay criteria.This paper is devoted to FSI simulation of flow around a cylinder equipped with hairy flaps.Structures like cylinder with many flaps, submitted to hydrodynamic pressure cross flows, requires fluid and structure to be in equilibrium at each time step, which can be performed by several ways.A first method consists in solving fluid and structure equations in a single linear or non-linear system using monolithic algorithm.This is a strong coupling process ensuring energy conservation of the full-coupled fluid structure system.However, this approach is often hard to set up for complex applications as it requires significant developments in fluid and structure solvers.This difficulty can be overcome by using a coupling procedure where fluid and structure are solved separately, the fluid is solved using an implicit or semi-implicit method, and the structure is solved using implicit time integration.In this paper a numerical simulation using the state of the art fluid structure capabilities is performed to simulate laminar and turbulent flow around a cylinder equipped with thin flexible flaps.A staggered coupling method approach is performed to solve a fluid-structure coupling problem.In the first chapter, full Navier-Stokes equation for the fluid with projection method is described.In a second part, structure equilibrium equations for nonlinear Mooney Rivlin constitutive model for flaps material is presented.In part 3 of the paper the coupling algorithm with the remeshing process of the fluid domain is developed.Part 4 of the paper is devoted to problem description and validation of the numerical results to experimental data presented in [5].For numerical modeling, we present numerical results and Strouhal number for two different Reynolds number Re=1.46 10 4 and 1.89 10 4 .For both Reynolds number experimental data is available.As reported in [5], for these Reynolds numbers, the tonal peak due to vortex shedding seems to be shifted to slightly higher frequencies, indicating that the tripping tape merely corresponds to a small increase of the effective cylinder diameter in this Reynolds number range.

EQUATION OF FLUID MOTION
For incompressible fluid, conservation of mass and momentum in the ALE frame of reference are represented by the Navier-Stokes equations: Mass conservation: Where is ρ is the density, σ the stress tensor.For Newtonian fluids the stress tensor σ may be defined as: . d p I where τ is the shear stress from the constitutive model, and p the dynamic pressure.w is the velocity of the reference geometry.Note that the Eulerian equations commonly used in fluid mechanics are derived by assuming that the velocity of the reference configuration is zero, 0 w = .In the Lagrangian formulation, mainly used in solid Mechanics, the relative velocity ( ) v w − between the material and the reference configuration is zero, therefore the material velocity v w = .The term in the relative velocity in (3) is usually referred to as the advective term, and accounts for the transport of material past the mesh.It is the additional term in the equations that makes solving the ALE equations much more difficult numerically than the Lagrangian equations, where the relative velocity is zero.
Using equations ( 2) and (3) with implicit time integration, the problem to be solve at each iteration is the following: It is well known in computational fluid dynamics, that the main difficulty arising in the numerical solution of the convection-diffusion equations is due to the presence of the convective term in equation (2) (3), which leads to numerical oscillations for high and moderate Reynolds numbers.The standard Galerkin method leads to no physical spatial oscillations when applied to the high convective case.There are different ways to stabilize the nonlinear term in the Navier Stokes equations ( 1) and (2).To preclude such anomalies, the most popular method being the use of upwind differencing on the convective term via Petrov-Galerkin methods.
Although these methods are precise and stable, we will use a 'split' method which is a simple mean to obtain a robust and effective formulation.The projection method, introduced initially by Chorin and Temam [6] and proposed by Gresho [7] is used.This velocity does not satisfy When solving equation ( 4) for the velocity 1 n v + , the incompressibility constraint given by Equation ( 1) is not satisfied.For the velocity to satisfy free divergence, the present scheme consists of three steps until the final velocity is achieved In the first step, a predictor velocity u is computed.Thus, in the second step the velocity u is projected into a space of divergence free vector field and a Poisson equation of pressure is obtained.This pressure will correct the velocity to get a final step.These algorithms are widely used in the previous papers and described in detail by authors of this paper in Idelson et al [8]. the following steps summarize the development in [9] for solving Navier Stokes equations ( 1) and (2).Each iteration inside the time step will involve three phases, each of them related to different steps of the time integration scheme: We first solve the momentum equation ( 5) for an intermediate velocity * v .Since the intermediate velocity * v is not divergence free, a projection type method is performed using pressure correction to solve Navier-Stokes equations ( 2) and (3).
-Phase 2 Poisson equation for Pressure The second step consists in deriving a Poisson equation for the new pressure 1 n P + .In fact, by taking the divergence of the equation (2) and using the incompressibility condition (1), the pressure is used as Lagrange multiplier to satisfy incompressibility condition.
-Projection.Compute velocity As the velocity * v does not yet satisfy the incompressibility condition (1), it is projected on a divergence free space to get an adequate approximation of the velocity.This is obtained from: . ( )

STRUCTURE EQUILIBRIUM EQUATIONS
The cylinder is equipped with hairy flaps made of a hyper-elastic rubber material.The equilibrium equation for the structure is a momentum equation similar to the one used in fluid mechanics with no advection term, since the Lagrangian framework is used for the structure, the velocity of the structure s u satisfies the following momentum conservation equation: ( ) where s ρ is the structure density, and s σ the Cauchy Stress Tensor for the structure material.
An hyper-elastic material is defined as an elastic material in which the stress at each point can be derived from a scalar function W(F) called the strain energy function.In order to meet the requirements of objectivity, the strain energy function must be invariant under changes in the observer frame of reference.It is well known that the Cauchy-Green deformation tensor is invariant under changes in the observer frame of reference.Thus, if the strain energy function can be written as a function of the right Cauchy-Green tensor s σ .The general stress- strain relationship is given by the formula (9): where S is the Piola-Kirchhoff stress tensor.The different models in the literature define the strain energy as a function of the strain field.We consider the Rivlin's theory, Rivlin (1948).Rivlin's theory of isotropic materials describes the energy as a function of the three Cauchy strain invariants I1, I2, and I3.Assuming material incompressibility, the stress can be obtained from the following Mooney-Rivlin form (10): The use of two terms in the series is sufficient to describe the elastic modulus in both the uniaxial and biaxial deformation modes.

(
) ( ) Here, c1 and c2 are the Mooney-Rivlin material coefficients, I1 and I2 are the invariants of the deviatoric part of the right Cauchy-Green deformation tensor, and I3=det(F) is the Jacobian of the deformation, K is the bulk modulus of the rubber material.When C2=0, this model reduces to an uncoupled version of the incompressible neo-Hookean constitutive model.

FLUID STRUCTURE COUPLING ALGORITHMS
In CFD simulation we consider rigid structure and solves for fluid pressure and velocity around or inside the structure.This can be acceptable if the structure is at its steady state.For general fluid structure applications, the nature of the flow can be affected by the deformation of the structure, which needs to be taken into consideration, mainly at the coupling interface.
At the fluid structure interface, the coupling method needs to satisfy the kinematic and dynamic conditions, continuity of normal velocity and normal stress. . .
The kinematic interface conditions (12) state that no interpenetration of fluid at the fluid structure interface, solid and fluid geometry must match at the interface but not necessarily the nodes.The dynamic interface condition (13) requires force equilibrium at the interface.
Structures submitted to hydrodynamic loading, require fluid and structure equilibrium equations to be solved at the same time.This can be performed by several ways.A first method consists in solving fluid and structure equations in a single linear or non-linear system using monolithic algorithm.This is a strong coupling process ensuring energy conservation of the full-coupled fluid structure system.However, this approach is often hard to set up for industrial applications as it requires significant developments in fluid and structure solvers.Previous research papers and mathematical analysis show that the full non-linear system to be solved is ill-conditioned with always non-defined positive matrices.This difficulty can be overcome by using a coupling procedure where the fluid is solved using an implicit or semiimplicit method, and the structure is solved using implicit time integration.This method is easier to set up and it allows independent model developments in both fluid and structure solvers.The procedure is iterative, and each iteration is made of three steps: This approach has a great flexibility due to its modularity.The partitioned procedure may rely on several kinds of time coupling schemes, explicit or implicit.
Most classical fluid structure coupling methods rely on the structure response from experimental measurements associated with convenient data processes.For applications where structure deformations are not important and are restricted to added mass effects, no strong coupling between structure and fluid motions can be considered, and the structure displacement is not supposed to strongly affect flow patterns.Thus, it is possible to solve flow and structure problems separately using explicit coupling algorithms.In Explicit coupling, fluid and structure computations are staggered in time.All explicit coupling algorithms inherently introduce energy because it is impossible to predict correctly the structural displacement inducing correct forces when solving the fluid problem.
On the coupling boundary the fluid velocity and the solid velocity should converge to the same value at each time step.An iterative process is needed at each time step for velocity convergence at the interface, Convergence occurs when the difference between velocities of successive iteration steps is less than the acceptable error.The iterative process is described in Figure 2. From a computational point of view, the problem to be solved is domain decomposition is used for parallel computing using shared of distributed memory.
In the SMP technology multiple CPUs share same memory.Since CPUs are accessing shared memory simultaneously, data supplied from memory to CPUs is likely to be slow, when the number of CPUs increases, which limits the number of CPUs that can be used efficiently.To overcome this problem, the MPP (Massive Parallel Computing, using distributed memory), version of crash codes has been developed.In MPP, each small group of CPUs has its own memory, and CPUs are interconnected via a high-speed networking system.When using the MPP version, prior to computer processing, domain decomposition needs to be performed to divide the problem.In domain decomposition the model is decomposed into subdomains, each subdomain is processed by a CPU and uses the memory dedicated to that CPU.To solve the fluid structure coupling where the fluid is meshed using few million elements, the MPP version of LSDYNA with 64 processors is performed for solving the problem is a reasonable time.Each problem takes 2 to 3 days, using massively parallel computing.

Moving mesh algorithm
The automatic remeshing algorithm is extensively described in the literature.In [9] Facundo et al. developed robust remeshing algorithms in LSDYNA, that has been validated and used with great success for several industrial applications.For small structure deformation a classical ALE method is performed based on mesh velocity solved using Laplacian equation see Souli et al. [10].
For high deformation, the remeshing process is a necessary step to remedy element distortions occurring during large structure motion.For high deformation of the flaps behind the cylinder, new parts of their boundaries come into contact which makes it difficult to fit a new mesh between flaps.In order for mesh solver to fit properly the new geometry with flaps into contact a new contact option needs to be developed using a gap as an offset to contact detection.A structural mesh would have been ideal for remeshing the new fluid domain, however, tetrahedrons are more flexible for morphing geometry in a dynamic environment.
The first step of remeshing is only to decide whether to re-mesh automatically or not.The decision is taken if one of the following is observed: 1-One element is overly distorted (the Jacobian of the transformation between the regular and the real element must be strictly positive).2-One edge of the boundary of the mesh is overly distorted (quadratic triangles are used and then the boundary edges are also quadratic).The curvature of an edge must be smaller than a prescribed value given by the user.3-The boundary node distribution does not permit a precise description of the contact area.4-Next, we define the geometric precision of a mesh in order to control and create the set of boundary nodes and boundary edges, i.e. the boundary mesh.

Structure Contact Algorithm
For the last two decades a large amount of research and efforts have been devoted to contact problems which reveals the importance of the phenomena in computational mechanics.
Various numerical methods and computational techniques have been proposed for different types of contact, including kinematic and penalty contacts.In this paper penalty contact is used between different parts edges.
The basic idea of penalty contact is to track the relative penetration at the contact interface.In contact algorithms, one surface is designated as a slave surface, and the second as a master surface.The nodes lying on both surfaces are also called slave and master nodes respectively.A contact force is computed proportional to the penetration vector, the amount the constraint is violated.Penalty method imposes a resisting force to the slave node, proportional to its penetration through the master segment, as shown in Figure 3 describing the contact process.This force is applied to both the slave node and the nodes of the master segment in opposite directions to satisfy equilibrium.The main difficulty in the coupling problem comes from the evaluation of the stiffness coefficient K in equation ( 14).The stiffness value is problem dependent, a good value for the stiffness should reduce energy interface in order to satisfy total energy conservation and prevent inter penetration.In this paper, the stiffness value is described in Hallquist [12].There are different methods to compute the spring stiffness K in the literature, the commonly used one is the spring stiffness k given by ( 14) in term of the bulk modulus K of the fluid element, the volume V of the fluid element and the average area A of the structure element connected to the structure node.testing but rather reduce the number of costly tests.The problem to be solved numerically consists to an experimental investigation performed by Kamps et al. [5] for acoustic and noise measurement of air flow around a cylinder equipped with hairy flaps described in Figure 1.
The structure shown in Figure 1 is an example of the complex geometry that needs to be simulated submitted to fluid loading.The goal of this experimental investigation is to reduce the frequency of vortex shedding and thus noise intensity.From this experimental investigation, the Strouhal number is measured for different fluid Reynolds numbers.These are the data we use for the validation of the numerical simulation.
As developed in [5], and described in Figure 1, the problem consists of an aerodynamically closed test section of approximately 1.5 m length that is mounted to the rectangular nozzle of the open jet wind tunnel, which has a contraction ratio of about 31.2 and an exit area of 0.23 m.The two lateral sides of the test section are covered by tensioned Kevlar that allows sound waves to propagate from the flow region inside the test section to the outside.The upper and lower walls of the test section are made from thick acrylic glass.Therefore, the test section provides a nearly two-dimensional flow.The cylinders are positioned approximately 0.2 m downstream from the nozzle and centered between the two Kevlar windows.
Silicone rubber material is used for flaps structure see [5].Silicone based materials offers a broad scope for adjusting the material properties, Young Modulus can be tuned depending on the mixing ratio of silicon to its curing agent as well as on the curing procedure.Furthermore, the properties of these materials do not change noticeably under temperature differences which facilitate the reproducibility of experiments.In the experimental investigation, several flap rows were individually cast as a ring with eight parallel flap structures of silicone rubber Elastosil RT 601, presenting a Young Modulus of 1.2 MPa.The thickness of the flaps is around 0.3 mm, and their length and spacing are 9 mm and 3.6 mm, respectively (see Figure 5).
From experimental studies, and the Reynolds number used in experimental tests, the flow can be turbulent.For turbulence model, the SST k-w turbulence model is a good candidate for this application as it is directly usable as a Low-Re turbulence model without any need for extra damping functions and it has shown good behavior in a wide range of industrial applications featuring strong adverse pressure gradients and separating flows [13], current implementation follows the guidelines of Mierka [14] and uses Dirichlet boundary conditions at the wall for k (k=0) and w (w=6.viscosity/((0.09) 1/2 .y 2 ) with y the distance to the nearest fluid volume node).6. NUMERICAL SETUP AND NUMERICAL Several new computational techniques have been used in order to improve the performance of numerical methods in term of CPU time, as domain decomposition and parallel computing.When solving complex coupling problems involving fluid structure coupling, turbulence model, contact between different parts of the structure, high performance computing can be very useful to solve problems in an acceptable engineering computing time.
In the simulation, material parameters and structure dimensions are the ones used for experimental tests and described in [5].Several experimental tests have been performed using different Reynolds number.For our simulation two different Reynolds number Re=1.38 10 4 and 1.89 10 4 have been simulated, using different mesh refinements.For these Reynolds number the average inflow velocity v = 7m/s and 10m/s, with a Poiseuille inflow velocity profile.The viscosity of the air is set In the flow direction, which corresponds to the horizontal direction in Figure 5, the fluid domain is five times larger than the diameter of the cylinder and 5 times in the transverse direction, which corresponds to the vertical direction.The model corresponds to an appropriate distance so that external boundary conditions, will not have significant influence on the flow around the cylinder and the flaps.
In order to capture contact between the flap tips, a finer mesh is used around the structure to resolve boundary layer issues as shown in Figure 6.As the structure experience high deformation due to the highly flexible material as shown in Figure 7, a remeshing process is required every at each time step, which increases the CPU time of the numerical simulation.High Performance Computing is performed using 64 processors with a good scaling, each simulation takes around 60 hours with 64 processors to reach termination time.
Lift force of the cylinder oscillates with the vortex shedding frequency, by using a classical Fourier transform of time history of the lift force described in Figures 11,13        As observed in Figures 9 and 10, the structure of the mean flow velocity with a plain rigid cylinder reaches stationary state with a fixed frequency of the vortex shedding as shown in Figure 15, where time history of the lift force and the frequency of vortex shedding are plotted in Figures 15 and 16.For flexible hairy cylinder flow recirculation between the flaps takes longer time to reach a stationary state.In Figures 11 and 12, the lift force and frequency of vortex shedding are plotted for Re=1.46 10 4 .For Re=1.89 10 4 , we plot time history of the lift force and the frequency of the vortex shedding, The Strouhal number is computed based on these frequencies and compared to experimental data presented in Tables 1 and 2.

CONCLUSION
In this paper, a numerical investigation of flow around a cylinder with flexible hairy-flaps at the aft-part of the cylinder is carried out.The complex model shows the capability of the software when different numerical algorithms are combined, fluid structure coupling, turbulence modeling, and structural contact between nonlinear materials.To the best of our knowledge and by searching in the literature through literature, there is no reported work on numerical simulation including full model with flexible hairy flaps, contact between flaps and coupling fluid structure.The software used in the paper is a Multiphysics software package using two robust solvers, a CFD with turbulence models and a nonlinear structure solver with efficient contact algorithms.The coupling method developed in the paper is based on the contact algorithms used in computational structural mechanics.The coupling method has been used successfully to solve several industrial applications in aerospace and automotive industry.The main difficulty in the present simulation is the meshing process when the flaps get into contact and the fluid volume between flaps is dramatically reduced.In this specific situation contact needs to be detected between flaps using an appropriate offset which allows the presence of a small fluid volume.The presence of the flaps modifies the Aerodynamics performances quantified by the Strouhal number, as well as the drag and the maximum lift.
In term of Strouhal number, numerical results are in good agreement with experimental data provided in [5]

Figure 1 :
Figure 1: Photograph of the cylinders used in the present study left aft body fur right hairy flaps [5] hairy flaps for simulation

1 -
computation of fluid forces acting on the structure by solving a CFD problem 2-estimation of structure displacement and velocity induced by fluid forces, solving a structure dynamic problem.3-actualization of the fluid domain according to the structure wall motion.

Figure 4 :
Figure 4: Description and dimensions of the cylinder with flaps

Figure 5 :
Figure 5: Experimental Setup from [5] a nonlinear material law for the flaps, in this case the hyper-elastic Mooney Rivlin material, this corresponds to a value of the coefficients 1 and 15, the first frequency of the Fourier transform corresponds to the frequency of the vortex shedding.Using the definition of the Strouhal number in eq.(15); frequency, D the diameter of the cylinder and U the mean inflow velocity.

Figure 6 :
Figure 6: Flow Velocity and high structure deformation at time t=1s

Figure 7 :Figure 8 :
Figure 7: Flow Velocity and high structure deformation at time t=1s

Figure 16 :
Figure 16: Fourier Transform of lift for plain cylinder Re=14600 . From experimental tests, for Reynolds number Re=1.46 10 4 and 1.89 10 4 , the Strouhal number increases to a value St=0.23 and St=0.26 for a cylinder with attached hairy flaps; however, for flow around a plain cylinder, Strouhal number for Re=1.46 104 is St=0.210.This increase of Strouhal number observed in experimental tests is in good correlation with our numerical results.

Table 1
presents experimental and numerical Strouhal numbers for the 2 Reynolds number used in the simulation for cylinder with hairy flaps.Table2presents experimental and numerical Strouhal number for Re=1.46 10 4 used in the simulation with a plain cylinder, as shown in Figure10

Table 1 :
Cylinder with hairy flaps experimental and numerical Strouhal number.

Table 2 :
Plain cylinder experimental and numerical Strouhal number.