Fixed-Point Fluid structure interaction analysis BASED ON geometrically exact approach

The fluid structure interaction analysis for structures exhibiting large deformations is carried out by using a strong coupling method, in which a fixed point method with Aitken’s dynamic relaxation is employed to accelerate convergence of the coupling iteration, and geometrically exact beam approach initiated by Simo is adopted to simulate the nonlinear flexible beam models. An improved implicit time integration algorithm is given to improve the computation accuracy of structural dynamics. To verify the validity of the fixed-point method in the compressible flows which is usually used in incompressible fluid, it is applied for flutter analysis of AGARD 445.6 wing in the transonic regime. The case of flow-induced vibration of a flexible beam demonstrates that the approach based on geometrically exact beam theory is suitable for the fluid structure interaction analysis and the fixed-point method with Aitken’s relaxation is of great efficiency and robustness in the FSI computation.

www.nature.com/scientificreports www.nature.com/scientificreports/ characteristics of the structures with large displacement. Block-Gauss-Seidel or fixed-point method 7-10 is a strong coupling algorithm and commonly applied to FSI coupling between incompressible flow and structures with large deformation. It can be regarded as an extension of the loose coupling method. For each time step the computation of the two domains and the communication of interface motions and traction forces are repeated until interface displacements converge to the allowable tolerance to ensure the conservations of kinematics and dynamics. Some relaxation means can be employed because the convergence might be slow if the coupling of fluid and structure is strong 10 , such as flexible wings. Among those means, Aitken's dynamic relaxation with robust stability and easy implementation is suggested in ref. 1 .
The geometrically exact intrinsic beam model initiated by Hodges 11 is widely adopted to compute the goemetric nonlinearity in aeroelasticity. Using this method, Hodges and Patil 12,13 investigated the aeroelasticity of the wings of HALE aircraft. They pointed out that the wings of HALE aircraft would be highly deflected in flight and the effect of geometric nonlinearities on the flutter behavior was significant. It is frustrating for the method that the number of independent variables increases exponentially as the number of discrete nodes increases. Hence the geometrically exact beam model initiated by Simo 14 is employed to study the nonlinear structural dynamics. In the numerical computation, the parameterization based on Euler angles is replaced by an implementation based on the use of quaternion parameters, which avoids singularity and minimizes storage requirements. In order to reduce the deviation when the integration time step is large, an improved Newmark algorithm based on implicit Simo-Newmark time stepping integration scheme is proposed.
The remainder of the paper is organized as follows: Firstly, the FSI coupling system, including the coupling conditions, coupling methods of time and space and the fixed-point method with Aitken's dynamic relaxation, is presented. Secondly, the introduction of the geometrically exact beam theory follows, which describes the characteristics of nonlinear dynamics of slender structures. Finally, two numerical examples, which include the flutter analysis of AGARD 445.6 wing to verify the fixed-point method and the flow-induced vibration of a flexible beam to test the feasibility of geometrically exact beam method in FSI problems, are given and discussed in detail.
Descriptions of FSI System. The interaction of the coupling system is essentially the communication of traction forces and motions on the interface of two subsystems 15 . Two partitions are processed by different programs with interaction effects treated as external vector inputs. To clarify the coupling of the problem, the Steklov-Poincare nonlinear equation 16 is introduced to represent the coupling interaction on the interface.
Let F be the Dirichlet-to-Neumann (D-t-N) nonlinear fluid operator, also called Steklov-Poincare operator, mapping fluid displacements u f to interface traction forces λ f as follows, (1) f f It represents the process of the deformation of the fluid domain and computation of flow field. Similarly, a Steklov-Poincare operator for the structural domain, relating structure displacements u s to interface traction forces λ s is express as Concerning the inverse of structure operator, we can define − S 1 as a Neumann-to-Dirichlet (N-t-D) operator associating the interface traction forces and displacements of the structure Generally the operators F and S do not have analytical expressions. The fluid domain is governed by compressible Euler/Navier-Stokes equations in aeroelasticity and the structural domain is solved by linear or nonlinear equations.
Coupling conditions. The interface matching conditions in the following are based on two classical mechanics principles: Eq.(4) represents the kinematic conservation which means the interface displacements of the two subsystems must be equivalent. Eq. (5) is the condition of dynamic conservation meaning that the resultant force on the interface must be zero. The fluid tractions are computed by where P, τ, n f represent the fluid pressure, viscid stress and normal vector on the fluid interface, respectively. Applying Eq. (1) and Eq. (2) to dynamic coupling condition Eq. (5), we can get f s With kinematic coupling condition Eq. (4), Eq. (7) can also be stated based on interface displacements www.nature.com/scientificreports www.nature.com/scientificreports/ Using the definition of the inverse operator S −1 , we can express Eq. (8) as Let G be an operator defined by: Eq. (9) can be reformulated using Eq. (10), resulting with the fixed-point formulation So the standard algorithm to solve problem (11) is based on fixed-point iterations. The fixed-point method is to find interface displacement u Γ defined on Γ. One cycle of FSI iteration is denoted by the operator G.
Time coupling. In FSI problems, the loose coupled or Conventional Serial Staggered method (CSS) 17 is broadly adopted, especially in computational aeroelasticity. At each time step, the D-t-N map F and N-t-D map − S 1 are applyed only once. The procedure is as follows in Algorithm 1.
Data: start and end time(t 0 ,t max ), time step Δt and initial interface displacement Structure solver: end Algorithm 1: CSS method As mentioned in Section 1, this method is only first-order accuracy in time. Generalized Serial Staggered (GSS) method 18 can improve the accuracy by prediction of the structure displacements at the beginning of the computation. The predictor is introduced The prediction can get first-order time accuracy at α 0 = 1 and α 1 =0, and second-order time accuracy at α 0 = 1 and α 1 = 1/2. The GSS algorithm is described as follows in Algorithm 2, Data: start and end time(t 0 , t max ), time step Δt and initial interface displacement Γ u 0 while < t t max do Predict displacements: = + Δ = + t t t n n , 1 end Algorithm 2: Generalized Serial Staggered (GSS) However, loose coupling strategy often leads to the so-called "added mass effect" 7 causing computation instability, which highly depends on the mass ratio between fluid and solid with ρ ρ / s f and the compressibility of the flow 9 . To avoid this drawback, the strong coupling algorithm called the block Gauss-Seidel (BGS) or fixed-point method is proposed. The standard fixed-point iteration with relaxation can be described as: for a given u k , com- and ω k is the relaxation coefficient. When ω k = 1.0, there is no relaxation.
The fixed-point method can be expressed in terms of Algorithm 3. Data: time step Δt, and initial t 0 , final time t ,max , initial interface displacement Γ u 0 while < t t max do Predict interface displacements: The key problem for this algorithm is how to determine the relaxation coefficient. One method called Aitken dynamic relaxation is suggested in Ref. 1 . To illustrate this method, it is of great benefit to consider a scalar fixed-point equation This equation can be solved through the secant method. The iteration formula is: given x k , + x k 1 is calculated by The Aitken relaxation coefficient is defined by Thus Eq. (14) becomes We compute the relaxation coefficient by Analogously, for the fixed-point Eq. (11), the interface residual is defined by Space coupling. Since different solvers are employed, the mesh on the interface of fluid and structural domains usually don't match with each other. So an appropriate transfer scheme of two domains on the interface is important for the accuracy of the whole simulation. The transfer scheme must strictly maintain the conservation of momentum and energy as well as possess high efficiency and accuracy. The equivalence of virtual work is considered to construct a conservative scheme 19 .
Introduce the coupling matrix H mapping the structure displacements to aerodynamic ones, And the relation of virtual displacements of fluid and structure mesh is achieved The structural traction λ s is calculated using Eq. (22) and Eq. (24).
Therefore, once the transfer matrix is determined, it can be used to calculate the aerodynamic nodes displacements and structural traction forces.
In the following, an interpolation scheme with Radial Basis Function (RBF) is presented. RBF is firstly proposed by Wendland and his coworkers 19,20 . The general form of interpolation based on RBFs is www.nature.com/scientificreports www.nature.com/scientificreports/ Coefficients a i and polynomials can be determined by the interpolation condition and the additional condition for all polynomials q with a degree deg(q) ≤ deg(p). The minimal degree of the polynomial depends on the choice of the basis function φ. A unique interpolation that satisfies both Eqs. (27) and (28) is determined if the basis function φ is conditionally positive definite of order m. Only linear polynomials are used here In FSI, the number of fluid and structure nodes N f and Ns on the interface are provided beforehand The points on the structure interface are taken as center points and their displacements are known as u s . Using the interpolation condition Eq. (27) and additional condition Eq. (28), we calculate the coefficients by solving the equations The displacements of fluid nodes on thze interface are interpolated by Combining Eqs. (31) and (32), u f is calculated as follows RBF interpolation can also be used to process mesh deformation. The procedure is similar to the transformation of interface displacements. Nodes on the moving boundaries are firstly selected as center points and the displacements of flow volume grid are interpolated. For large scale of data, a greedy algorithm is adopted to reduce the number of center points.

Geometrically exact Beam theory
To investigate the geometric nonlinearity of slender flexible body, structure nonlinearities solver based on geometrically exact beam theory has been developed by Institute of Mechanics, Chinese Academy of Sciences for Structure dynamics computations. In this section, an introduction of geometrically exact beam theory is presented. For detailed descriptions, refer to the researches 21,22 . Throughout this section, S refers to the arc-length parameter of the rigid beam, and the prime (·)′ represents the derivative with respect to S. The dot ⋅ ⋅ ( ) denotes the derivative with respect to time t. Let For the sake of simplicity, the centerline of the undeformed beam is taken to be a straight line which coincides with the axis of e 3 . Out of all possible configurations, the undeformed state is chosen to be the reference one. In this configuration, choose the triad 3 with d 1 0 and d 2 0 directed along the principal axes of the cross-section and d 3 0 tangent to the curve ϕ 0 . But it is not necessary that the vector field d 3 be tangent to the deformed centroid. Thus such model can describe the shear deformation.
The orthonormal frame (3) represents the special orthogonal group or the rotation group. Compared with standard Cosserat beam theory, elements in SO(3) are regarded as the basic variables and make no further reference to directors field. Hence, a Cosserat beam configuration is completely determined by a pair of curves: . The set Q of all possible configurations of the beam is described by The strain measures are calculated by www.nature.com/scientificreports www.nature.com/scientificreports/ Γ is translational strain measure denoting the shear and extension strain, while K is rotational strain measure denoting the bending and torsional strain.
The constitutive law of the linear elastic material is expressed by the linear relations between the material contact force Λ = n N T and material contact moment Λ = m M T , and the strain measures Γ and K, E and G are the Young's modulus and the shear modulus. A is the cross-section area of the beam, and A 1 , A 2 are referred to the effective areas considering the cross-section distortion. EA represents the elastic axial stiffness of the cross-section, The terms GA 1 ,GA 2 , EI 1 , EI 2 and GJ are regarded as the shear stiffness, the elastic bending stiffness relative to principal axes d 1 , d 2 , and elastic torsional stiffness, respectively.
The weak form of the equilibrium equation is obtained by taking the dot product of Eq. (37) with an arbitrary admissible variation η ν ( , ), which uses direct spatial integration by parts schemes over the domain [0, L]. As a result, it can be written as where G int is the virtual work done by the stress caused by the elastic deformation of material, and given by int L 0 G dyn is the virtual work produced by the inertia of the beam, and computed by dyn L 0 and G ext is the virtual works yielded by the externally applied loads and boundary stress resultants and obtained from the following expression: The procedure of the linearization of Eq. (43) can be found in 14 . And the interpolation of the configuration of beam is discussed in the paper of Zhong et al. 21 .
Improved newmark integration algorithm. In the analysis of nonlinear structural dynamics, implicit schemes are usually taken in the time integration. In terms of computational efficiency, a large enough time step should be adopted. The Newmark integration algorithm is the most widely used method in the solving of structural dynamics. It is described as follows , the Newmark method is second-order accurate in time and unconditionally stable. Note that the construction of the schemes are based on the linear dependence of U n+1 on U n ,  U n , U n . For the centerline of the beam, there is nothing different with the classical scheme. However, the motion of the cross-section of the beam is described by a curve defined in the rotation group SO(3).
SO (3) is not linear space and the rotation matrix Λ + n 1 cannot be estimated by a linear combination of Λ n , Λ  n , Λ̈n. Simo et al. 23 proposed an implicit time stepping algorithm based on the discrete forms of the exponential map and parallel transport in the special orthogonal group SO(3). It performs well in the computation with a small time step but may result in considerable deviations from the exact solution with a large time step. Therefore, a modified Newmark algorithm is presented to improve the accuracy of the scheme given by Simo.
Firstly, Simo-Newmark algorithm is reviewed. For every rotation vector θ, it can be mapped by the exponential map on the rotation matrix Λ ∈ SO (3) www.nature.com/scientificreports www.nature.com/scientificreports/ where θ is the skew matrix associating the rotation vector θ. This is known as Rodrigues' formula. Subscript n represents the temporal discrete approximate of time-varying variables at time t n , u n and θ n denote the incremental displacement and rotation vector from time t n to time t n+1 The basic procedure of the discrete time stepping update is described by: given a configuration Q ( , ) , corresponding accelerations (a n+1 , α n+1 ), in a manner that is consistent and stable with Eq. (43). The implicit time stepping method is represented in Table 1. For convenience, the relations of material and spatial variables are given in Table 2.
It is easy to get for the centerline the update procedure of iterations is described in Table 3. However, Makinen 24 pointed out that there are some flaws in the Simo-Newmark scheme. It needs some modifications to improve the accuracy of time integration. Here we skip the derivation of the correction to the Simo-Newmark scheme and give the results directly. Let Θ Θ = denote the Euclidean norm of the vector Θ and an operator is defined by

Material Spatial Relation
Relative rotation vector Λ Λ Θ = exp( ) 2 1  The modified implicit time stepping algorithm is described in Table 4, where Θ + T( )   Table 3. Update procedure of the implicit time stepping algorithm.

Translation Rotation
Predictor phase Corrector phase for k = 0,1,… www.nature.com/scientificreports www.nature.com/scientificreports/ where U is a conservative variable defined as, are convective fluxes and viscous fluxes, define by , v g , P, τ and µ f represent the fluid density, the flow velocities in Cartesian coordinate system, the velocities of the nodes of the grid, the pressure, the viscous stress tensor and the fluid viscosity. Cp is the specific heat and I is a 3rd-order unit tensor. Specially, for inviscid flows the viscous fluxes disappear and the equations are reduced to Euler equations. Before solving of the equations, the new fluid grid have to be generated due to the moving boundaries. In SU2, the volumetric deformation is based on a classical spring method which may be computational costly and inefficient. So a fast multivariate interpolation method based on RBF with data reduction is applied to do the volumetric mesh deformation. www.nature.com/scientificreports www.nature.com/scientificreports/ where Φ i is the mode shape.
Introducing nondimensional mass ratio µ, wing flutter boundary can be expressed by the reduced velocity ⁎ V as where ∞ u , b and ω α are the free stream velocity, the half-chord length of the root and the first torsion modal frequency, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ The FSI was performed through fixed-point iteration with Aitken's relaxation with the time step size of × − 5 10 4 s. Initial analysis provided the steady-state flow field with the rigid wing in it. The time histories of the generalized displacements in all four modes that developed in response to aerodynamic forces are shown in Fig. 3 with different reduced velocity V * at Ma = 0.96 and the angle of attack of 0 degree. With the different ∞ u , several free stream static pressure and static temperature would be changed at a given Mach number. It can be seen that the all modal displacement amplitudes are in a decreasing trend when ⁎ V is smaller than the flutter boundary. With the value of ⁎ V close to the critical value, the aeroelastic equilibrium response is achieved. For ⁎ V bigger than the flutter boundary, the all modal displacement amplitudes show an increasing trend to the diverging response. At Ma = 0.96, the predicted flutter boundary value is 0.3287, which is close to the experimental result 0.3059 provided by ref. 27 . The flutter speed for Mach 0.499 to 1.072 for the current method and the experimental results are shown in Fig. 4. On the whole, the numerical results at subsonic and transonic regime agree well with the experimental data. However, the flutter boundary values are higher in the computational solutions than the experimental ones in the supersonic regime. This difference between the predicted and measure results has also been noted in the previous studies, which is commonly believed that there exist measurement errors in the experiment of supersonic regime.

Flow-induced vibration of a flexible beam.
To demonstrate the ability of geometrically exact beam to deal with flexible structure in the FSI application, we consider a slender elastic beam attached to a square rigid body immersed in the flow. The geometry and the boundary conditions are shown in 0 35 for Poisson ratio coefficient. The first natural frequency of the problem f 1 = 3.03 Hz.
The fluid domain was solved through unsteady compressible Navier-Stokes equations with dual time method. For the time integration of the structural computation, the improved Newmark method described in Section 2 was adopted and SU2 was also used as the fluid solver.
This application originally carried out in 28 has been used as a benchmark for FSI problem. A hybrid mesh was generated for CFD computation as shown in Fig. 6. To capture the vortex shedding accurately, refined  www.nature.com/scientificreports www.nature.com/scientificreports/ quadrilateral elements were generated in the boundary layers and triangular elements are generated for the rest. The number of total points was 18984 and the number of total cells was 23872. The flexible beam was modeled by geometrically exact beam and discretized into 41 nodes and 40 finite element cells.
The FSI problem was solved with dynamic Aitken's relaxation besides using the second-order predictor to estimate the interface displacements. The Aitken relaxation technique was used with an initial value of 0.5, which rapidly increased to the range of values from 0.8 to 1.0. The convergence tolerance of the interface displacement residual, which is depicted as TOL in Algorithm 3, is set to − 10 7 . The comparisons of the displacement at the free end of the beam and corresponding Fourier analysis at the time steps Δt 1 = 0.005 s and Δt 2 = 0.0075 s are plotted in Fig. 7.
The following observations were made from the Fig. 7. The deviation of the responses obtained for Δt 1 and Δt 2 is very small. At the beginning of the simulation, the flow field was in a symmetric state and the cantilever was stationary. Then the flow exhibited an unsteady behavior with vortices shedding from the corner of the square body to induce oscillations of the flexible beam. After 3 seconds, the whole system reached to an almost harmonic response dominated by the vibration frequency of 3.2 Hz, close to the first eigenfrequency of the structure and the amplitude of the response is about 1.05. The amplitude and the dominant frequency of the response obtained both agree well with the previous results from the literature 26,28,29 .
The flow pressure and velocity distributions at several instants during stable oscillations are displayed in Fig. 8. The deformed shapes of the flexible beam given by Fig. 8 reveal the vibration dominated by the first mode.
The diagrams in Fig. 9 display the evolution of the tip displacement of the beam in time considering different coupling methods. Figure 9(a) shows that the responses obtained by loose coupling and by strong coupling deviate only by a small amount at the time step Δt 1 = 0.005 s. Contrary to Δt 1 = 0.005 s, the simulation from the loose coupling diverges due to added mass effect after hundreds of time steps though the convergence of the tip  www.nature.com/scientificreports www.nature.com/scientificreports/ displacement obtained by the strong coupling is observed be good at Δt 2 = 0.0075 s. The results confirm the fact that the loose coupling method is more sensitive to added mass effect and diverges more easily with the increase of the time step size. Thus the use of the strong coupling method is justified and efficient to solve the coupled problem with the larger time step size.
The final study is considered with respect to different relaxation methods. The number of each FSI iterations is shown in Fig. 10(a). It can be seen that no more than 8 cycles in FSI iterations are required to reach the tolerance TOL = − 10 7 with Aitken's dynamic relaxation, nevertheless, methods with a fixed relaxation parameter ω = . 0 5 and no relaxation ω = .
1 0 lead to more iterations, especially a noticeable increase in fixed relaxation. Figure 10(b) presents the residual convergence for different relaxation methods. With constant ω, a constant decrease of the residual with a low order is observed. With Aitken's relaxation, the convergence exhibits a less smooth trend, with a faster decrease of the residual. Hence, Aitken's dynamic relaxation is of great efficiency to speed up the FSI computation.   www.nature.com/scientificreports www.nature.com/scientificreports/ conclusions In the work, we consider the fluid-structure interaction as the model problem, and discuss its solution procedure from the point view of finding the solution of the interface Steklov-Poincare equation. The computation is performed by the fixed-point method with Aitken's dynamic relaxation. The open source CFD code SU2 is extended with the fixed-point method and integrated with two different structure solvers: the linear modal solver and the nonlinear dynamics solver based on geometrically exact beam theory of Simo. Numerical examples are provided to give evidence of the method effectiveness and efficiency. The case of AGARD 445.6 wing is computed by modal approach and proves the applicability of the fixed-point method in the analysis of transonic flow. The problem of flow-induced vibration of a flexible beam is handled by the fixed-point with Aitken's relaxation in the strong coupling scheme. The strategy exhibits very good convergence properties in FSI iterations, as well as a sufficiently robust performance, especially in a relatively large time step size.