Preliminary Study on Mesh Stiffness Models for Fluid-structure Interaction Problems

. One of the challenges in modern computational engineering is the simulation of fluid-structure interaction (FSI) phenomena where one of the crucial issues in the multi-physics simulation is the choice of stiffness model for mesh deformation. This paper focuses on the application of iteratively implicit coupling procedure on two transient FSI cases of vortex induced-vibration (VIV) that manifest oscillating flexible structures. The aim is to study various mesh stiffness models in the Laplace equation of diffusion employed within the arbitrary Lagrangian-Eulerian (ALE) methodology to handle the moving mesh. In the first case where a laminar flow interacted with a flexible splitter, it was demonstrated that a near FSI boundaries increased-stiffness model prevails to manage a large deformation of the moving structure as compared to a near volume increased-stiffness model. However, the potential technique could not be exploited to the second FSI configuration, where the effect of the turbulence of flow was included. It was found that the mesh topology near the FSI interface was collapsed. Instead of utilizing the same approach, a mesh stiffness based on a wall distance was found to be auspicious. Thus, the mesh stiffness model in the FSI simulation is case-dependent.


Introduction
Interaction between external or internal flow and deforming or moving structure, so-called as fluid-structure interaction (FSI), is omnipresent in a multitude of engineering systems in aero-acoustics, biomedical engineering, civil engineering, mechanical engineering such as reed valve (in musical instrument and refrigeration compressor), heart valve, Tacoma bridge, and wind turbine. To simulate the FSI devices, in principle various discretization and coupling techniques available can be exercised. Within the discretization strategy, one can distinguish the arbitrary Lagrangian Eulerian (ALE) formulation, the Immersed Boundary (IB) method, and the fictitious domain. Furthermore, depending on the level of deformation in the moving structure, either the monolithic approach, the explicit-partitioned method or the implicit-partitioned scheme can be applied for the coupling technique.
In this paper, laminar and turbulent FSI cases of vortex induced-vibration (VIV) demonstrating oscillating pliable structures are numerically studied [1,2]. Such FSI phenomena play an important role in numerous engineering systems in which their occurrences can impact on the life span of the devices primarily determined by the reciprocal interaction between the neighbouring flow and the maximum extent of structural vibrations. The laminar-FSI configuration was to demonstrate their numerical solution strategy for the FSI simulation. The geometry consists of a fixed square body with an elastic splitter attached at the back of the body and is submerged in an incompressible laminar flow. After Hübner et al. [1], De Nayer et al. [2], Dettmer [3], Lee and You [4], De Nayer [Error! Reference source not found.] studied the benchmark for the evaluation of their numerical solution approaches. In the study of Hübner et al. [1] the monolithic procedure and a decreased flow velocity at inlet were used.
Moreover, a variant of the IB method proposed by Mittal et al. [6,7] was applied by Lee and You [4] in the first FSI problem. For the latter case, the turbulent-FSI benchmark is composed of a fixed circular cylinder with an attached rubber at the rear of the cylinder and is immersed in an incompressible turbulent flow at a sub-critical Reynolds number. A prior study by De Nayer et al. [Error! Reference source not found.] utilized a Smagorinsky subgrid scale model for the turbulence simulation over the test case. The effect of turbulence on the mean flow is considered to be sufficient, this study thereby employs an industrial turbulence model for the simulation of the FSI benchmark.
In contrast with Hübner et al. [1], the iteratively staggered coupling is exploited in this study. The procedure does not require more resource and expertise from a software engineering point of view as in the monolithic approach, according to Hou et al. [8]. In addition, the coupling strategy is usually more robust than the explicit-partitioned method and demonstrates an improved convergence. Furthermore, the ALE methodology to be used on boundary-fitted grids is utilized in this study; the similar aspect as in some studies [1][2][3]5]. The ALE method which potentially can manage complicated mesh movements in the fluid domain is not based on an interface capturing approach the basis of the IB method. In the ALE approach, on the other hand, the fluid mesh is moved in every time step thus making it be more straightforward than the IB method with respect to the treatment of boundary conditions and their ramifications on the accuracy and conservation properties of the numerical scheme [9]. For a more detailed description of the IB method, the readers are referred to Mittal's research [6,7]. Additionally, The ALE approach is not limited to high Reynolds number flows, which is the weakness of the fictitious domain counterpart of Baaijens.
A pivotal issue within the ALE formulation concerns with a reliable technique to handle the moving mesh, especially for the turbulent flow where the mesh motion can be intricate owing to the turbulence effect. In this situation, the strategy should be able to preserve the quality of mesh while the topology of the moving mesh is changing in time and at the same time to capture the motion of structure accurately. From literatures, various methods are available for the treatment of the moving mesh update such as mesh stiffness models in the Laplace equations of diffusion, interpolation techniques such as transfinite interpolation (TFI) of Gordon and Hall as used by Kneißl et al. [10] for example, elliptic method of Spekreijse, and pseudo-structure approximation based-model of De Nayer [Error! Reference source not found.]. Throughout this paper, various mesh stiffness models contained in the Laplace equation of diffusion are investigated to study their potentials for the FSI simulations.

Arbitrary Lagrangian-Eulerian (ALE) formulation
In the case of FSI, the governing equations of fluid and structure dynamics written respectively in the Eulerian and Lagrangian frameworks require modifications. The alterations are achieved by means of the ALE formulation. As the fluid mesh deforms in time because of the interrelations between forces passed by the fluid to the structure and displacement sent back by the solid to the fluid, this can cause the destruction of the structural mesh when the large grid distortions in the solid (provoked by the large deformations from the fluid) occur. The Lagrangian description thus is not ideal for FSI problems. Conversely, when the Eulerian approach is applied to the solid for the FSI cases its accuracy is lessened given that the basis of the Eulerian system is to monitor the motion at a fixed location in space. The ALE method was used in pioneering works, where its capability to handle various FSI problems with great distortions of the continuum was proven.
Considering the fluid mesh deformation, in the ALE approach an observer can move arbitrarily, meaning that the viewer is neither at a fixed position in space nor moves together with a material point. To implement this idea, the fluid integral conservation equations, i.e. the conservation equations of mass and momentum are modified by applying the Leibnitz rule. This results in a set of new conservation equations in their ALE formulations. (2) In equations (1) and (2), the velocity of time-dependent surface of the control volume (CV) or the velocity of the moving grid g u is introduced to transfer the fluid velocity into a moving reference system. It is obvious that when the grid velocity 0  g u Equation (1) and Equation (2) will turn into the Navier-Stokes equations in the Eulerian description and into the Lagrangian formulation when the grid velocity u u g  . Moreover, a relative velocity g u u  for the calculation of convective fluxes at the cell faces of CV is used in the equations. As the faces of CV moves in the FSI simulation, this introduces additional mass fluxes, which implies that the mass balance is not necessarily ensured. Therefore, the space conservation law (SCL) must be imposed to avoid the mass imbalance in the modified conservation equations, which reads as: Equation (3) In the FSI simulation, the transfer of the forces f from the fluid part and of the displacement δ from the solid part takes place on the interface.

The Laplace equation of diffusion
) (x  in equation (6) can be defined for near boundary increased-stiffness, near volume increased-stiffness, and independent variable based-stiffness models such as wall distance. The expression of the stiffness models is given in equations (7), (8), and (9) for the near boundary increased-stiffness model, the near volume increased-stiffness model, and the wall distance based stiffness model, respectively.
In the equations, Lref is the reference length, d is the nearest boundary, Vref is the reference volume, V is the size of CV, and dw is the wall distance. Cstiff is the stiffness model exponent of which its value can be adjusted.

Implicit-partitioned approach
The implicit partitioned scheme to couple a structural solver and fluid solver is illustrated in Figure 1 adapted from Du [1111]. As can be seen in Figure 1, the interfacial boundary conditions are exchanged several times between the fluid and the structural computations within a global time step. In this study, only forces and displacements are exchanged on the FSI interface where each displacement and force transfer signifies the inception of a new FSI or staggered iteration. In each FSI iteration, the inner fluid and structural iterations based on the updated boundary conditions and the previously determined solution are carried out until their individual convergence criteria are reached. Additionally, to accelerate and stabilize the convergence of the FSI predictions the implicit coupling between both solvers requires an underrelaxation of the transferred variables, i.e. underrelaxation of the force and displacement components.

Transport equations of turbulence variables
The k- Shear Stress Transport (SST) of Menter is used in this study. This industrial turbulence model have two transport equations of turbulence variables, including turbulence kinetic energy k [J kg -1 ] and specific dissipation rate [1 s -1 ]. The k-and turbulence transport equations are given as: In addition, a new definition of the turbulence eddy viscosity vt is given as in Equation (9), where the shear strain rate tensor S replaces the vorticity tensor  .

Numerical setups
The computational domains of two FSI problems being studied is described in   In the two cases, inlet, outlet, slip walls, symmetries at the lateral sides, no-slip walls at the bodies, and FSI interfaces at the flexible plates were prescribed for boundary conditions of the fluid parts. In the structural domains, boundary conditions were set to the Dirichlet conditions for the upstream surfaces of the splitters attached on the rear of the cylinders to force zero degree of freedom, the Neumann conditions for the FSI interface surfaces (upper, lower, and end surfaces) of the elastic parts, and frictionless supports for the two lateral sides of the plates. For the laminar FSI configuration, an inlet velocity of 31.5 cm s -1 was specified, corresponding to a Reynolds number of 204 while an inlet velocity of 1.385 cm s -1 was defined for the turbulent FSI benchmark giving a sub-critical Reynolds number of 30 470. Inlet turbulence parameters were also introduced for the turbulent FSI case.
The FSI computations were performed with CFX flow solver, Mechanical APDL structural solver, and MFX multiphysics solver of ANSYS 15.0.7. The MFX was to couple the fluid and structural solvers with the implicit coupling algorithm illustrated in Figure 1. In this preliminary study, a timestep size t of 0.001 s was chosen for the laminar FSI configuration following one of Hübner et al. [1] and the turbulent FSI benchmark. The chosen timestep size was the same value with one of Hübner et al. [1] and a half of the change value after time sensitivity studies. The fluid meshes for the laminar and turbulent FSI cases were around 12 000 and 54 000 control volumes, respectively where unstructured meshes with one cell in the spanwise direction were exploited. For this preliminary study, the fluid mesh resolutions consider the sizes used by Hübner et al. [1] and Ali [12] for the first and second cases, respectively. Structured meshes of 120 and 400 second higher orderbrick elements with aspect ratios of around 1 were designed for the flexible plates in the laminar and the turbulent FSI problems, respectively. The structural mesh for the first FSI configuration considers the size used by Hübner et al. [1]. Nevertheless, the structural mesh resolution for the second FSI geometry was finer than one. Near FSI boundaries and near volume increased-stiffness models with the stiffness exponent of 2 were utilized for the first FSI benchmark of Hübner et al [1] while the near boundaries increased-model and the wall distance based-stiffness equation with the exponent of 2 × 10 6 were applied for the second FSI configuration of De Nayer et al. [2]. Figure 4. gives the result of the fluid mesh deformation modeled with the near small volume increased-stiffness captured at a sufficiently long period of simulation where the large deformation of the flexible part has appeared. As shown in the figure (see the right part for a clearer view), it was demonstrated that the given stiffness model is not safe to continue with larger deformations in the elastic structure. It is seen at the upper trailing edges of the oscillating plate there were several mesh elements with potentially folded topologies to occur. This situation is not expected as it can fail the FSI simulation. Orthogonal angle, expansion factor, and aspect ratio of the mesh can be monitored at every coupling iteration.  Unlike the first situation, the near FSI boundaries increased-stiffness model was proven to be more prospective than the first approach to manage the fluid mesh movement with the large deformation. It is seen that there were no collapsed mesh topologies around the trailing edges of the moving part. This is shown in Figure 5 (see the right part). Therefore, the near FSI boundaries increased-stiffness model was employed instead of the first mesh stiffness formulation and the simulation was extended for a longer time to allow a quasiperiodic oscillating motion of the elastic splitter with larger deformations. Figure 6 displays the results of larger deforming mesh in the fluid due to the larger deformations of the moving structure already in the quasi-periodic stage of the oscillation. From the figure, it was found that the near boundary increased-stiffness method succeeds to control the larger mesh deformation with no entangled mesh elements around the trailing sides of the moving body as shown in the right pictures. To this FSI problem, Hübner et al [1], Lee and You [4], and De Nayer [Error! Reference source not found.] successfully exercised various moving mesh strategies, i.e. geometric balance equation based mesh motion algorithm, artificial elastic continuum scheme with element-wise constant stiffness distribution, pseudo-structure approximation based remeshing module, and overlay mesh motion procedure of Campbell [9] based on the Laplace equation, respectively. Further discussions on the mesh update techniques can be found from Dettmer [3].  . demonstrates the result of the mesh deformation modelled with the near boundaries increased-stiffness approach for the FSI simulation of the second case. It is seen that the same mesh stiffness model as applied for the first FSI problem failed to preserve the mesh quality near the FSI interface when the large deformation caused by a turbulent flow occured. The mesh topologies close to the FSI boundaries were collapsed. Such a situation is not desired for correct calculations of the wall shear stress which affects accurate structural oscillations as the grids near the surface of the moving splitter did not satisfy the wall y + requirement for the chosen turbulence model. As a solution of this problem, a mesh stiffness scheme based on the wall distance from equation (6) was applied. The mesh stiffness model which is a hyperbolic approach to control the mesh stiffness in the fluid domain depending on the wall distance dw gives an exponential decrease of the mesh stiffness as the distance from the nearest boundary increases. The potential of the wall distance based-stiffness approach is shown in Figure 8. It is demonstrated that the chosen mesh stiffness model now can protect the grids near the wall of the moving structure from the violation of the y + requirement.

Conclusion
Various mesh stiffness models in the Laplace equations of diffusion have been studied through the FSI simulation of two vortex induced vibration (VIV) test cases available in the literature. From this preliminary study, it is demonstrated that the choice of prospective mesh stiffness models is case-dependent. Du also confirmed that the choice depends on mesh type and the geometry as there is no universal optimal mesh stiffness model [11]. Numerous mesh movement strategies are nowadays available where one can exploit straightforward the mesh update techniques provided in various multiphysics codes. Moreover, as instability inherently exists when the implicit partitioned method is used on the FSI problems with the incompressible fluids and the slender structures strong constant under-relaxation factors for the transferred forces at the FSI interfaces were used to stabilize the FSI simulations. The instability is also known as artificial added mass effect formulated by Förster et al. [13].