Analysis of some partitioned algorithms for fluid-structure interaction

Published in ''Engineering Computations'' Vol. 27 (1), pp. 20-56, 2010<br />
doi: 10.1108/02644401011008513
==Abstract==

The field of fluid-structure interaction (FSI) has important applications ranging from the bio-medical to the civil and aeronautical engineering fields. Despite the different approaches developed over recent years the general solution of FSI problems remains however far from being closed. This paper addresses the theoretical analysis of a several partitioned solution schemes  for FSI. This is achieved by choosing a simple but representative model problem which allows to obtain analytical results. A loose coupling partitioned procedure is presented first. An iterative version of the same algorithm is presented next together with an analysis of its convergence properties and a brief discussion of a possible acceleration strategy. Finally a new iterative strategy based on solving a modified equation  is described. Some examples of application of the partitioned algorithms considered are presented including the aeroelastic analysis of a bridge.

'''Keywords''': Fluids, Structures, Flow, Fluid dynamics, Programming and algorithm theory, Iterative methods

==1 Introduction==

The coupling between different physical phenomena constitute a  challenge in many applications of practical interest. Empirical or analytical techniques which classically provide a satisfactory description of coupled problems at the asymptotic limits, do not guarantee the necessary guidance when the coupling becomes strong. The field of FSI (fluid-structure interaction) represents an important and active area of research due  its relevance in many problems in science and engineering.

Until recent times the numerical simulation of strongly coupled FSI problems was considered to be unviable. Experimental wind-tunnel testing were for long  the only approach for the design of complex aeroelastic interaction-sensitive structures such as bridge decks or slender towers. Although well established, wind tunnel testing is often expensive and time-consuming. In addition some concerns on the quality of the predictions arise when severe scale reduction is needed to fit the model in the test facility.

The maturity of numerical finite element (FE) techniques in the fields of computational fluid dynamics (CFD) and structural dynamics, together with the increasing computational power start nowadays to make competitive the purely numerical approach for solving FSI problems.

From a mathematical point a wide class of coupled problems can be represented symbolically as

<span id="eq-1"></span>
<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_S(x,y) = 0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|-
| style="text-align: center;" | <math>     f_F(x,y) = 0       </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}
|}

where <math display="inline">f_S(x,y)</math> and <math display="inline">f_F(x,y)</math> are functions that denote the behaviour of the two coupled systems and <math display="inline">x,y</math> represent some significative variables for each of the systems. Typically subscripts <math display="inline">S</math> and <math display="inline">F</math> denote structure and fluid contributions, respectively. Eqns.([[#eq-1|1]]) and ([[#eq-2|2]])  can be linearized  to give a system in the form

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{pmatrix}\frac{\partial f_S}{\partial x}  &  \frac{\partial f_S}{\partial y}  \\          \frac{\partial f_F}{\partial y}  &  \frac{\partial f_F}{\partial y}      \end{pmatrix}     \left(\begin{array}{c}dx \\ dy \end{array} \right)= \left(     \begin{array}{c}-f_S(x^k,y^k) \\ -f_F(x^k,y^k) \end{array}     \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where the off-diagonal terms of the tangent (or Jacobian) stiffness matrix express the dependency of one field on the variations of the second. This form guarantees theoretically optimal quadratic convergence properties. In practice  the two symbolic operators <math display="inline">f_S(x,y)</math> and <math display="inline">f_F(x,y)</math> are non  linear and in most cases are not explicitly known. Hence  the computation of the derivatives in Eqn([[#eq-3|3]]) may be  difficult or undesiderable for other reasons.

Even if theoretically possible, the linearization of the fluid field is usually not performed, as it is often preferred to rely on fixed point type iterations to solve the non linearities. This implies that the computation of the jacobian in Eqn([[#eq-3|3]]) would involve a redefinition of the fluid strategy which is not desiderable from the point of view of software modularity. This makes full Newton schemes generally unattractive for solving FSI problems. Further the resulting “monolithic” matrix, would inherit the same conditioning problems from the corresponding fluid counterpart. Thus making its inversion difficult using any iterative technique.

Two main approaches exist for solving this impasse. The purely algebraic one focuses on the iterative solution of Eqns([[#eq-1|1]]) and ([[#eq-2|2]]) approaching as much as possible the optimal behavior of the Newton system in Eqn([[#eq-3|3]]). The second one focuses on partitioned (sometimes called “explicit”) coupling techniques where convergence to the coupled solution is guaranteed by a careful adjustment of the predicted variable and of the data transferred between the interacting domains.

The first approach has its simpler expression in the so called Block-Jacobi or Block-Seidel iterations which  involve the following steps

* <math display="inline"> f_S(x^{k+1},y^k) = 0 \rightarrow y^{k+1}</math>
* <math display="inline"> f_F(x_\omega ^{k+1},y^{k+1})= 0 \rightarrow x^{k+1} </math>

The above scheme can be eventually modified to include a relaxation step such as <math display="inline">x_w^{k+1} = \omega x^{k+1} + (1-\omega ) x^{k}</math>. This technique is very effective for some problems. Nevertheless in some cases  may not be computationally competitive versus the partitioned technique. The choice of the optimal relaxation factor is still matter of research. For instance in <span id='citeF-1'></span>[[#cite-1|[1]]],  <span id='citeF-2'></span>[[#cite-2|[2]]]  the  Aitken accelerator is used to improve the convergence of the scheme. Alternatively a number of quasi-Newton techniques were presented over the years in an effort for achieving quadratic convergence when solving Eqn([[#eq-3|3]]). An interesting description can be found in <span id='citeF-3'></span>[[#cite-3|[3]]], <span id='citeF-4'></span>[[#cite-4|[4]]] where iterative solvers are exploited for solving the coupled problem without building explicitly the Jacobian matrix.

This matrix-free approach is optimal from the point of view of software modularity, and naturally tends to the solution of the monolithic formulation from which it inherits stability properties. The  disadvantage is the high computational cost due to the need for solving the implicit system to a high degree of accuracy.

In addition, the partitioned approach relies on a careful design of the prediction and correction phases in order to minimize the spurious energy introduced in the system. Pioneering work on the subject was reported by <span id='citeF-5'></span>[[#cite-5|[5]]] who developed partitioned solution strategies for large scale parallel computing applications in aeronautics<span id='citeF-6'></span> [[#cite-6|[6]]].

Given the absence of a common mathematical formulation for the two coupled fields the investigation of the stability properties of the different methods is extremely challenging. A 1D case is investigated in <span id='citeF-6'></span>[[#cite-6|[6]]] (for a compressible flow) while in <span id='citeF-5'></span>[[#cite-5|[5]]] a simple test is advocated to assess the stability of a partitioned scheme. Even if the proposed method  has no general mathematical validity it was shown to discriminate successfully stable methods from unstable ones. The techniques proposed were tested on large aeronautical examples <span id='citeF-7'></span>[[#cite-7|[7]]] and provided satisfactory results in describing experimental flutter envelopes on real aircrafts. They are also believed to be effective in dealing with linear and non-linear problems as long as the displacements are small. Some concern arises however regarding their robustness in dealing with FSI problems for highly flexible structures and a truly incompressible fluid. Recent results however have proven their applicability to general cases when using compressible flow solvers.

The application of space-time FE based approaches can be found in <span id='citeF-8'></span>[[#cite-8|[8]]] and <span id='citeF-9'></span>[[#cite-9|[9]]] with reference to both civil and aeronautical applications. Finally some examples of applications to large flexible civil engineering structures can be found in <span id='citeF-10'></span>[[#cite-10|[10]]],<span id='citeF-11'></span>[[#cite-11|[11]]].

The objective of this work is to discuss several algorithms for solving coupled problems of the type of Eqn([[#eq-2|2]]) by considering  a loose coupling  partitioned approach and two possible improvements by iteration. The effort is to provide a framework to study the properties of the coupling procedures  by evaluating the impact of different time integration schemes. The layout of the paper is the following. First a simple but representative 1D model problem is chosen. A partitioned solution strategy is identified and applied symbolically to the model problem. This allows us to obtaining  an error matrix which depends on the  time integration scheme chosen. Two possibilities of iterative improvement of the partitioned method is discussed and an estimate of the error evolution with the iterations is given. Finally an enhanced iterative strategy based on the solution of a modified equation is described. The partitioned schemes studied are applied to a number of FSI problems including to the solution of a real problem in bridge aerodynamics for which wind tunnel results are available for comparisons.

==2 The Model Problem==

Partitioned (or staggered)  techniques are appealing for solving coupled FSI problems as they allow the use of state-of-the-art single-field solvers. As such they are widely used in many fields of application. The disadvantage is often connected to a lack of robustness in dealing with some category of FSI problems, typically involving truly incompressible fluids or requiring long-term analyses. Strong or fully coupling techniques on the other hand guarantee a convergence to the monolithic solution by providing a completely “implicit” coupling between the different fields. In practice however, they tend to be costly and may experience convergence difficulties when dealing with complex problems.

The aim of the paper is to investigate the reasons behind this difficulties. The same effort was attempted in different works, see for example <span id='citeF-5'></span>[[#cite-5|[5]]] for a discussion of partitioned methods or <span id='citeF-12'></span>[[#cite-12|[12]]] for an analysis of a semi-implicit technique. Both papers present interesting results and define satisfactory algorithms for many applications.

This work follows a slightly different path by considering a further requirement: an ideal coupling scheme should be robust upon changes in the single-field solvers. In other words, its properties should be independent from the particular solver chosen for each of the coupling fields. As a consequence, the particular features of the solvers  should not be exploited in discussing a coupling algorithm, unless the interest focuses exclusively on a very particular choice. In order to highlight our approach a simple linearized model problem is chosen and the coupling algorithm proposed is applied to its solution.

The choice for such a simple model problem is unfortunately non univocal in FSI, and needs to be justified. To do so we assume  that in the “common” engineering practice the fluid solution is often avoided and the interaction between a structure and the surrounding fluid is often described by the exchange of a pressure force between the fluid and the structure <span id='citeF-13'></span>[[#cite-13|[13]]]. Such pressure force is often modeled as

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p = A(\omega )  \dot{x} + B(\omega )   {x}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where  coefficients <math display="inline">A</math> and <math display="inline">B</math> generally depend on the  frequency <math display="inline">\omega </math> of the motion. A  careful choice of such constants for aeroelastic analysis (generally from a good fit  to wind tunnel data) can lead to an excellent description  of the  most relevant features of the wind action on a structure. This implies  that this simple model can reproduce correctly  the coupling between the structure and the fluid.

A load defined by Eqn([[#eq-4|4]]), which is naturally suited for the frequency domain analysis, generally needs to be transformed in the time domain. Such transformation is not trivial (see for example <span id='citeF-14'></span>[[#cite-14|[14]]]) and will not be discussed here. The interesting result is that Eqn([[#eq-4|4]]) can be substituted in the time domain by a relation of the type

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p = K_{f} x +  C_{f}  \dot{x} +  M_{f}  \ddot{x}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

eventually generalized including history dependent terms. In Eqn([[#eq-5|5]]) index <math display="inline">f</math> denotes terms in the fluid domain computed at the fluid-structure interface. This implies that the instantaneous force exerted by the fluid on the structure, depends directly on the acceleration, velocity and displacement of the structure itself. The dependence on the acceleration in particular is crucial and the corresponding “mass-like” term is known  as “added mass”.

Interestingly enough a similar dependence of the interaction forces on the structural motion can be obtained  from the analysis of the continuum problem <span id='citeF-12'></span>[[#cite-12|[12]]] or from a simple algebraic manipulation of the discrete Navier Stokes problem. The important point is that a simple model as that of  Eqn([[#eq-5|5]]) captures satisfactorily the important features of the coupling, and is thus appropriate to define a test problem.

The idea we will exploit here is conceptually simple: let us consider a 1D problem representing the interaction between a given structure and the surrounding fluid. Without loss of generality we can assume that the abstract problem defined by Eqns([[#eq-1|1]]) and ([[#eq-2|2]]) can be rewritten as:

<span id="eq-6"></span>
<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_S(x,p) := M_s  \ddot{x} +  C_s  \dot{x} + K_s  x - p=0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|-
| style="text-align: center;" | <math>          f_F(x,p) := p - K_{f} x + C_{f}  \dot{x} +  M_{f}  \ddot{x} =0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}
|}

where indices <math display="inline">s</math> and <math display="inline">f</math> denote terms emanating from the solid and fluid domains ''computed at the fluid-structure interface'', <math display="inline">p</math> is the pressure force that originates from the coupling and it plays the role of the “coupling variable” in the study of the interaction. Eqns([[#eq-6|6]]) and ([[#eq-7|7]]) represent the simplest abstract form of the coupled equations for a FSI problem. The performance of a given coupling algorithm can be therefore assessed by applying it  to the solution of Eqns([[#eq-6|6]]) and ([[#eq-7|7]]).

It is  intuitive that procedures that “behave well” on a test problem with the form of Eqn([[#eq-6|6]]) and ([[#eq-7|7]]) for any physically relevant choice of the problem  parameters can describe correctly the effect of the fluid on the structure. In other words procedures that solve accurately the test problem can be expected to be well suited for the simulation of the FSI interaction for  cases where a model giving the pressure force by an expression of the form of Eqns([[#eq-4|4]]) is acceptable for design purposes.

Clearly the time discretization of Eqns([[#eq-6|6]]) and ([[#eq-7|7]]) implicitly contains an error connected to the  choice for the time integration algorithm. This error has little interest as it does not depend on the way we enforce the coupling. Our interest is rather connected to the error between the approximate and exact discrete solutions of the coupling problem. In other words we assume that no error cancelation manifests due to the interaction of the different solution schemes. Note that this is a strong assumption. Different approaches, in particular the work of Farhat and Piperno <span id='citeF-5'></span>[[#cite-5|[5]]] , exploit exactly the error cancelation to tune the features of their coupled solvers. Not considering the error cancelation is however mandatory if the requirement of “independence from the single-field solvers” needs to be enforced.

As a final comment we point out that necessary features for a “good” coupled time integrator are (Felippa <span id='citeF-15'></span>[[#cite-15|[15]]])

<ol>

<li>preserve the stability of stable mathematical models     </li>
<li>manifest the instability of unstable mathematical models. </li>

</ol>

Effort will be devoted in this work to studying the stability properties of the coupling algorithms analysed.

===2.1 A simple form for the analysis of the discrete coupling problem===

The system described by Eqn([[#eq-6|6]]) and ([[#eq-7|7]]) can be rewritten as a reduced monolithic problem  simply by substituting the second equation into the first. To simplify the subsequent developments we will express the mass and stiffness coefficient for the fluid (<math display="inline"> K_{f} </math>  and <math display="inline"> M_{f}</math> ) as a ratio to the correspondent structural quantities. This gives

<span id="eq-8"></span>
<span id="eq-9"></span>
<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K_s = \omega ^2 M_s\, ; \,  K_{f} = \beta   K_s; </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|-
| style="text-align: center;" | <math>       C_s = 2\xi \omega \, ; \,  C_{f} = 2\xi _{f}\omega </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|-
| style="text-align: center;" | <math>            M_s = 1 \, ; \,  M_{f} = \alpha   M_s </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}
|}

where the structural mass (<math display="inline">M_s</math>) is assumed to have a unit value without loss of generality.

Substituting the value of <math>p</math> in Eqn([[#eq-7|7]]) into Eqn([[#eq-6|6]]) and using Eqns([[#eq-8|8]])&#8211;([[#eq-10|10]]) yields  the monolithic system

<span id="eq-11.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(M_s - M_f) \ddot{x} + (C_s -C_f)  \dot{x} + (K_s-K_f) {x} = {\boldsymbol 0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11.a)
|}

or

<span id="eq-11.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(1-\alpha \right) \ddot{x} + 2\omega \left(\xi  - \xi _{f} \right) \dot{x} + \left(1-\beta \right)\omega ^2  x = {\boldsymbol 0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11.b)
|}

Eqn([[#eq-11.b|11.b]]) or ([[#eq-11.a|11.a]]) are equivalent to Eqns([[#eq-6|6]])  and ([[#eq-7|7]]) and express synthetically the dynamic behaviour of the coupled problem. We are mostly interested in the case where the presence of the fluid influence the transient behaviour of the structure. We shall however assume that even in presence of a very important coupling the mathematical behaviour of the coupled system is preserved, or, in other words that the overall mass and stiffness are both non negative. Note however that  a negative damping <math display="inline"> C_{f} </math> (leading to a divergent structural behaviour) can be considered. These considerations allows us to assume that

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha < 1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta < 1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

In general we can expect the coupled system to be heavier than the single system, which suggests that the stricter condition <math display="inline">\alpha{<0}</math> will hold.

Up to this point no approximation was introduced and Eqn([[#eq-11.a|11.a]]) or ([[#eq-11.b|11.b]]) is still completely equivalent to the system described in Eqns([[#eq-6|6]]) and ([[#eq-7|7]]). To proceed further we perform the discretization in time. Without loss of generality we shall consider that a  general class of time integrators can be expressed symbolically as

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{y}_{n+1} = \mathbf{A} \mathbf{y}_n + \mathbf{L}_n     p_n + \mathbf{L}_{n+1}p_{n+1}\, ;</math> 
|-
{| style="text-align: left; margin:auto;width: 100%;" 
| style="text-align: center;" | <math>\mathbf{y}_n = \left(     \begin{array}{c}x_n \\  \dot{x} _n  \end{array}\right)\, ; \, \mathbf{y}_{n+1} = \left(     \begin{array}{c}x_{n+1} \\  \dot{x} _{n+1}  \end{array}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where the  choice of the linear operators <math display="inline">\mathbf{A}</math>, <math display="inline">\mathbf{L}_n</math> and <math display="inline">\mathbf{L}_{n+1}</math> defines the time integrator used and the effect of external forces. Such general form allows us to predict the evolution of a system between the time stations <math>n</math> and <math>n+1</math> when subjected to a force varying linearly between the values <math display="inline">p_n</math> and <math display="inline">p_{n+1}</math>.

By defining the auxiliary vectors

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{G} = \begin{pmatrix}K & C \end{pmatrix}, \, \, \, \, \mathbf{G}_{f} = \begin{pmatrix}K_{f} & C_{f} \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

we express the model problem of Eqn([[#eq-11.b|11.b]]) in the simple matrix form

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\ddot{x} + \mathbf{G} \mathbf{y} = \mathbf{G}_{f} \mathbf{y} +     M_{f}  \ddot{x}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

and the “pressure force” (the term exchanged between the first and second equation in Eqns([[#eq-6|6]]) and ([[#eq-7|7]])) at the initial and final time stations as 

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{n+1}=  M_{f}\ddot x_{n+1}+  \mathbf{G}_{f}\mathbf{y}_{n+1}, \,     \, \, p_n =  M_{f} \ddot x_{n} + \mathbf{G}_{f}\mathbf{y}_{n} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

This  suffices to express in matrix form both the exact solution and the approximate one.

===2.2 Discrete solution of the exact coupling problem===

The evaluation of the exact coupled solution can be performed preserving the matrix notation introduced in the previous section. To do so it is necessary to express the link between the pressure and the structural motion. Dynamic equilibrium (Eqn([[#eq-6|6]])) gives

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\ddot x_{n+1}=  p_{n+1}- \mathbf{G}\mathbf{y}_{n+1}      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

substituting back into Eqn([[#eq-17|17]]) and using the definitions in Eqns([[#eq-8|8]])&#8211;([[#eq-10|10]]), we recover the relation

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(1-\alpha \right) p_{n+1}= \left(\mathbf{G}_{f}     - \alpha \mathbf{G}\right)\mathbf{y}_{n+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

On the other hand the structure is advanced in time following the time integration rule of Eqn([[#eq-14|14]]). Substituting Eqn([[#eq-19|19]]) into Eqn([[#eq-14|14]]) and collecting terms we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{pmatrix}\mathbf{I} & \mathbf{-L_{n+1}} \\         \alpha \mathbf{G} - \mathbf{G}_{f} & 1-\alpha      \end{pmatrix}     \begin{pmatrix}\mathbf{y}_{n+1} \\  p_{n+1} \end{pmatrix}     =     \begin{pmatrix}\mathbf{A} & \mathbf{L}_{n} \\         0 & 0     \end{pmatrix}     \begin{pmatrix}\mathbf{y}_{n} \\  p_n  \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

which allows us to write synthetically

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1} = \mathbf{A}_{ex}\mathbf{z}_n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{A}_{ex}:=     \begin{pmatrix}\mathbf{I} & \mathbf{-L}_{n+1} \\         \alpha \mathbf{G} -  \mathbf{G}_{f} & 1-\alpha      \end{pmatrix}^{-1}     \begin{pmatrix}\mathbf{A} & \mathbf{L}_{n} \\         0 & 0     \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1} := \begin{pmatrix}\mathbf{y}_{n+1} \\  p_{n+1}\end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

Eqn([[#eq-21|21]]) yields the amplification form of the exact solution for the coupled problem.

===2.3 A fractional step-like algorithm===

In  the previous section we expressed the exact solution in an amplification form, relating pressures, velocities and displacements at two consecutive time steps.

Unfortunately in real life not all the information needed is available at the same time. Therefore an approximate coupling strategy needs  to be devised. In this section we propose and analyse a fractional step-like partitioned approach including the following steps:

''Predictive step'':
* Compute the fluid forces acting on the structure,
* Advance in time the structure under to the predicted forces.
<br />

''Correction step'':
* Compute the new fluid force,
* Solve for the fluid domain according to the structural prediction,
* Update the  position of the structure according to the new fluid forces.
<br />

This simple strategy can be  improved by successive corrections leading to the definition of an iterative strategy which minimizes the coupling error inside each time step. For such a strategy our interest focuses on the convergence properties of the iteration sequence.

In the following lines we detail the prediction and correction phases.

''Prediction step.''  The pressure acting on the structure is calculated as

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{n+1}^I = \eta  p_n + \phi  p_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

Then we ''advance in time'' the solution for the structure under the predicted pressure. Taking into account the form assumed for the structural integrator we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{y}_{n+1}^I = \mathbf{A} \mathbf{y}_n + \mathbf{L}_n p_n + \mathbf{L}_{n+1}p_{n+1}^I </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

grouping the two equations we get

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{pmatrix}\mathbf{I} & \mathbf{-L_{n+1}} \\ 0 & 1   \end{pmatrix}     \begin{pmatrix}\mathbf{y}_{n+1}^I \\  p_{n+1}^I  \end{pmatrix}     =     \begin{pmatrix}\mathbf{A} & \mathbf{L}_{n} \\ 0 & \eta \end{pmatrix}     \begin{pmatrix}\mathbf{y}_{n} \\  p_n  \end{pmatrix}     +     \begin{pmatrix}0 & 0 \\        0 & \phi    \end{pmatrix}     \begin{pmatrix}\mathbf{y}_{n-1} \\ p_{n-1}    \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

By introducing the auxiliary matrices

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{A}^I_1 := \mathbf{C}\begin{pmatrix}\mathbf{A} & \mathbf{L}_{n} \\ 0 & \eta      \end{pmatrix} \; \mathbf{A}^I_2 := \mathbf{C}\begin{pmatrix}\mathbf{0} & \mathbf{0} \\ 0     & \phi \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{C} :=     \begin{pmatrix}\mathbf{I} & \mathbf{-L_{n+1}} \\ 0 & 1   \end{pmatrix}^{-1}     =     \begin{pmatrix}\mathbf{I} & \mathbf{L}_{n+1} \\    0 & 1   \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

we can express the prediction for the fluid and structural variables as

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}^I = \mathbf{A}^I_1 \mathbf{z}_{n} + \mathbf{A}^I_2     \mathbf{z}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

''Correction step.''  The first step of the correction step is the computation of the corrected pressure field. To include this in our simplified model we need to calculate the predicted accelerations. Dynamic equilibrium (Eqn([[#eq-18|18]]) gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\ddot x}_{n+1}^I = p_{n+1}^I - \mathbf{G} \mathbf{y}_{n+1}^I </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Using the definitions in Eqn([[#eq-8|8]]) the corrected  pressure value in the fluid can be expressed as

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{p}_{n+1}^{II} = \mathbf{M}_{f_{s_i}} {\ddot x}_{n+1}^I + \mathbf{G}_{f_{s_i}}\mathbf{y}_{n+1}^I = \alpha p_{n+1}^I + \left(\mathbf{G}_{f_{s_i}} - \alpha \mathbf{G}\right)\mathbf{y}_{n+1}^I      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

The final correction step for the structure has the usual form

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{y}_{n+1}^{II} = \mathbf{A} \mathbf{y}_n + \mathbf{L}_n p_n + \mathbf{L}_{n+1} p_{n+1}^{II}      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

Eqns([[#eq-31|31]] and  [[#eq-32|32]]) can be grouped in an amplification form as

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}^{II} = \mathbf{B}_0 \mathbf{z}_{n+1}^I + \mathbf{B}_1 \mathbf{z}_{n} + \mathbf{B}_2 \mathbf{z}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{B}_0 := \mathbf{C} \begin{pmatrix}0 & 0 \\ \mathbf{G}_{fsi} - \alpha \mathbf{G} & \alpha    \end{pmatrix} =\begin{pmatrix}\mathbf{I} & \mathbf{L}_{n+1} \\    0 & 1   \end{pmatrix}     \begin{pmatrix}0 & 0 \\ \mathbf{G}_{fsi} - \alpha \mathbf{G} & \alpha    \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{B}_1 := \mathbf{C} \begin{pmatrix}\mathbf{A} & \mathbf{L}_n \\ 0 & 0     \end{pmatrix} =\begin{pmatrix}\mathbf{I} & \mathbf{L}_{n+1} \\    0 &     1   \end{pmatrix} \begin{pmatrix}\mathbf{A} & \mathbf{L}_n \\ 0 &     0  \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{B}_2 := \mathbf{C} \begin{pmatrix}\mathbf{0} & \mathbf{0} \\ 0 & 0     \end{pmatrix} =   \begin{pmatrix}\mathbf{0} & \mathbf{0} \\ 0 & 0     \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Note that matrix <math display="inline">\mathbf{B}_2</math> is empty and was included only for generalization purposes. Substituting Eqn([[#eq-29|29]]) into Eqn([[#eq-33|33]]) we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}^{II} = \mathbf{B}_0 \left(\mathbf{A}^I_1 \mathbf{z}_{n} + \mathbf{A}^I_2 \mathbf{z}_{n-1}\right)+ \mathbf{B}_1 \mathbf{z}_{n} + \mathbf{B}_2 \mathbf{z}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

which allows us to obtain the final amplification form for the fractional step procedure as

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}^{II} = \left(\mathbf{B}_0     \mathbf{A}^I_1 + \mathbf{B}_1 \right)\mathbf{z}_{n} + \left(\mathbf{B}_0 \mathbf{A}^I_2 + \mathbf{B}_2 \right)\mathbf{z}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

This form will be used  to study the accuracy and stability  of the proposed scheme in the next sections.

===2.4 Analysis of the truncation error===

Assuming that the exact solution is known in  steps 1 to <math>n</math>, the truncation error of the fractional step algorithm explained in the previous section can be obtained by difference between the approximate and exact solutions.

The amplification form obtained in the previous section Eqn([[#eq-38|38]]) allows us the evaluation of this error. To perform this operation we observe that <math display="inline">\mathbf{y}_n</math> and <math display="inline">\mathbf{p}_n</math> can not be prescribed independently as the exact solution is used as a starting point. Typically only the displacements and velocities are prescribed and the pressure is calculated subsequently. A dynamic equilibrium relation necessarily holds relating pressures displacements and velocities in the initial steps.  The definition of <math display="inline">\mathbf{z}_n</math> implies

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_n = \mathbf{D} \mathbf{y}_n \quad \hbox{with}\quad  \mathbf{D} := \begin{pmatrix}\mathbf{I} \\ \frac{1}{1-\alpha }\left(\mathbf{G}_{f} - \alpha \mathbf{G} \right) \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

The definition of the exact amplification matrix Eqn([[#eq-22|22]]) on the other hand gives

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_n^{ex} = \mathbf{A}_{ex} \mathbf{z}_{n-1}^{ex} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

This together with Eqn([[#eq-39|39]]) allows us to express the exact solution at step <math display="inline">n+1</math> in the form

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}^{ex}_{n+1} = \mathbf{A}_{ex}\mathbf{A}_{ex} \mathbf{D} \mathbf{y}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

Following the same rationale, the approximate solution at the given step becomes

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}^{ap}_{n+1} = \left(\left(\mathbf{B}_0 \mathbf{A}^I_1     + \mathbf{B}_1 \right)\mathbf{A}_{ex} + \left(\mathbf{B}_0 \mathbf{A}^I_2 + \mathbf{B}_2     \right)\right)\mathbf{D} \mathbf{y}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

By subtracting Eqns([[#eq-41|41]]) and ([[#eq-42|42]]) it is possible to define a “reduced” rectangular error matrix relating the pressure and the structural variables at  step <math display="inline">n+1</math> with the structural variables at the former step, i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}^{ex} - \mathbf{z}_{n+1}^{ap} = \mathbf{E}^{r} \mathbf{y}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

with

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}^{r} := \left(\mathbf{A}_{ex}\mathbf{A}_{ex} - \left(\mathbf{B}_0 \mathbf{A}^I_1 + \mathbf{B}_1 \right)\mathbf{A}_{ex} - \left(\mathbf{B}_0 \mathbf{A}^I_2 + \mathbf{B}_2 \right)\right)\mathbf{D} \mathbf{y}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

This last equation contains all the information needed concerning the order of accuracy. The fact that <math display="inline">\mathbf{E}^{r}</math> is not square indicates that for the exact solution only the structural variables can be prescribed, as  pressure can be calculated as a dependent variable.

===2.5 Newmark and exact time integration===

The results up to this point are still general and hold for any time integration scheme with the form of Eqn([[#eq-14|14]]). This allowed us to describe “symbolically” an approximate fractional-step coupling procedure and to evaluate an expression for the truncation error matrix.  In the  following we will consider the “exact” time integrator  and the Newmark time integrator and evaluate  the behaviour of the error  for the two cases.

We highlight that the so called “exact” integrator should  guarantee exact results once the initial conditions and the values of the force at the beginning and end of the step are provided. Applying the ''exact'' time integrator within an ''approximate'' solution procedure  does not guarantee an exact result. The interesting conclusion is that the coupling error obtained with both the Newmark and the exact schemes is very similar and  depends on the coupling procedure rather than on the time integration procedure.

The amplification matrix <math display="inline">\mathbf{A}</math> and the force vectors <math display="inline">\mathbf{L}_n</math> and <math display="inline">\mathbf{L}_{n+1}</math> for the Newmark time integration scheme are widely known <span id='citeF-16'></span>[[#cite-16|[16]]]. In the most general case they take the form

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{A} := \mathbf{A}_1^{-1}\mathbf{A}_2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>     \mathbf{A}_1 :=     \begin{pmatrix}         1+h^2\beta _nK & h^2\beta _nC \\ h\gamma _nK & 1+h\gamma _nC     \end{pmatrix} </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>     \mathbf{A}_2 :=     \begin{pmatrix}         1-\frac{h^2}{2}(1-2\beta _n)K &         h(1-h(1-2\beta _n)\frac{C}{2}) \\ -h(1-\gamma _n)K & 1 -         h(1-\gamma _n)C     \end{pmatrix} </math>
|}
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{L}_n  := \mathbf{A}_1^{-1} \left(     \begin{array}{c}         \frac{h^2}{2}(1-2\beta _n) \\  h(1-\gamma _n)     \end{array} \right)\quad ,\quad  \mathbf{L}_{n+1} :=  \mathbf{A}_1^{-1} \left(     \begin{array}{c}         \frac{h^2}{2}(2\beta _n) \\ h\gamma _n     \end{array} \right) </math>
|}
|}

where <math display="inline">h</math> represents the time step size. We will consider in the following exclusively the second order accurate, non dissipative case.

By assuming that the forces vary linearly within the time step, it is possible to work out a similar result corresponding to the ''exact time integrator''. The amplification matrix becomes now

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{A} := \begin{pmatrix}     cos(\omega h) & \frac{sin(\omega h)}{\omega } \\     -sin(\omega h)\omega & cos(\omega h) \end{pmatrix} </math>
|}
|}

and the force vectors are

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{L}_{n+1} := \begin{pmatrix}     \frac{\omega h-sin(\omega h)}{\omega ^3 h} \\      -\frac{-1+cos(\omega h)}{\omega ^2 h} \end{pmatrix}\quad ,\quad  \mathbf{L}_{n} := \begin{pmatrix}     \frac{-(-sin(\omega h)+cos(\omega h) \omega h)}{\omega ^3 h} \\      \frac{(sin(\omega h) \omega h-1+cos(\omega h))}{(\omega ^2 h)} \end{pmatrix} </math>
|}
|}

Substituting above expressions into Eqn([[#eq-44|44]]) is straight-forward but implies some heavy algebraic calculations. Performing this operation in Maple or a similar Computer Algebra System (CAS) program presents no particular difficulty even for the general case. It is  interesting to focus on the undamped case (no structural damping) as this leads to simpler results. We note that the estimates obtained hold as well for the general damped case

An important (although not unexpected) result is that consistency requires the condition <math display="inline">1 = \eta + \phi </math> for the “free” parameters in the prediction step ( Eqn([[#eq-24|24]]) ). The two choices <math display="inline">\eta=1,\phi=0</math> and <math display="inline">\eta=2,\phi=-1</math> are appealing and correspond to first and second order accuracy for the predicted force, respectively. The order of this initial prediction determines that for the overall approximate time integration scheme for the coupled problem.

Table [[#table-1|1]] contains the first non-zero terms for the Taylor expansion of the error matrix. The analysis of these results shows how the same error estimates hold for the Newmark and exact time integration schemes. This feature is desiderable as it suggests that different choices for the time integration scheme can be taken without affecting the coupling error.

The properties of the coupled fractional step depend however on the values assumed by the aeroelastic mass and the damping term. For <math display="inline">\alpha{>1}</math> or <math display="inline">\beta{>1}</math> the mathematical structure of the problem changes as the overall coupled problem assumes a negative mass matrix. By assuming <math display="inline">\alpha{<1}</math> or <math display="inline">\beta{<1}</math> (see remark [[#eq-12|12]] and [[#eq-13|13]]) we obtain the inequalities

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>     \frac{\alpha ^2}{\alpha{-1}} < 0 \quad ,\quad    \frac{\beta{-1}}{\alpha{-1}} > 0 </math>
|}
|}

which give some insight on the damping properties of the proposed fractional step procedure. 
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>'''Table. 1''' Error matrix for the fractional step coupling procedure as a function of the time integrator, the time step size <math>h</math> and the prediction order
|-
|                    
'''FIRST ORDER PRESSURE PREDICTION''' <math>\eta = 1</math> <math>\phi = 0</math>
|-
| Exact time integration             
|-

|-
| style="text-align: center;" | <math>                 \mathbf{E}_1^{r} :=                 \frac{\alpha ^2}{\alpha{-1}}\frac{\left(\alpha{-\beta}\right)}{\alpha }                 \begin{pmatrix}                     \frac{1}{4}\frac{\beta{-1}}{\alpha{-1}}\omega ^4 h^4     &        -\frac{1}{6}\omega ^2h^3 \\                 \frac{3}{4}\frac{\beta{-1}}{\alpha{-1}}\omega ^4 h^3     &        -\frac{1}{2}\omega ^2h^2 \\                 \frac{3}{2}\frac{\beta{-1}}{\alpha{-1}}\omega ^4 h^2     &        -\omega ^2h \\                 \end{pmatrix}             </math>
|-
|Newmark time integration             
|-
| style="text-align: center;" | <math>             \mathbf{E}_1^{r} :=             \frac{\alpha ^2}{\alpha{-1}}\frac{\left(\alpha{-\beta}\right)}{\alpha }                 \begin{pmatrix}                     \frac{3}{8}\frac{\beta{-1}}{\alpha{-1}}\omega ^4 h^4     &        -\frac{1}{4}\omega ^2h^3 \\                     \frac{3}{4}\frac{\beta{-1}}{\alpha{-1}}\omega ^4 h^3     &        -\frac{1}{2}\omega ^2h^2 \\                     \frac{3}{2}\frac{\beta{-1}}{\alpha{-1}}\omega ^4 h^2     &        -\omega ^2h \\                 \end{pmatrix}             </math>
|-
|
'''SECOND ORDER PRESSURE PREDICTION''' <math display="inline">\eta = 2</math> <math display="inline">\phi = -1</math>
|-
| Exact time integration             
|-
|-
| style="text-align: center;" | <math>             \mathbf{E}_{2}^{r} :=             \frac{\alpha ^2}{\alpha{-1}}\frac{\left(\alpha{-\beta}\right)}{\alpha }             \frac{\beta{-1}}{\alpha{-1}}             \begin{pmatrix}                 \frac{1}{6}\omega ^4 h^4     &        -\frac{1}{6}\omega ^4h^5 \\                 \frac{1}{2}\omega ^4 h^3     &        -\frac{1}{2}\omega ^4h^4 \\             \omega ^4 h^2                            &        -\omega ^4h^3 \\             \end{pmatrix}             </math>
|-
|Newmark time integration             
|-
|-
| style="text-align: center;" | <math>            \mathbf{E}_{2}^{r} :=             \frac{\alpha ^2}{\alpha{-1}}\frac{\left(\alpha{-\beta}\right)}{\alpha }             \frac{\beta{-1}}{\alpha{-1}}                 \begin{pmatrix}                 \frac{1}{6}\omega ^4 h^4     &        -\frac{1}{6}\omega ^4h^5 \\                 \frac{1}{2}\omega ^4 h^3     &        -\frac{1}{2}\omega ^4h^4 \\                 \omega ^4 h^2                            &        -\omega ^4h^3 \\                 \end{pmatrix}             </math>
|}


The exactness of these results was tested using the Newmark scheme to verify that the theoretical predictions match the numerical values. The numerical approach allows as well the direct assessment of the  accuracy order. Figures [[#img-1|1]] and [[#img-2|2]] are obtained for the general case (including non zero structural damping) and contain a log-log plot of the error for the first and second order fractional step schemes for diminishing time step size. It can be seen how the error on the pressure dominates the error. This is not surprising as the pressure term depends on the acceleration. <div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-error_order1_black.png|500px|First order fractional step. log-log error plot.]]
|- style="text-align: center; font-size: 85%;"
| colspan="1" | '''Figure 1:''' First order fractional step. log-log error plot.
|}

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-error_order2_black.png|500px|Second order fractional step. log-log error plot.]]
|- style="text-align: center; font-size: 85%;"
| colspan="1" | '''Figure 2:''' Second order fractional step. log-log error plot.
|}

It is of interest that these estimates hold for the general case. When the aeroelastic mass is negligible the performance of the fractional step procedure is  enhanced and the overall order of accuracy is three. Using the exact time integrator we obtain the error matrix

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}^r := \begin{pmatrix}-\frac{1}{6}\xi \omega ^5(-1+\beta )(\beta{+4}\xi ^2)h^5   &    -\frac{1}{3}\xi ^2\omega ^4(2\beta{+4}\xi ^2-1)h^5 \\     -\frac{1}{2}\xi \omega ^5(-1+\beta )(\beta{+4}\xi ^2)h^4   &    -\xi ^2\omega ^4(2\beta{+4}\xi ^2-1)h^4 \\     -\xi \omega ^5(-1+\beta )(\beta{+4}\xi ^2)h^3 &     -2\xi ^2\omega ^4(2\beta{+4}\xi ^2-1)h^3 \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

===2.6 Stability issues===

To complete the discussion we study the stability of the scheme proposed. Let us consider for simplicity the one-step (first order) fractional step scheme (<math display="inline">\mu =1, \phi =0</math> in Eqn([[#eq-24|24]])). We already know that the exact solution of the coupled problem (including both the pressure and the structural variables) can be written as

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}^{ex}=\mathbf{A}_{ex} \mathbf{z}_{n}^{ex} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

while the approximate solution, obtained via some coupling procedure, can be expressed as

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}=\mathbf{A}_{app} \mathbf{z}_{n} + \mathbf{E} \mathbf{z}_{n}     \, ; \, \mathbf{E} := \mathbf{A}_{app} - \mathbf{A}_{ex} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

where we have introduced a new  error matrix <math display="inline">\mathbf{E}</math> relating the structural variables and the pressure. Some confusion might arise concerning the amplification matrices used. <math display="inline">\mathbf{A}_{ex} \, , \mathbf{A}_{app} </math> are here <math display="inline">3\times 3</math>  matrices and should not be confused with the “reduced” forms (rectangular) used for assessing the accuracy (see [[#table-1|Table 1]]).

A classical discussion of  stability is not trivial due to the complexity of the approximate amplification matrix. We know however that stability is governed by the spectral radius of the amplification matrix. For the one-step fractional step procedure, the eigenvalues can be computed analytically. For the general case the results are  very lengthy and difficult to handle. Some useful information can however be obtained by expressing the eigenvalues as functions of the auxiliary variable <math display="inline"> \Omega :=\omega h </math> which is a non dimensional measure of the  number of steps for integrating of one structural period. By expanding  in Taylor series the expressions obtained for the eigenvalues , the leading  terms of these series (which represent the results for the time step  size <math display="inline">h</math> tending to zero) give:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e_I =     \frac{\alpha{-1}+i\Omega \sqrt{1-\alpha }{\sqrt{\alpha{-2}\beta{+1}}}}{\alpha{-1}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e_{I} = abs\left({\alpha }\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e_{III} =     \frac{\alpha{-1}-i\Omega \sqrt{1-\alpha }{\sqrt{\alpha{-2}\beta{+1}}}}{\alpha{-1}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

where the index represents the eigenvalue number. A necessary requirement for stability is that at the limit for <math display="inline">\Omega =\frac{\omega }{h}</math> tending to zero, the spectral radius <math display="inline">\rho </math> is lower or equal to one. By taking the limit in the simplified expressions we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho _{max} = max\left(\| \alpha \| ,1\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

which suggests that stability is problem dependent and the proposed procedure is ''unconditionally unstable'' for <math display="inline">\alpha < -1</math>.

This definition of stability is however not completely appropriate for the coupling problem of interest. For the FSI coupling both divergent and convergent solutions are possible. This implies that the exact amplification matrix may be characterized by eigenvalues of value greater than one. When, as in the case we are treating here, no damping is considered, the target value for the spectral radius is exactly one. A value greater than 1 (but tending to 1 for diminishing step size) may be accepted, even if it corresponds to a numerically divergent solution. When the system is characterized by <math display="inline">\alpha < -1</math>, even for little time steps the solution is divergent and the situation can not be improved by reducing the time step. This leads to explosive instabilities.

===2.7 Case of negligible aeroelastic mass===

It is interesting to consider the behaviour of the algorithm for <math display="inline">\alpha </math> tending to zero. When the aeroelastic mass disappears the equations simplify and one order of accuracy can be gained for the second order prediction case.

When the aeroelastic damping is zero the equations simplify even further and the stability conditions can be investigated analytically (for the first order prediction case). Under these conditions the procedure is unconditionally stable if <math display="inline">\beta{>}-1</math> and conditionally stable otherwise. For <math display="inline">\beta{<}-1</math> the stability condition is <math display="inline">h \leq \frac{2}{\omega \sqrt{-1-\beta }}</math>.

The general damped case, on the other hand, leads to very long expressions which are not well suited for analytical treatment. The spectral radius can now be expressed as <math display="inline"> \rho (n_{div},\xi _{fsi},\beta ) </math> where <math display="inline">n_{div} = \frac{2\pi }{\Omega }</math> represents the number of time steps used for discretizing the natural period of the structure alone. Useful information concerning the behaviour of the procedure of interest can be obtained by plotting  the spectral radius versus <math display="inline"> \xi _{fsi},\beta </math> for an increasing number of divisions of the structural period.

Figure [[#img-3|3]] shows  how the approximate spectral radius (<math display="inline">\rho _{approx}</math>) is always less than one when the aeroelastic damping constant is negative (which implies a positive damping for the system). It can be also seen how the stability region extends to the positive damping range for very low values of <math display="inline">n_{div}</math>. This suggests that the algorithmic damping is strong when the number of divisions is  low. For high values of the number of divisions the isolines of <math display="inline">\rho _{approx}</math> become straight and parallel which testifies how the additional algorithmic damping vanishes.

The plots in Figure [[#img-4|4]] show a measure of the algorithmic damping versus the numerical one. The scale of the graphs changes as the error vanishes very quickly. However the plots indicate the areas of the <math display="inline">\xi _{fsi},\beta </math> plane where the algorithmic damping is greater or lesser than that for the Newmark scheme applied to the exact case.

<div id='img-3a'></div>
<div id='img-3b'></div>
<div id='img-3c'></div>
<div id='img-3d'></div>
<div id='img-3e'></div>
<div id='img-3f'></div>
<div id='img-3g'></div>
<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-rho_ndiv2.png|275px|Plot of ρ<sub>approx</sub> - 1  n<sub>div</sub> = 2  ]]
|[[Image:Draft_Samper_231103336-rho_ndiv4.png|275px| n<sub>div</sub> = 4   ]]
|- style="text-align: center; font-size: 75%;"
| (a) <math> n_{div} = 2</math> 
| (b) <math display="inline"> n_{div} = 4  </math> 
|-
|[[Image:Draft_Samper_231103336-rho_ndiv10.png|275px| n<sub>div</sub> = 10  ]]
|[[Image:Draft_Samper_231103336-rho_ndiv20.png|275px| n<sub>div</sub> = 20   ]]
|- style="text-align: center; font-size: 75%;"
| (c) <math display="inline"> n_{div} = 10  </math>
| (d) <math display="inline"> n_{div} = 20  </math> 
|-
|[[Image:Draft_Samper_231103336-rho_ndiv40.png|275px| n<sub>div</sub> = 40  ]]
|[[Image:Draft_Samper_231103336-rho_ndiv80.png|275px| n<sub>div</sub> = 80   ]]
|- style="text-align: center; font-size: 75%;"
| (e) <math display="inline"> n_{div} = 40  </math>
| (f) <math display="inline"> n_{div} = 80  </math> 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 3''' Plot of <math>\rho_{approx} - 1</math>
|}

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-a05b0.png|550px|Error plot: Aitken Acceleration - red line,     “Consistent Acceleration” - green line. α=0.5 \, ,     \, β= 0.0]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Error plot: Aitken Acceleration - red line,     “Consistent Acceleration” - green line. <math>\alpha=0.5 \, ,     \, \beta = 0.0</math>
|}

==3 Iterative enhancement of the fractional step scheme==

An appealing possibility for improving the fractional step procedure previously explained is to attempt improving the properties of the method by iterating through the different domains. Interestingly the analysis framework introduced can be applied directly to such an iterative improvement and allows us to describe analytically the behaviour of the error during the iterations.

This is done in symbols by replacing the predicted value in Eqn([[#eq-33|33]]) with the last known approximation of the system variables, which gives

<span id="eq-53"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{z}_{n+1}^{i+1} = \mathbf{B}_0 \mathbf{z}_{n+1}^{i} + \mathbf{B}_1 \mathbf{z}_{n} + \mathbf{B}_2 \mathbf{z}_{n-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

where the upper index “<math>i</math>” identifies the iteration count. Note that  using the last known values corresponds to a Jacobi type iteration and may not converge in some situations.

Following the theoretical framework proposed above, the error matrix corresponding to the (unrelaxed) iterative application of the correction step can be directly assessed. The result is simple and interesting. By identifying as <math display="inline">\mathbf{E}^{r}</math> the “reduced” error matrix containing the leading terms in the series expansion for the error of the first order fractional step, the following relations hold

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}_{0}^{r} = \frac{1}{\alpha } \mathbf{E}^{r} \, \, \hbox{(explicit prediction)} </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}_{1}^{r} =  \mathbf{E}^{r} \, \, \hbox{(end of basic fractional step procedure)} </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}_{2}^{r} = \alpha \mathbf{E}^{r} </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}_{3}^{r} = \alpha ^2 \mathbf{E}^{r} </math>
|}
|}

This result was verified using a CAS (Computer Algebra System) for the first 4 iterations. By induction we can assume

<span id="eq-54"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}^r_{i} = \alpha ^{i-1} \mathbf{E}^r </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

This result is important as it states that the simple Jacobi iteration will not converge when <math display="inline">\alpha < -1</math>. This corresponds to a case in which the fluid is much heavier than the structure or more correctly a case for which the added mass effect is crucial. Under these conditions an accelerator is needed to ensure convergence. This behavior is well known in many practical FSI applications.

It is also important that for values of <math display="inline">\alpha \approx -1</math>, even on the “safe side”, the simple iterative strategy will converge very slowly. Fast convergence will  be achieved for cases in which the added mass effect is not too important. The conclusion is  that the basic iterative strategy works well only when “explicit” schemes are effective, which indeed reduces its practical importance.

In practice a relaxation factor <math display="inline">\omega </math> is generally introduced to guarantee convergence. The optimal value for <math display="inline">w</math>  can be found by trial and error or eventually using an acceleration technique. Mok and Wall <span id='citeF-1'></span>[[#cite-1|[1]]],<span id='citeF-2'></span>[[#cite-2|[2]]] for example advocate the use of the Aitken accelerator for FSI problems as its the performance appears to be comparable to the best gradient accelerated techniques. No theoretical justification is however given for this behavior. 

Interestingly, the optimal value of the relaxation factor depends on the aeroelastic mass. To show this let us consider first the definition of error

<span id="eq-55"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{y}^i_{n+1}-\mathbf{y}^{ex}_{n+1} \approx \alpha ^{i-1}\mathbf{E} \mathbf{y}_{n} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

by subtracting the solution between successive iterations and pre-multiplying by a given vector <math display="inline">\mathbf{v}^T</math>  we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{v}^T\left(\mathbf{y}^i_{n+1}-\mathbf{y}^{i-1}_{n+1}\right)\approx      \alpha ^{i-2}\left(\alpha{-1}\right)\mathbf{v}^T\mathbf{E} \mathbf{y}_{n} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{v}^T\left(\mathbf{y}^{i-1}_{n+1}-\mathbf{y}^{i-2}_{n+1}\right)\approx      \alpha ^{i-3}\left(\alpha{-1}\right)\mathbf{v}^T\mathbf{E} \mathbf{y}_{n} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

The ratio between the two above equations gives finally an estimate for the aeroelastic mass coefficient as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha = \frac{\mathbf{v}^T\left(\mathbf{y}^i_{n+1}-\mathbf{y}^{i-1}_{n+1}\right)}{\mathbf{v}^T\left(\mathbf{y}^{i-1}_{n+1}-\mathbf{y}^{i-2}_{n+1}\right)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

This estimate is useful for computing the “optimal” relaxation parameter. By definition the relaxation takes the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{y}^* = \omega \mathbf{y}^{i} + \left(1-\omega \right)\mathbf{y}^{i-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

Using Eqn([[#eq-55|55]]) and simplifying we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{y}^* \approx \mathbf{y}^{ex} + \left(\omega +     \left(\omega{-1}\right)\alpha \right)\mathbf{E} \mathbf{y}_n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|}

which allows us to obtain the optimal relaxation factor as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\omega = \frac{\alpha }{1+\alpha }=  \frac{\frac{\mathbf{v}^T\left(\mathbf{y}^i_{n+1}-\mathbf{y}^{i-1}_{n+1}\right)}{\mathbf{v}^T\left(\mathbf{y}^{i-1}_{n+1}-\mathbf{y}^{i-2}_{n+1}\right)}}{ 1 + \frac{\mathbf{v}^T\left(\mathbf{y}^i_{n+1}-\mathbf{y}^{i-1}_{n+1}\right)}{\mathbf{v}^T\left(\mathbf{y}^{i-1}_{n+1}-\mathbf{y}^{i-2}_{n+1}\right)}}     =\frac{\mathbf{v}^T\left(\mathbf{y}^i_{n+1}-\mathbf{y}^{i-1}_{n+1}\right)}{ \mathbf{v}^T\left(\mathbf{y}^{i-1}_{n+1}-\mathbf{y}^{i-2}_{n+1}\right)+   \mathbf{v}^T\left(\mathbf{y}^i_{n+1}-\mathbf{y}^{i-1}_{n+1}\right)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

or simplifying.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\omega = \frac{\mathbf{v}^T\left(\mathbf{y}^i_{n+1}-\mathbf{y}^{i-1}_{n+1}\right)}{ \mathbf{v}^T\left(\mathbf{y}^{i}_{n+1}-\mathbf{y}^{i-2}_{n+1}\right)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

This expression allows us to reduce to zero the leading term of the error (on the SDOF test system) which leads to a higher order of accuracy in time.

This estimation procedure delivers results which are very similar to the so called Aitken acceleration (see for example Figure([[#img-4|4]])) which anticipates the success in its application to FSI.

===3.1 On the use of dissipative algorithms for the structural integration===

All the results so far shown refer to an exact or non-dissipative (Newmark) time integration scheme. In many problems of practical interest it is important the use of a dissipative solver to diminish the noise connected to the high frequencies, or simply to guarantee the stability of the time integrator for non linear problems. Unfortunately the assumption in Eqn([[#eq-14|14]]) does not hold for arbitrary dissipative schemes making the analysis not valid for such cases.

Undesiderable  surprises may arise due to a perverse combination of errors between the coupling iterations and the structural solution. This situation can be reproduced easily on the test model problem by choosing a rhs in the form <math display="inline">RHS=\alpha \ddot x</math> and using a Bossak time integration scheme. Even if this aspect is not studied in detail in this work we believe this is a interesting problem to be addressed in the future.

Fig. [[#img-5|5]] shows results of the iterative coupled scheme vs the ones obtained using the same integrator and a monolithic approach. It is clear how the iterative solution exhibits some spurious damping which is not present when a non-dissipative algorithm is used. This spurious damping appears to be proportional to the aeroelastic mass and disappears as the time step tends to zero. As a consequence, the importance of this effect is reduced as long as the mass of the structure is sensibly greater than that of the surrounding fluid. However the phenomena may become of major concern in dealing with heavy fluids.

We would like to point out that a very simple modification allows to use the fractional step procedure in combination with a dissipative algorithm. The accuracy can be recovered by using a undamped solver for the structural prediction while keeping the numerically damped solver for the  correction step.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-Bossak01a05.png|550px|Comparison between exact and iterative coupled     solutions using the Bossak algorithm.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Comparison between exact and iterative coupled     solutions using the Bossak algorithm. Bossak parameter = -0.1,     period divided in 50 steps - Monolithic solution is dashed while the     solution of iterative coupling is represented with a continuum line
|}
{|
|- style="text-align: center; font-size: 75%;"
|'''Notes:''' Bossak parameter = -0.1,     period divided in 50 steps - Monolithic solution is dashed while the     solution of iterative coupling is represented with a continuum line
|}

==4 An improved iterative coupling procedure via a modified equations==

As discussed above the simple iterative strategy described in Section [[#3 Iterative enhancement of the fractional step scheme|3]] may encounter convergence difficulties or even diverge depending on the features of the flow. Even relaxed iterative strategies may fail for complex FSI problems when the aeroelastic mass term becomes very important.

This can be justified from an heuristic point of view  by considering that most relaxation strategies perform a relaxation on the displacements term which in turn reflects on some form of control over the acceleration dependent contributions. Given the high dependency of the acceleration on small variations of the displacements, this turns out to be ineffective to control the (aeroelastic) mass-induced oscillations.

A natural solution for this problem is to control the variation of the acceleration in order to improve the convergence. The idea is better expressed in mathematical terms. The differential equations governing the dynamic equilibrium, follow a basic iterative scheme of the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M_s  \ddot{x} _{n+1}^{i+1} + D_s  \dot{x} _{n+1}^{i+1} + K_s     x_{n+1}^{i+1} = p_{n+1}^{i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
|}

This equation can be modified symbolically to give

<span id="eq-64"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M_s  \ddot{x} _{n+1}^{i+1} - \alpha M_s      \ddot{x} _{n+1}^{i+1}+ D_s  \dot{x} _{n+1}^{i+1} + K_s     x_{n+1}^{i+1} =  p_{n+1}^{i}- \alpha M_s  \ddot{x} _{n+1}^{i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|}

where <math display="inline">\alpha </math> is a scalar parameter. Clearly as <math display="inline">\left\| \ddot{x} ^{i+1}- \ddot{x} ^i \right\| \rightarrow  0</math> the solution of this modified equation model tends to the original solution. Nevertheless the parameter <math display="inline">\alpha </math> affects the convergence of the modified problem converges to the coupled solution.

Once again the simple model problem  help us  in investigating the convergence properties for this modified equation. We already know that the case <math display="inline">p = M_{f} \ddot{x}  </math> is the main responsable of the convergence difficulties for the iterative strategy. We will therefore develop our analysis on the simpler undamped case

<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M_s  \ddot{x} _{n+ 1}^{i+1} + K_s     x_{n+1}^{i+1} = \alpha ^{ex}  M_s \ddot{x} _{n+1}^{i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|}

where <math display="inline">\alpha ^{ex}</math> is the exact value of <math display="inline">\alpha </math> giving the exact aeroelastic mass (Eqn([[#eq-65|65]])) has the same convergence properties as the more general case. The modified problem becomes

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(M_s -\alpha  M_s\right) \ddot{x} _{n+1}^{i+1} + K_s     x_{n+1}^{i+1} = \left(\alpha ^{ex}-\alpha \right)     M_s \ddot{x} _{n+1}^{i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}

by dividing by <math display="inline">1-\alpha ^i</math> this can be recasted as

<span id="eq-67"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M_s  \ddot{x} _{n+1}^{i+1} + K_s     x_{n+1}^{i+1} = \frac{\alpha ^{ex}-\alpha }{1-\alpha }     M_s \ddot{x} _{n+1}^{i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

were  we modified accordingly the stiffness term (which does not play a role in the stability). The system assumes now the form of Eqn([[#eq-65|65]]), of which it inherits the same properties. The convergence condition becomes

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left|\frac{\alpha ^{ex}-\alpha }{1-\alpha }\right|< 1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
|}

As observed above (see Eqn([[#eq-12|12]])) <math display="inline">\alpha _{ex}<1</math>. As <math display="inline">\alpha </math> is an approximation of <math display="inline">\alpha _{ex}</math> we can therefore that <math display="inline">\alpha{<1}</math>. Under this assumption the unequality can be solved to give

<span id="eq-69"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha <     \frac{1}{2}\alpha ^{ex}+\frac{1}{2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
|}

Basing on the previous results  we  conclude that for a system in the form of Eqn([[#eq-67|67]]) the error varies with iterations as <math display="inline">\mathbf{E}^r_{i} = \left(\frac{\alpha ^{ex}-\alpha }{1-\alpha }\right)^{i-1}\mathbf{E}^r_0</math>. It follows that the optimal choice for the modified equation method is given by <math display="inline">\alpha = \alpha ^{ex}</math>, for which the spurious oscillations are immediately damped out.

The  parameter <math display="inline">\alpha </math> affects sensibly the condition for stability, hereby  increasing the area of convergence. Any choice of the type <math display="inline">\alpha < \alpha ^{ex}</math> guarantees stability and monothonic convergence (<math display="inline">0<\frac{\alpha ^{ex}-\alpha }{1-\alpha }<1</math>) while the choice <math display="inline">\alpha ^{ex}<\alpha <\frac{1}{2}\alpha ^{ex}+\frac{1}{2}</math> guarantees stability but gives oscillatory convergence. Unfortunately the exact value for <math display="inline">\alpha </math> is not known. However Eqn([[#eq-69|69]]) ensures that the convergence is robust and that a rough estimation of the optimal value of <math display="inline">\alpha </math> is sufficient to ensure convergence.

Let us make now some heuristic considerations on the properties of the modified equation approach proposed. The crucial point is that at every iteration we solve a problem that is different from the original one and which solution becomes increasingly similar to that of the original problem as the iterative process converges. Even at convergence, an error exists which implies that the solution found corresponds to a slightly different problem. The choice of an error tolerance therefore tell us how closely the modified equations will “mimic” the real problem. The key is however to guarantee ''“a priori”'' the stability of the coupled procedure while improving the accuracy as needed. Note that this is not the case for a general relaxation scheme.

===4.1 Estimating alpha===

The real interest of any acceleration scheme is of course in application  to real (non 1D) models. For such problems we can observe that the interface force stemming from the pressure contribution should be oriented as the normal to the interface and proportional to the surface area, while an heuristic consideration tells us that only the acceleration normal to the interface will play a role in the coupling. This suggests that <math display="inline">\mathbf{M}_{f}</math> should take, for each interface node, a form of the type <math display="inline">\mathbf{M}_{f}=\alpha \mathbf{A}_n \mathbf{A}_n^T</math> where <math display="inline">\mathbf{A}_n</math> represents the normal at the interface node multiplied by the corresponding influence area.

Several possibilities exist for the practical estimation of a coefficient <math display="inline">\alpha </math> to be used in improving the convergence of the coupled solver. A good estimation could be derived by physical reasoning or by considering the actual choice of the flow solver.

An alternative procedure, which makes use exclusively of the data at different iterations is derived next for the non-scalar case. A general iterative scheme for the estimation of such parameter may lead to the use of a different estimate of <math display="inline">\alpha </math> at each iteration “<math>i</math>”. If we assume however the use of a constant <math display="inline">\alpha ^i</math> through the iteration process, the application of Eqn([[#eq-67|67]]) to two successive iterations gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M}\mathbf{ \ddot{x} }_{n+1}^{i+1} = \frac{\alpha ^{ex}-\alpha }{1-\alpha ^i}     \mathbf{M}\mathbf{ \ddot{x} }_{n+1}^{i} - \frac{1}{1-\alpha ^i}\mathbf{K}\mathbf{x}_{n+1}^{i+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M}\mathbf{ \ddot{x} }_{n+1}^{i} = \frac{\alpha ^{ex}-\alpha }{1-\alpha ^i}     \mathbf{M} \mathbf{ \ddot{x} }_{n+1}^{i-1} - \frac{1}{1-\alpha ^i}\mathbf{K}\mathbf{x}_{n+1}^{i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
|}

Note that we are  considering here the real, multidimensional problem, it follows that <math display="inline"> \ddot{x} </math>, <math display="inline">\mathbf{x}</math> and <math display="inline">M</math> refer to the acceleration, displacement and mass associated to the interface degrees of freedom.

Subtracting the above equations and pre-multiplying by <math display="inline">\mathbf{M}^{-1}</math> gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{ \ddot{x} }_{n+1}^{i+1} - \mathbf{ \ddot{x} }_{n+1}^{i} =     \frac{\alpha ^{ex}-\alpha }{1-\alpha ^i}     \left(\mathbf{ \ddot{x} }_{n+1}^{i}-\mathbf{ \ddot{x} }_{n+1}^{i}\right)-     \frac{1}{1-\alpha ^i}\mathbf{M}^{-1}\mathbf{K}\left(\mathbf{x}_{n+1}^{i+1}-\mathbf{x}_{n+1}^{i}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
|}

For small time steps the accelerations governs the phenomena (a small variation in displacements in a small time induces a high acceleration to rise). This allows to neglect the term <math display="inline">\mathbf{M}^{-1}{\boldsymbol K} \left(\mathbf{x}_{n+1}^{i+1}-\mathbf{x}_{n+1}^{i}\right)</math>. Under this assumption, premultiplying both sides by (<math display="inline"> \mathbf{ \ddot{x} }_{n+1}^{i}-\mathbf{ \ddot{x} }_{n+1}^{i}</math>) we obtain the scalar relation

<span id="eq-73"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>c = \frac{\alpha ^{ex}-\alpha }{1-\alpha ^i}     \quad \hbox{and}\quad \alpha ^{ex}=c(1-\alpha ^i)+\alpha  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>c := \frac{\left(\mathbf{ \ddot{x} }_{n+1}^{i+1} -     \mathbf{ \ddot{x} }_{n+1}^{i}\right)^T\left(\mathbf{ \ddot{x} }_{n+1}^{i+1} -     \mathbf{ \ddot{x} }_{n+1}^{i}\right)}{\left(\mathbf{ \ddot{x} }_{n+1}^{i+1} -     \mathbf{ \ddot{x} }_{n+1}^{i}\right)^T\left(\mathbf{ \ddot{x} }_{n+1}^{i} -     \mathbf{ \ddot{x} }_{n+1}^{i-1}\right)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
|}

Eqn([[#eq-73|73]]) can be solved for <math display="inline">\alpha ^{ex}</math> and used for calculating <math display="inline">\alpha ^{i+1}:=\alpha ^{ex}</math>

This procedure is only one of the many possibilities for estimating <math display="inline">\alpha </math>. In many cases it may be attractive to use a value obtained by other sort or reasoning. In any case, once <math display="inline">\alpha </math> has been estimated (by any means) the relaxation method can be easily implemented with a minor modification to the structural solution code.

===4.2 A simple example to verify the implementation===

We choose a simple example  to verify the correctness of the implementation. This is obtained by taking a 1m cube of elastic material subjected to a variable pressure on one side.

The external pressure is made dependent on the acceleration following a law of the type <math display="inline">p = M_{f} \ddot{x} _{norm}</math>. For this numerical setting, the proposed method is expected to recover exactly the aeroelastic matrix and, apart for an initial transient phase, to reproduce the analytical result in just one iteration (two iterations are actually needed for convergence to be detected).

Taking <math display="inline">\rho = 1000 \frac{Kg}{m^3}</math>, <math display="inline">E = 5e9 \frac{N}{m^2}</math>, <math display="inline">\nu = 0</math>, <math display="inline">t = 1m</math> and using a lumped mass matrix the actual system to be solved is equivalent to the scalar problem

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{1}{2} M_s  \ddot{x} + \frac{EA}{L} x =  M_{f}  \ddot{x}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
|}

Assuming <math display="inline">M_{f} = -1000</math>, for which the standard staggered approach would not converge, the system takes the following form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>500  \ddot{x} + 5000000 x = 1000  \ddot{x}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
|}

Fig. [[#img-6|6]] shows a comparison between the exact and analytical results. Convergence is obtained in two iterations except during the first three time steps. A plot of the error norm <math display="inline">\left\| \ddot{x} ^{i+1}- \ddot{x} ^i\right\|</math> is given in Fig. [[#img-7|7]]

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-AnalyticalVsNumerical.png|600px|Comparison between analytical and numerical solution     for the cube model problem]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Comparison between analytical and numerical solution     for the cube model problem
|}

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-errorplot.png|550px|Acceleration error for the simple cube benchmark     scheme]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Acceleration error for the simple cube benchmark     scheme
|}

==5 Restrictor Flap==

This test was proposed by Wall and Mok <span id='citeF-1'></span>[[#cite-1|[1]]]  and Mok <span id='citeF-2'></span>[[#cite-2|[2]]], as a challenging test for the stability of solution algorithm procedures. The geometry is described in Fig. [[#img-8|8]]. The data used for the analysis are <math display="inline">\rho _F = 956 Kg/m^3</math>, <math display="inline">\mu = 0.145 Kg/ms</math>, <math display="inline">\rho _S = 1500 Kg/m^3</math>, <math display="inline">E = 2.3 MPa</math>, <math display="inline">\nu = 2.3 MPa</math>. The resulting pressure histories are shown in Fig. [[#img-9|9]] and directly compared with the results of the literature. The frames shown in Fig. [[#img-10|10]] show the evolution of the coupled fluid-structure deformation and compare well with the results of the same example as reported in Mok <span id='citeF-2'></span>[[#cite-2|[2]]]. A second order fractional step solver is used for solving the problem. The resulting pressure histories are shown in Fig. [[#img-9|9]] and directly compared with the results of the literature

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-mesh_restrictor.png|600px|Restrictor flap - mesh of 7928 elements]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Restrictor flap - mesh of 7928 elements
|}

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-pressure_comparison.png|400px|Pressure vs time at mid-height (A) and at the tip     (B)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Pressure vs time at mid-height (A) and at the tip     (B)
|}
{|

<div class="center" style="font-size: 75%;">
'''Notes:''' Results of pressure vs time are shown on the top of the curves found in the literature
</div>

<div id='img-10a'></div>
<div id='img-10b'></div>
<div id='img-10c'></div>
<div id='img-10d'></div>
<div id='img-10e'></div>
<div id='img-10f'></div>
<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-restrictor5.png|450px|t =     5s]]
|- style="text-align: center; font-size: 75%;"
| (a) <math>t = 5s</math> 
|-
|[[Image:Draft_Samper_231103336-restrictor10.png|450px|t = 10s]]
|- style="text-align: center; font-size: 75%;"
|(b) <math>t = 10s</math> 
|-
|[[Image:Draft_Samper_231103336-restrictor15.png|450px|t =     15s]]
|- style="text-align: center; font-size: 75%;"
| (c) <math>t =     15s</math>
|-
|[[Image:Draft_Samper_231103336-restrictor20.png|450px|t =     20s]]
|- style="text-align: center; font-size: 75%;"
| (d) <math>t =     20s</math>
|-
| colspan="1"|[[Image:Draft_Samper_231103336-restrictor25.png|450px|t =     25s]]
|- style="text-align: center; font-size: 75%;"
|  colspan="1" | <math>(e) t =     25s</math>
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10''' Analysis of a restrictor flap
|}
<div class="center" style="font-size: 75%;">
'''Notes:''' Countours of velocity. Symmetry boundary conditions is used on the top boundary
</div>

==6 Flag flutter==

An interesting validation example for the iterative solution scheme of Section [[#5 Restrictor Flap|5]] is that of a square bluff body followed by a flexible cantilever plate as shown in Fig. [[#img-11|11]]. The prescribed distribution of velocities is shown in the same figure. The numerical setup is reproduced following <span id='citeF-8'></span>[[#cite-8|[8]]]. The fluid characteristics are <math display="inline">\rho _f = 1.18\frac{Kg}{m^3}, \mu=1.8e-5\frac{Kg}{m s} </math> the structure has a elastic modulus of <math display="inline">200000\frac{N}{m^2}</math> and a density <math display="inline">\rho _{s}=2000\frac{Kg}{m^3}</math> with the first three natural frequencies at 0.61,3.8 and 10.6 Hz.

By applying a constant velocity at the inflow <math display="inline">U_{inflow}=0.315\frac{m}{s}</math>, the bluff body sheds Von Karman vortices with a frequency  of 3.7 Hz. (on the rigid system) close to the second natural frequency of the structure. In the rigid case the pressure oscillations are relatively small. Consequently only small displacement motion is expected due to the vortex shedding.

Given the proximity between the shedding frequency and the structural frequency we can expect a large amplitude “resonant” motion. After a first phase of growth, the amplitude of vibration is however expected to reach a stable solution with a peak amplitude of around 0.08m.

Loose coupling procedures based on high order prediction steps (see <span id='citeF-5'></span>[[#cite-5|[5]]]) are known to be unstable for this example. On the other hand, procedures based on a first order prediction step enforce poorly the energy conservation at the interface, leading to an important numerical positive damping (unless the time step is very small). This example was in fact originally proposed as an argument for the necessity of implicit procedures for non-linear aeroelasticity problems. The results shown  were however obtained using the loosely coupled fractional step procedure described in this work. The Fourier transform of the displacement history is given in Fig. [[#img-12|12]]. The leading peak is found at a frequency of 3.04Hz in good agreement with the value of 3.1Hz published in <span id='citeF-8'></span>[[#cite-8|[8]]]. The displacement history of the tip is reported in Fig. [[#img-13|13]] on top of the history reported in <span id='citeF-9'></span>[[#cite-8|[8]]]. The excellent agreement found proves the effectiveness of our scheme  for challenging simulations.


<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-flag_view.png|550px|Flexible plate behind a square bluff body - view of     the domain]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Flexible plate behind a square bluff body - view of     the domain
|}
<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-fftanalysis.png|550px|Flexible plate, Fourier analysis of tip displacement]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Flexible plate, Fourier analysis of tip displacement
|}

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-tipdisp.png|503px|Flexible plate, tip displacement history - Numerical results]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Flexible plate, tip displacement history - Numerical results
|}

==7 Application to the aeroelastic analysis of a bridge==

We propose a third application for the simplest (and fastest) fractional step algorithm proposed in the first part of the paper. The goal of this example is twofold: from one side our interest is in showing that the algorithm actually works. Also prove that the coupling error is dominated by the fluid solution phase rather than by the coupling procedure.

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-A2_A3_H2_H3.png|600px|Fitting of experimental aeroelastic derivatives,     using a first order accurate solver for the fluid domain]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Fitting of experimental aeroelastic derivatives,     using a first order accurate solver for the fluid domain
|}
<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-a2_a3_h2_h3_2nd.png|600px|Free rotation fitting of experimental aeroelastic     derivatives using a second order accurate solver for the fluid     domain]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Free rotation fitting of experimental aeroelastic     derivatives using a second order accurate solver for the fluid     domain
|}

To do so we will focus our attention on the calculation of the flutter derivatives for a long span bridge deck (The Great Belt bridge section) under a constant wind speed. First we present briefly the engineering theory, then we review two methods used for the experimental determination of the coefficients needed, and finally we will show that in comparison to experimental and analytical results the fractional step procedure delivers satisfactory results.

A well established theory of bridge aerodynamics <span id='citeF-17'></span>[[#cite-17|[17]]] was developed in the early 1970s and is still used in engineering applications. In a modern revision <span id='citeF-18'></span>[[#cite-18|[18]]] the theory  is based on the assumption of a mathematical model in the form

<span id="eq-77"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L(t) = \rho \,{U}^{2}B \left({\frac{K{     H_1}\,\dot y}{U}}+{\frac{K{H_2}\,B{ \dot      \alpha }}{U}}+{K}^{2}{H_3}\,\alpha{+}{\frac{{K}^{2}{     H_4}\,y}{B}} \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
|}

<span id="eq-78"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M(t) = \rho \,{U}^{2}{B}^{2} \left({\frac{K{     A_1}\,\dot y}{U}}+{\frac{K{A_2 }\,B{\dot      \alpha }}{U}}+{K}^{2}{A_3}\,\alpha{+}{\frac{{K}^{2}{     A_4}\,y }{B}} \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
|}

where <math display="inline">K=\frac{\omega B}{U}</math>. Eqns([[#eq-77|77]]) and ([[#eq-78|78]]) express the forces exerted by the fluid on the structure due to the motion of the latter. The coefficients in these equations are known as “aeroelastic derivatives” and are obtained experimentally through the analysis of sectional models. Interestingly enough this model can be shown to be equivalent to the model <math display="inline">A(\omega )\mathbf{v} + B(\omega )\mathbf{x}  </math> which we assumed for the analysis of our coupling algorithms.

The  most frequent approach in the numerical study of the bridge aerodynamic problems is focused on the direct determination of the flutter limit <span id='citeF-19'></span><span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-22'></span>[[#cite-19|[19,20,21,22]]]. This is generally measured through the direct assesment of the system's energy, or possibly by the analysis of the displacement history. Unfortunately many factors come into play in this process making difficult the identification of the error sources. Even if in general these errors are on the numerical side it is interesting that even the experimental results may be affected by factors such as  the  correspondence of the numerical Reynolds to the real one, or  the blockage ratio for the model tested.

In  practice, at least two techniques exist for the experimental determination of the flutter derivatives. One based on the prescription of the section motion and on the analysis of the history of the forces applied from the fluid to the structure. The second one is based on the analysis of the displacement history of the section model mounted on springs. Details on the procedure used for the identification of the aeroelastic derivatives in the two cases are given in the Appendix.

Both of those experimental procedures can be reproduced numerically. The forced vibration approach does not involve any coupling error and can be considered as a reference solution for evaluating  the equivalent coupled approach. Its performance is governed exclusively by the properties of the fluid solver and by the mesh used for the computation.

In the hypothesis that the same initial solution, mesh and boundary conditions are used, the aeroelastic derivatives computed using a “free vibration” coupled FSI approach can be directly compared to the ones obtained from the forced analysis. This comparison allows us to isolate the coupling error from the other error sources thus reaching the desired objective.

A similar result could be achieved by direct comparison with the experimental results once the ability of the fluid solver to describe correctly the phenomena is guaranteed. This would require fine meshes greatly increasing the computational cost, thus making the testing procedure unsuitable for the study of different coupling strategies. On the other hand by comparing two numerical results obtained under the same condition we can directly asses the impact of the coupling on the final results, which is the objective of current work.


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>'''Table. 2''' Description of the validation procedure
|-
|          

* Construct the model and compute a starting solution by keeping the section fixed
* Prescribe the motion of the section and identify the aeroelastic derivatives for different frequences
* Suspend the section on springs calculating the spring so to be in the range of reduced velocities of interest.
* Calculate the aeroelastic derivatives in free motion
* Evaluate the coupling properties by comparing the two curves
|}

===7.1 Numerical setup===

The identification techniques described above were used for the identification of the aeroelastic derivatives of a bridge cross-section. Due to the availability of validation data the cross-section of the great belt suspension span in Denmark was chosen for the benchmark exercise. This section, was studied numerically and experimentally <span id='citeF-20'></span><span id='citeF-23'></span><span id='citeF-24'></span>[[#cite-20|[20,23,24]]]. The experimental plots used for the comparison were taken from <span id='citeF-25'></span>[[#cite-25|[25]]].

The analysis domain, as represented in Fig.[[#img-17|17]], was chosen so to keep around 6% the ratio between the section chord and the vertical dimension of the domain in order to minimize the blockage effect. An unstructured mesh composed of 9197 nodes and 17862 triangular 3-noded elements with a minimum element size (around the bridge deck) of 0.15m was used.

Only pressure forces are transferred to the structure, negleting the contribution of viscous forces (which would be poorly approximated by a mesh as the one described). This approximation is commonly accepted for bridges which behave as bluff bodies.

The actual dimensions of the cross-section were used with a Reynolds number that was consequently far greater than the one reproduced in wind tunnel testing. Even if this could theoretically lead to important discrepancies with the experimental data, it is currently believed that for the case of interest the derivatives are fairly independent of the Reynolds number of the flow. The inflow velocity was therefore taken as 10m/s corresponding to a <math display="inline">(Re = \frac{v_{\infty }B}{\nu }\approx  2*10^7</math>).

No turbulence model was included in the simulation however a subgrid scale stabilization method was used [].

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_231103336-mesh_picture.png|600px|Cross section of Great Belt suspension bridge. View of the mesh used in the computations. ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Cross section of Great Belt suspension bridge. View of the mesh used in the computations. 
|}

The structural displacements were kept small with a maximum pitching angle of 3deg and a maximum motion in heave of 1.57m in accordance to the values reported in <span id='citeF-24'></span><span id='citeF-17'></span>[[#cite-24|[24,17]]]. The reduced frequency <math display="inline">\frac{U_{inflow}}{nB}</math>, was varied by changing the value of the frequency <math display="inline">n</math> of the prescribed sectional motion. This allowed to start all the simulations exactly from the same initial conditions. An initial velocity was prescribed to the structure in order to study the decay or growth of the initial solution. This velocity was chosen so to provide a maximum displacement “in vacuo” corresponding to the one for the prescribed motion setting. This was enforced by considering that in absence of fluid the system is conservative, from which follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{1}{2}Mv_{in}^2 = \frac{1}{2}Kx_{max}^2 \rightarrow      v_{in} = \sqrt{\frac{K}{M}}x_{max} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
|}

This setup was used for the assessment of the partitioned coupling procedures implemented in Kratos, a Multi-Physics code developed at CIMNE (www.cimne.com/kratos) . Both first and second order solvers were used, the latter being vastly superior in reproducing the experimental curves. For the first order case, showed in Fig. [[#img-15|15]]  the mesh and step size are not sufficiently fine to reproduce the experimental curves. The results of interest are however clearly met as the identified derivatives follow a similar behavior using both identification techniques, proving that the error is governed by the fluid solution and not by the coupling error. This allows in practice to use a much coarser mesh than  that normally be needed to eliminate the errors in the fluid solution, which in turn provides an advantage for the quick assessment of the performance of the coupling procedure.

We remark that the use of a second order solver for the fluid solution (on the same computational setup) reproduces accurately the experimental curves using the free vibration procedure (see Fig. [[#img-16|16]]). This confirms the good performance of the coupling procedure used.

==8 Conclusions==

A simple model problem has been presented and used to predict the behaviour of some partitioned coupling algorithms for FSI problems. In particular a fractional step type algorithm and a coupling strategy based on the modification of the dynamic equilibrium equations have been analysed. Although a simple one, the model problems allows us to justify the stability or instability of the coupled strategies analyzed and to justify the success of some accelerators proposed in the literature. The effectiveness of the different algorithms  is proved first by comparing their performance with results from the literature and later by showing their effectiveness for the solution of a problem of bridge aeroelasticity. This latter case is treated using free and forced vibration approaches, so as to assess the error induced by the FSI coupling. The results confirm that the coupling error is negligible when compared to the errors in the fluid solution.

==Acknowledgements==

This work was partially supported under the auspices of the Beatriu de Pinos Programme of the Generalitat de Catalunya and by the SEDUREC project of the Consolider Programme of the Ministerio de Educacion y Ciencia of Spain. The authors thank Dr. Lazzari, Prof. Saetta and all the group of Prof. Vitaliani in Padova for  many helpful discussions on the topics analysed.

==REFERENCES==

<div id="cite-1"></div>
[[#citeF-1|[1]]] U. Küttler, W.A. Wall. (2008) Fixed-point fluidstructure interaction solvers with dynamic relaxation, Volume 43 ''Computational Mechanics'' 61-77 

<div id="cite-2"></div>
[[#citeF-2|[2]]] D.P. Mok. (2001) Partitionierte lungsansze in der strukturdynamik und der fluid-strukturinteraktion. Institut fr Baustatik, Fakult Bauingenieur- und Umweltingenieurwissenschaften

<div id="cite-3"></div>
[[#citeF-3|[3]]] H.G. Matthies, J. Steindorf. (2002) Partitioned strong coupling algorithms for fluid– structure interaction. Technical University Braunschweig

<div id="cite-4"></div>
[[#citeF-4|[4]]] H.G. Matthies, J. Steindorf. (2002) ''Strong Coupling Methods''. Technical University Braunschweig

<div id="cite-5"></div>
[[#citeF-5|[5]]] S.Piperno, C.Farhat. (1995) Partitioned Procedures for the transient solution of coupled aeroelastic problems - Part2 - energy transfer analysis and three dimensional applications, Volume 190. ''Computer Methods in Applied Mechanics and Engineering'' 3147-3170

<div id="cite-6"></div>
[[#citeF-6|[6]]] S.Piperno, C.Farhat, B. Larroutorou. (1995) Partitioned Procedures for the transient solution of coupled aeroelastic problems - Part1 - Model problem, theory and two dimensional applications, Volume 124. ''Computer Methods in Applied Mechanics and Engineering'' 79-112

<div id="cite-7"></div>
[[#citeF-7|[7]]] C. Farhat, P. Geuzaine, G. Brown. (2003) Application of a three-field non-linear fluid-structure formulation to the prediction of the aeroelastic parameters of an F-16 fighter, Volume 32. ''Computers and Fluids'' 3-29

<div id="cite-8"></div>
[[#citeF-8|[8]]] B. Hubner, E. Walhorn, D. Dinkler. (2004) A monolithic approach to Fluid-Structure interaction using space-time finite elements. CMAME

<div id="cite-9"></div>
[[#citeF-9|[9]]] E. Walhorn, A. Kolke, B. Hubner and D. Dinkler. (2005) Fluid-Structure Coupling within a monolithic model involving free surface flows, Volume submitted. ''computers and structures'' 2100-2111

<div id="cite-10"></div>
[[#citeF-10|[10]]] A. Halfmann,E. Rank, M. Gluck and others. (2000) A partitioned solution approach for the fluid-structure interaction of wind and thin walled structures. TU Munchen, Universitat Erlangen-Nurnberg

<div id="cite-11"></div>
[[#citeF-11|[11]]] E. Rank, D. Scholz, A.Halfmann and others. (2005) Fluid-Structure interaction in Civil Engineering. Second MIT conference on Computational Fluid And Solid Mechanics

<div id="cite-12"></div>
[[#citeF-12|[12]]] Miguel Angel Fernandez, Jean-Frederic Gerbeau, Celine Grandmont. (2007) "A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid". INRIA

<div id="cite-13"></div>
[[#citeF-13|[13]]] R.Rossi. (2005) Ligh-weight structures - Numerical Analysis and Coupling Issues. University of Padova, Italy

<div id="cite-14"></div>
[[#citeF-14|[14]]] M. Lazzari. (2005) Time domain modelling of aeroelastic bridge decks: a comparative study and an application, Volume 62. ''International Journal For Numerical Methods in Engineering''

<div id="cite-15"></div>
[[#citeF-15|[15]]] C. Felippa. (2004) FSI course: chapters 1 - 11 http://www.colorado.edu/engineering/CAS/courses.d/FSI.d/Home.html

<div id="cite-16"></div>
[[#citeF-16|[16]]] T.J.R. Hughes. (2000) The Finite Element Method. Dover

<div id="cite-17"></div>
[[#citeF-17|[17]]] R. Scanlan, J.J. Tomko. (1971) Airfoil and bridge deck flutter derivatives, Volume EM6. ''J. Mechanical Division''

<div id="cite-18"></div>
[[#citeF-18|[18]]] A. Larsen. (1998) Advances in aeroelastic analyses of suspension and cable stayed bridges, Volume 74. ''Journal of Wind Engineering and Industrial Aerodynamics'' 73-90

<div id="cite-19"></div>
[[#citeF-19|[19]]] E. Briand, S.Piperno. (2002) Validacion du code NSI3FS sur des ecoulements turbulents autour d'un tablier de pont elementaire. INRIA

<div id="cite-20"></div>
[[#citeF-20|[20]]] J.B.Frandsen. (2001) "Simultaneous pressures and Accelerations measured full-scale on the great belt east suspension bridge", Volume 89. Journal of Wind Engineering and Industrial Aerodynamics 95-129

<div id="cite-21"></div>
[[#citeF-21|[21]]] R.P.Selvam, S.Govindaswamy. (2001) Aeroelastic Analysis of Bridge Girder Section using computer modeling. University Of Arkansas internet

<div id="cite-22"></div>
[[#citeF-22|[22]]] S.Piperno, E.Bournet. (1999) Numerical Simulation of wind effects on flexible civil engineering structures. INRIA

<div id="cite-23"></div>
[[#citeF-23|[23]]] A.Larsen, J.H.Walther. (1998) Discrete Vortex Simulation of flow around five generic bridge deck sections, Volume 77. ''Journal of Wind Engineering and Industrial Aerodynamics'' 591-602

<div id="cite-24"></div>
[[#citeF-24|[24]]] A.Larsen, J.H.Walther. (1997) Aeroelastic Analysis of bridge girder sections based on discrete vortex simulation", Volume 67. ''Journal of Wind Engineering and Industrial Aerodynamics'' 253-265

<div id="cite-25"></div>
[[#citeF-25|[25]]] D. Cobo del Arco and A.C. Aparicio Bengoechea. (1999) An Analysis of Wind Stability.Improvements to the Response of Suspension Bridges. CIMNE

==Appendix A. Computation of aeroelastic derivatives for bridge analysis==

===A.1 Forced motion procedure===

One first approach is based on the application of a prescribed sinusoidal motion to the bridge section which is then subjected to the fluid flow. A description of the testing procedure can be found in <span id='citeF-23'></span>[[#cite-23|[23]]] or, with a similar approach in <span id='citeF-21'></span>[[#cite-21|[21]]]. The “experimental” setup can be expressed as follows:

* fix the horizontal and vertical motion of the cross section
* prescribe a sinusoidal pitching motion (a prescribed rotation) with a known angular frequency
* measure the forces and moments on the section
* perform a (linear) least-square fitting of the measured forces with an equation in the form <math display="inline">x(t) = Asin(\omega t + \phi )</math> both for the moment history and lift history.
* use the fitted values for <math display="inline">A  </math> and <math display="inline">  \phi </math> to identify the coefficients as described below

The identification is performed through a (linear) least square fit of the recorded moment and lift histories using  an expression of the form <math display="inline">Asin(\omega t + \phi )</math>. The simple  trigonometric identity

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>Asin(\omega t + \phi ) = Bsin(\omega t) + Ccos(\omega t)) \, ;     \, A = \sqrt{B^2+C^2} \, , \, \phi = tan^{-1}\left(     \frac{C}{B} \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (80)
|}

together with the auxiliary vectors (N is the number of measured time steps, <math display="inline">\delta t</math> the stepsize)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{v}_s =    \left(         \begin{array}{c}sin(\omega \delta t) \\ ... \\ sin(\omega N\delta t)  \end{array}     \right)\; \mathbf{v}_c =    \left(         \begin{array}{c}cos(\omega \delta t) \\ ... \\ cos(\omega N\delta t)  \end{array}     \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
|}

allows to rewrite the fitting problem as

<span id="eq-82"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M} = B_m\mathbf{v}_s + C_m \mathbf{v}_c \, ; \, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (82)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{L} = B_l\mathbf{v}_s + C_l \mathbf{v}_c \, ; \, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (83)
|}

where the moment and lift histories are written in vector form as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M} =  \left(         \begin{array}{c}M( \delta t) \\ ... \\ M( N\delta t)  \end{array}     \right)\; \mathbf{L} =  \left(         \begin{array}{c}L( \delta t) \\ ... \\ L( N\delta t)  \end{array}     \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
|}

by multiplying each of them once by <math display="inline">\mathbf{v}_s^T</math> and then by <math display="inline">\mathbf{v}_c^T</math> we obtain the systems

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{pmatrix}\mathbf{\mathbf{v_s}} \bullet \mathbf{\mathbf{v_s}}  &  \mathbf{\mathbf{v_s}} \bullet \mathbf{\mathbf{v_c}}  \\          \mathbf{\mathbf{v_s}} \bullet \mathbf{\mathbf{v_c}}  &  \mathbf{\mathbf{v_c}} \bullet \mathbf{\mathbf{v_c}}      \end{pmatrix}     \left(\begin{array}{c}B_m \\ C_m \end{array} \right)=     \left(\begin{array}{c}\mathbf{\mathbf{\mathbf{v}_s}} \bullet \mathbf{\mathbf{M}}  \\  \mathbf{\mathbf{\mathbf{v}_c}} \bullet \mathbf{\mathbf{M}}      \end{array}  \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{pmatrix}\mathbf{\mathbf{v_s}} \bullet \mathbf{\mathbf{v_s}}  &  \mathbf{\mathbf{v_s}} \bullet \mathbf{\mathbf{v_c}}  \\          \mathbf{\mathbf{v_s}} \bullet \mathbf{\mathbf{v_c}}  &  \mathbf{\mathbf{v_c}} \bullet \mathbf{\mathbf{v_c}}      \end{pmatrix}     \left(\begin{array}{c}B_l \\ C_l \end{array} \right)=     \left(\begin{array}{c}\mathbf{\mathbf{\mathbf{v}_s}} \bullet \mathbf{\mathbf{L}}  \\  \mathbf{\mathbf{\mathbf{v}_c}} \bullet \mathbf{\mathbf{L}}      \end{array}  \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
|}

which can be solved for the desired fitting coefficients. By considering further that the measured load histories were obtained by prescribing a sinusoidal motion (pitching) of the type

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha = \alpha _0 sin(\omega t) </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
|-
| style="text-align: center;" | <math>     \dot \alpha = \alpha _0 \omega  sin(\omega t) </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (88)
|-
| style="text-align: center;" | <math>     \ddot \alpha = - \omega ^2 \alpha  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
|}
|}

and substituting this in Eqn([[#eq-78|78]]) we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{L} = B_l \mathbf{v}_s + C_l \mathbf{v}_c = (\rho      B^3\omega ^2H_3\alpha _0)\mathbf{v}_s+ (\rho B^3\omega ^2H_2\alpha _0)\mathbf{v}_c </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M} = B_m \mathbf{v}_s + C_m \mathbf{v}_c = (\rho      B^4\omega ^2A_3\alpha _0)\mathbf{v}_s+ (\rho B^4\omega ^2A_2\alpha _0)\mathbf{v}_c </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
|}

which links the fitting coefficients to the aeroelastic derivatives. By proceeding in a totally analogous way prescribing a displacement history in the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>y = y_0 sin(\omega t) </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (92)
|-
| style="text-align: center;" | <math>     \dot y = y_0 \omega  sin(\omega t) </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (93)
|-
| style="text-align: center;" | <math>     \ddot y = - \omega ^2 y </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (94)
|}
|}


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>'''Table. 3''' Relations for computing the aeroelastic derivatives:     the upper index indicates the prescribed sectional motion, the     lower index     the equation used in fitting the load history
|-
|
* Prescribed heave:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>             H_4 = {\frac{{             B^{heave}_{lift}}}{\rho \,{B}^{2}{w}^{2}{ y_0}}} \, ;             \, H_1 = {\frac{{             C^{heave}_{lift}}}{\rho \,{B}^{2}{w}^{2}{ y_0}}}         </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>             A_4 = {\frac{{             B^{heave}_{mom}}}{\rho \,{B}^{3}{w}^{2}{ y_0}}} \, ;             \, A_1 = {\frac{{             C^{heave}_{mom}}}{\rho \,{B}^{3}{w}^{2}{ y_0}}}         </math>
|}
|}
* Prescribed rotation: 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>             H_3 = {\frac{{             B^{pitch}_{lift}}}{\rho \,{B}^{3}{w}^{2}{ \alpha _0}}}             \, ; \, H_2 = {\frac{{             C^{pitch}_{lift}}}{\rho \,{B}^{3}{w}^{2}{ \alpha _0}}}         </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>             A_3 = {\frac{{             B^{pitch}_{mom}}}{\rho \,{B}^{4}{w}^{2}{ \alpha _0}}}             \, ; \, A_2 = {\frac{{             C^{pitch}_{mom}}}{\rho \,{B}^{4}{w}^{2}{ \alpha _0}}}         </math>
|}
|}
|}

which represents a heaving motion, it is possible to identify the remaining four derivatives. Table ([[#table-3|3]]) summarizes  the relevant formulae for computing the aeroelastic derivatives.

===A.2 Free vibration method===

The second procedure for the determination of the aeroelastic coefficients is based on the analysis of the displacement history of the sectional model mounted on a set of springs of known stiffness. This testing procedure is more challenging for the interaction as it involves a completely coupled analysis, which introduces a further error source. An excellent example of application to bridges can be found in <math display="inline"><span id='citeF-19'></span><span id='citeF-22'></span>[[#cite-19|[19,22]]]</math>. In these works however the free vibration procedure is used to assess directly the flutter speed by evaluating the energy variation.

The original identification procedure as proposed by Scanlan <span id='citeF-17'></span>[[#cite-17|[17]]], is based on a three step analysis in which the cross section mounted on springs is first restrained in heave, then in rotation and finally left completely free. This was needed as the available setup allowed to record exclusively the displacement history but not the restrain forces. The approach described in the following is slightly different and it involves only two sets of simulations by making use of both displacement histories and forces records.

The dynamical equation of motion in rotation for a bridge section subjected to small displacements takes the form

<span id="eq-95"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I \ddot \alpha + C_\alpha \dot \alpha +     K_\alpha \alpha = M(t) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (95)
|}

In order to identify the first set of aeroelastic derivatives it is convenient to restrain the motion in heave while leaving  the rotation free. Under this assumption, substituting Eqn( [[#eq-78|78]]) into Eqn( [[#eq-95|95]]) gives

<span id="eq-96"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I\ddot \alpha + C_\alpha \dot \alpha + K_\alpha  \alpha = \rho{B}^{3}w \left(B{ A_2}{ \dot \alpha }+wB{ A_3}\alpha  \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (96)
|}

which provides a description of the motion of the section in its interaction with the fluid. Eqn([[#eq-96|96]])can be formally rewritten as

<span id="eq-97"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I \ddot \alpha + 2 I \omega \xi \dot \alpha +     I\omega ^2 \alpha = 0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (97)
|}

with

<span id="eq-98"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>2 I \omega \xi = C_\alpha - \rho \omega B^4 A_2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (98)
|}

and

<span id="eq-99"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>I\omega ^2 = K_\alpha - \rho \omega ^2 B^4 A_3 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (99)
|}

The differential equations in the form described in Eqn([[#eq-97|97]]) have an analytical solution of the type.

<span id="eq-100"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _{teo}(t) = A e^{-\xi \omega t} sin(\omega      t + \phi ) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (100)
|}

Assuming that the mathematical model proposed by Scanlan describes correctly the real structural behaviour, it is possible to fit the “experimental” solution with a mathematical expression given by Eqn([[#eq-100|100]]). A non-linear fitting of the motion history allows therefore computing the coefficients <math display="inline"> A,\xi ,\omega , \phi </math>. which can be used for the determination of <math display="inline">A_2, A_3</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>A_2 = - \frac{-C_\alpha + 2I\xi \omega }{\rho B^4\omega } \, ;     \, A_3 = - \frac{-K_\alpha + I\omega ^2}{\rho B^4\omega ^2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (101)
|}

The derivatives <math display="inline">H_2</math> and <math display="inline">H_3</math> can be identified making use of the recorded lift forces using a technique that is very similar to the one performed in the previous section. The lift history can  be seen as generated from an imposed displacement in pitch described by an equation such as Eqn([[#eq-97|97]]). The first step is therefore to fit the lift history with the following function

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L_{fitted} = B^{pitch}_{lift}e^{-\xi \omega t}sin(\omega t +     \phi ) + C^{pitch}_{lift}e^{-\xi \omega t}cos(\omega t + \phi ) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (102)
|}

The second step is to substitute Eqn([[#eq-100|100]]) into the moments equation to obtain a theoretical expression in the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L_{teo} = -{\rho {B}^{3}{\omega }^{2}{  A} {e^{-\xi \omega t }}     \left({  H_2} \xi \sin \left(\omega t +\phi \right)-{  H_2}     \cos \left(\omega t +\phi \right)-{  H3} \sin \left(\omega      t +\phi \right) \right)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (103)
|}

by equating the last equations, collecting the “cos” and “sin” terms and equating them to zero, we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>H_2 = \frac{C^{pitch}_{lift}}{\rho B^3 \omega ^2} \, ; \, H_3     = \frac{B^{pitch}_{lift}+C^{pitch}_{lift}\xi }{\rho B^3     \omega ^2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (104)
|}

which is the desired result.

===A.3 Free heave===

The identification of the remaining four derivatives follows a similar path. Once the non linear fitting of the lift equation is performed, two of the aeroelastic parameters can be identified as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>H_1 = - \frac{-C_y + 2m\xi \omega }{\rho B^2\omega } \, ; \,     H_4 = - \frac{-K_y + m\omega ^2}{\rho B^2\omega ^2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (105)
|}

The remaining two parameters, after fitting the moment equation as described above, take the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>A_1 = \frac{B^{heave}_{mom} + C^{heave}_{mom}\xi }{\rho B^3     \omega ^2} \, ; \, A_4 = \frac{C^{pitch}_{lift}}{\rho B^3     \omega ^2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (106)
|}

The identification of the parameters can be performed by non linear fitting techniques. The nonlinear fitting routines provided in scipy a scientific library included in python (see www.python.org) were used in the computations. The convergence of the fitting process was ensured by taking as initial value for the fitting process an estimate for all the relevant parameters.