Fluid structure interaction (FSI) problems, which couple fluids to structures, are confronted in many fields such as civil engineering, especially in aeroelasticity. With an increased interest in the use of high-altitude-long-endurance (HALE) aircraft characterized by the long and slender wing structures, large displacements may occur at nominal operating conditions and as a result, new methods need to be developed to study the aeroelastic behaviors of such structures. The numerical simulations in the FSI computation of aeroelastic problems are mainly oriented towards two schemes: monolithic approach and the partition approach. The former solves the aerodynamic forces and structural motions simultaneously in the same manner with no need to distinguish between the fluid and structure part. It is possible to be implemented for simple structural problems with a few of degrees of freedom. The latter solves the aerodynamic forces and structural responses as separate modules in a different manner with fluid-structure interface to pass information between these domains. Because of using the existing methods for both parts independently, the partition approach is a favorite choice in current aeroelastic analysis.

The Dirichlet-Neumann partitioning^{1} using the existing solvers is one of partitioning techniques. In this approach, the nonlinear equations of fluid and structure domains are analysed by different solvers of subsystems. Partition approach usually adopts two different coupling strategies: loose coupling and strong coupling which are both commonly used in previous studies in computational aeroelasticity.

The loose coupling procedure is described as follows: (a) transfer the displacements on the interface boundary of the structure to the fluid system, (b) update the flow mesh accordingly and advance the fluid system to the next time step and compute the fluid traction forces, (c) convert the fluid interface forces into structural loads and advance the structure system to the next time step. Although this approach is easily realized and broadly used in the work^{2,3,4}, it can obtain only first-order accuracy in time. To obtain a higher order accuracy in time, some other improved loose coupling algorithms which can get second-order accuracy in time are given, but the methods don’t enforce the synchronization of the two subsystems.

Besides the disadvantage of low order time accuracy, the loose coupling also shows instability and may introduce spurious solutions in simulation with large time step. Therefore, the strong coupling method is suggested. Yang *et al*.^{5} and Melville *et al*.^{6} described a strongly coupled algorithm in which aerodynamic loads were computed for each pseudo-time step and then structural response produced by those loads were computed, alternately. Though the fluid and structure systems can be perfectly synchronized at each time step with enough iterations, the computation is of high cost and the unsteady flow solver also needs to be modified.

Significant amount of research on aeroelasticity has focused on the structural motions governed by linear dynamic equations. But the analysis procedures based on small deformation assumption fail to analyze characteristics of the structures with large displacement. Block-Gauss-Seidel or fixed-point method^{7,8,9,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 ({mathscr{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,

$${{boldsymbol{lambda }}}_{f}={mathscr{F}}({{boldsymbol{u}}}_{f})$$

(1)

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

$${{boldsymbol{lambda }}}_{s}={mathscr{S}}({{boldsymbol{u}}}_{s})$$

(2)

Concerning the inverse of structure operator, we can define ({{mathscr{S}}}^{-1}) as a Neumann-to-Dirichlet (N-t-D) operator associating the interface traction forces and displacements of the structure

$${{boldsymbol{u}}}_{s}={{mathscr{S}}}^{-1}({{boldsymbol{lambda }}}_{s})$$

(3)

Generally the operators ({mathscr{F}}) and ({mathscr{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:

$${{boldsymbol{u}}}_{f}={{boldsymbol{u}}}_{s}={{boldsymbol{u}}}_{Gamma },{rm{on}},Gamma $$

(4)

$${{boldsymbol{lambda }}}_{f}+{{boldsymbol{lambda }}}_{s}=0,,{rm{on}},Gamma $$

(5)

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

$${{boldsymbol{lambda }}}_{f}=-Pcdot {{boldsymbol{n}}}_{f}+tau cdot {{boldsymbol{n}}}_{f}$$

(6)

where *P*, (tau ), *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

$${mathscr{F}}({{boldsymbol{u}}}_{f})+{mathscr{S}}({{boldsymbol{u}}}_{s})=0$$

(7)

With kinematic coupling condition Eq. (4), Eq. (7) can also be stated based on interface displacements

$${mathscr{F}}({{boldsymbol{u}}}_{Gamma })+{mathscr{S}}({{boldsymbol{u}}}_{Gamma })=0$$

(8)

Using the definition of the inverse operator ({mathscr{S}})^{−1}, we can express Eq. (8) as

$${{mathscr{S}}}^{-1}({mathscr{F}}({{boldsymbol{u}}}_{Gamma }))+{{boldsymbol{u}}}_{Gamma }=0$$

(9)

Let ({mathscr{G}}) be an operator defined by:

$${mathscr{G}}={{mathscr{S}}}^{-1}circ (-{mathscr{F}})$$

(10)

Eq. (9) can be reformulated using Eq. (10), resulting with the fixed-point formulation

$${mathscr{G}}({{boldsymbol{u}}}_{Gamma })={{boldsymbol{u}}}_{Gamma }$$

(11)

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 ({mathscr{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 ({mathscr{F}}) and N-t-D map ({{mathscr{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 ({{boldsymbol{u}}}_{varGamma }^{0})

**while** (t < {t}_{max }) do

Fluid solver: ({{boldsymbol{lambda }}}_{f}^{n+1}={mathscr{F}}({{boldsymbol{u}}}_{Gamma }^{n}))

Structure solver: ({{boldsymbol{u}}}_{varGamma }^{n+1}={{mathscr{S}}}^{-1}(,-,{{boldsymbol{lambda }}}_{f}^{n+1}))

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

$${mathscr{P}}({{boldsymbol{u}}}_{Gamma }^{n})={{boldsymbol{u}}}_{Gamma }^{n}+{alpha }_{0}Delta t{dot{{boldsymbol{u}}}}_{Gamma }^{n}+{alpha }_{1}Delta t({dot{{boldsymbol{u}}}}_{Gamma }^{n}-{dot{{boldsymbol{u}}}}_{Gamma }^{n-1})$$

(12)

The prediction can get first-order time accuracy at ({alpha }_{0}) = 1 and ({alpha }_{1}) =0, and second-order time accuracy at ({alpha }_{0}) = 1 and ({alpha }_{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 ({{boldsymbol{u}}}_{Gamma }^{0})

while (t < {t}_{max }) do

Predict displacements: ({{boldsymbol{u}}}_{Gamma ,p}={mathscr{P}}({{boldsymbol{u}}}_{Gamma }^{n}))

Fluid solver: ({{boldsymbol{lambda }}}_{f}^{n+1}={mathscr{F}}({{boldsymbol{u}}}_{Gamma ,P}))

Structure solver: ({{boldsymbol{u}}}_{Gamma }^{n+1}={{mathscr{S}}}^{-1}(,-,{{boldsymbol{lambda }}}_{f}^{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 ({rho }_{s}/{rho }_{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}, compute ({{boldsymbol{u}}}_{k+1}={{boldsymbol{u}}}_{k}+{omega }_{k}({bar{{boldsymbol{u}}}}_{k}-{{boldsymbol{u}}}_{k})),where ({bar{{boldsymbol{u}}}}_{k}={mathscr{G}}({{boldsymbol{u}}}_{k})) 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 ({{boldsymbol{u}}}_{Gamma }^{0})

**while** (t < {t}_{max }) do

Predict interface displacements: ({{boldsymbol{u}}}_{varGamma ,p}^{n+1}={mathscr{P}}({{boldsymbol{u}}}_{varGamma }^{n}))

*k* = 0

**while** (k < {k}_{max }) do

FSI iteration: ({bar{{boldsymbol{u}}}}_{Gamma ,k}^{n+1}={mathscr{G}}({{boldsymbol{u}}}_{Gamma ,k}^{n+1}))

Compute interface residual: ({r}_{k}^{n+1}={bar{{boldsymbol{u}}}}_{Gamma ,k}^{n+1}-{{boldsymbol{u}}}_{Gamma ,k}^{n+1})

Update the interface displacements ({{boldsymbol{u}}}_{Gamma ,k+1}^{n+1}={{boldsymbol{u}}}_{Gamma ,k}^{n+1}+{{rm{omega }}}_{k}{r}_{k}^{n+1})

if (Vert {r}_{k}^{n+1}Vert < TOL) then

Break

else

*k* = *k* + 1

end

end(t=t+Delta t,n=n+1)end

**Algorithm 3** Fixed-point method with relaxation

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

$$r(x)=f(x)-x=0,xin R$$

(13)

This equation can be solved through the secant method. The iteration formula is: given ({x}_{k}), ({x}_{k+1}) is calculated by

$${x}_{k+1}={x}_{k}-{r}_{k}frac{{x}_{k}-{x}_{k-1}}{{r}_{k}-{r}_{k-1}}$$

(14)

The Aitken relaxation coefficient is defined by

$${omega }_{k}=frac{{x}_{k}-{x}_{k-1}}{{r}_{k}-{r}_{k-1}}$$

(15)

Thus Eq. (14) becomes

$${x}_{k+1}={x}_{k}+{omega }_{k}{r}_{k}$$

(16)

We compute the relaxation coefficient by

$${x}_{k+1}={x}_{k}+{r}_{k},{omega }_{k+1}=-{omega }_{k}frac{{r}_{k}}{{r}_{k+1}-{r}_{k}}$$

(17)

Analogously, for the fixed-point Eq. (11), the interface residual is defined by

$${r}_{k}^{n+1}={mathscr{G}}({{boldsymbol{u}}}_{Gamma ,k}^{n+1})-{{boldsymbol{u}}}_{Gamma ,k}^{n+1}$$

(18)

and the corresponding relaxation coefficient is calculated by

$${omega }_{k+1}=-{omega }_{k}frac{{r}_{k}^{n+1}cdot ({r}_{k+1}^{n+1}-{r}_{k}^{n+1})}{{Vert {r}_{k+1}^{n+1}-{r}_{k}^{n+1}Vert }^{2}}$$

(19)

For the first FSI cycle, it’s suggested to use the ({omega }^{n}) of the previous time step with a constrained parameter^{1}

$${omega }_{0}^{n+1}=,max ({omega }^{n},{omega }_{max })$$

(20)

Therefore, the interface displacements of the next step are

$${{boldsymbol{u}}}_{Gamma ,k+1}^{n+1}={{boldsymbol{u}}}_{Gamma ,k}^{n+1}+{omega }_{k}{{boldsymbol{r}}}_{Gamma ,k}^{n+1}$$

(21)

### 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}.

$$delta W=delta {{boldsymbol{u}}}_{f}cdot {lambda }_{f}=delta {{boldsymbol{u}}}_{s}cdot {lambda }_{s}$$

(22)

Introduce the coupling matrix **H** mapping the structure displacements to aerodynamic ones,

$${{boldsymbol{u}}}_{f}={bf{H}}cdot {{boldsymbol{u}}}_{s}$$

(23)

And the relation of virtual displacements of fluid and structure mesh is achieved

$$delta {{boldsymbol{u}}}_{f}={bf{H}}cdot delta {{boldsymbol{u}}}_{s}$$

(24)

The structural traction *λ*_{s} is calculated using Eq. (22) and Eq. (24).

$${lambda }_{s}={{bf{H}}}^{T}cdot {lambda }_{f}$$

(25)

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

$$s(x)=mathop{sum }limits_{i=1}^{N}{a}_{i}phi (Vert x-{x}_{i}Vert )+p(x)$$

(26)

where (phi (,cdot ,)) is a specific kind of RBF, *x*_{i} is coordinates of center points, *p*(*x*) is a polynomial.

Coefficients *a*_{i} and polynomials can be determined by the interpolation condition

$${d}_{a}=mathop{sum }limits_{i=1}^{N}{a}_{i}phi (Vert {x}_{a}-{x}_{i}Vert )+p({x}_{a})$$

(27)

and the additional condition

$$mathop{sum }limits_{i=1}^{N}{a}_{i}q(Vert {x}_{i}Vert )=0$$

(28)

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 (phi ). A unique interpolation that satisfies both Eqs. (27) and (28) is determined if the basis function (phi ) is conditionally positive definite of order *m*. Only linear polynomials are used here

$$p(x)={beta }_{0}+{beta }_{1}x+{beta }_{2}y+{beta }_{3}z$$

(29)

In FSI, the number of fluid and structure nodes *N*_{f} and *Ns* on the interface are provided beforehand

$$begin{array}{c}X={begin{array}{cccc}{x}_{1}, & {x}_{2,} & cdots & {x}_{{N}_{s}}end{array}},\ Y={begin{array}{cccc}{y}_{1}, & {y}_{2,} & cdots & {y}_{{N}_{f}}end{array}}end{array}$$

(30)

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

$$[begin{array}{c}{{boldsymbol{u}}}_{s}\ 0end{array}]=[begin{array}{cc}{{bf{C}}}_{ss} & {{bf{P}}}_{s}\ {{bf{P}}}_{s}^{T} & 0end{array}],[begin{array}{c}alpha \ beta end{array}]$$

(31)

where ({{bf{C}}}_{ss}=[begin{array}{cccc}{phi }_{{s}_{1}{s}_{1}} & {phi }_{{s}_{1}{s}_{2}} & cdots & {phi }_{{s}_{1}{s}_{{N}_{s}}}\ {phi }_{{s}_{2}{s}_{1}} & {phi }_{{s}_{2}{s}_{2}} & cdots & {phi }_{{s}_{2}{s}_{{N}_{s}}}\ vdots & vdots & ddots & vdots \ {phi }_{{s}_{{N}_{s}}{s}_{1}} & {phi }_{{s}_{{N}_{s}}{s}_{2}} & cdots & {phi }_{{s}_{{N}_{s}}{s}_{{N}_{s}}}end{array}]), ({{bf{P}}}_{s}=[begin{array}{cccc}1 & {x}_{{s}_{1}} & {y}_{{s}_{1}} & {z}_{{s}_{1}}\ 1 & {x}_{{s}_{2}} & {y}_{s2} & {z}_{{s}_{2}}\ vdots & vdots & vdots & vdots \ {1}_{{s}_{{N}_{s}}{s}_{1}} & {x}_{{s}_{{N}_{s}}} & {y}_{{s}_{{N}_{s}}} & {z}_{{s}_{{N}_{s}}}end{array}]), ({phi }_{{s}_{i}{s}_{j}}=phi (Vert {x}_{{s}_{i}}-{x}_{{s}_{j}}Vert ))

The displacements of fluid nodes on thze interface are interpolated by

$${{boldsymbol{u}}}_{f}=[begin{array}{cc}{{bf{M}}}_{fs} & {{bf{P}}}_{f}end{array}][begin{array}{c}alpha \ beta end{array}]$$

(32)

where ({{bf{M}}}_{fs}=[begin{array}{cccc}{phi }_{{f}_{1}{s}_{1}} & {phi }_{{f}_{1}{s}_{2}} & cdots & {phi }_{{f}_{1}{s}_{{N}_{s}}}\ {phi }_{{f}_{2}{s}_{1}} & {phi }_{{f}_{2}{s}_{2}} & cdots & {phi }_{{f}_{2}{s}_{{N}_{s}}}\ vdots & vdots & ddots & vdots \ {phi }_{{f}_{{N}_{s}}{s}_{1}} & {phi }_{{f}_{{N}_{s}}{s}_{2}} & cdots & {phi }_{{f}_{{N}_{s}}{s}_{{N}_{s}}}end{array}]), ({{bf{P}}}_{f}=[begin{array}{cccc}1 & {x}_{{f}_{1}} & {y}_{{f}_{1}} & {z}_{{f}_{1}}\ 1 & {x}_{{f}_{2}} & {y}_{{f}_{2}} & {z}_{{f}_{2}}\ vdots & vdots & vdots & vdots \ 1 & {x}_{{f}_{{N}_{f}}} & {y}_{{f}_{{N}_{f}}} & {z}_{{f}_{{N}_{f}}}end{array}]), ({phi }_{{f}_{i}{s}_{j}}=phi (Vert {x}_{{f}_{i}}-{x}_{{s}_{j}}Vert ))

Combining Eqs. (31) and (32), *u*_{f} is calculated as follows

$${{boldsymbol{u}}}_{f}=[begin{array}{cc}{{bf{M}}}_{fs} & {{bf{P}}}_{f}end{array}]{[begin{array}{cc}{{bf{C}}}_{ss} & {{bf{P}}}_{s}\ {{{bf{P}}}_{s}}^{T} & 0end{array}]}^{-1}[begin{array}{c}{{boldsymbol{u}}}_{s}\ 0end{array}]=hat{{bf{H}}}[begin{array}{c}{{boldsymbol{u}}}_{s}\ 0end{array}]$$

(33)

Thus the matrix **H** is built by the first *N*_{f} rows and *Ns* columns of (hat{{bf{H}}})

$${bf{H}}=hat{{bf{H}}}[begin{array}{cc}1:{N}_{f}, & 1:{N}_{s}end{array}]$$

(34)

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.