J Serrano

P Álvarez

J Lewis

INTER-LABORATORY SYNCHRONIZATION FOR THE CNGS PROJECT

CERN will start sending a neutrino beam to Gran Sasso National Laboratory in Italy in July 2006. This beam will cover a distance of around 730 km through the crust of the earth from an extraction line in CERN's SPS to dedicated detectors in Gran Sasso. This paper describes the technological choices made to fulfill the specification of inter-laboratory synchronization in the region of 100 ns, as well as some preliminary results. The common time standard is UTC as disseminated by the GPS system, and the techniques are similar to those used by national metrology laboratories for the manufacturing of UTC itself. In addition, real-time messages sent through the Internet allow the detectors in Gran Sasso to go into calibration mode when no beam is being sent. Data concerning the delay and determinism of this international network link is also presented.

INTRODUCTION

The CNGS (CERN Neutrinos to Gran Sasso) project aims at delivering a high intensity neutrino beam though the Earth's crust from CERN (Geneva, Switzerland) to the Gran Sasso National Laboratory (LNGS) located 120 km east of Rome in Italy. This beam will contain exclusively neutrinos of the muon type, and the detectors in LNGS will try to find tau-neutrinos resulting from the oscillation of muon-neutrinos as they travel though the 730 km straight line separating CERN and LNGS.

In order to correlate the events observed in LNGS with the beam pulses sent from CERN's SPS, both of them will be time-tagged using Universal Coordinated Time (UTC) as a reference, therefore allowing discrimination of events caused by CERN neutrinos with respect to spurious events. The local time bases in LNGS and CERN are provided by GPS Disciplined Oscillators (GPSDO), delivering Pulse-Per-Second (PPS) outputs and 10 MHz clocks to break each second in 10 million ticks. Both GPSDO devices contain Rubidium oscillators to ensure good short-term stability and rely on the numerous Cesium clocks driving the GPS constellation of satellites for their long-term stability. As a result, the quality of the 10 MHz output is excellent, with typical Allan variances [START_REF] Levine | Introduction to time and frequency metrology[END_REF] of 10 -11 for an averaging time of 1 second and 10 -12 for 100 seconds. The limiting factor to achieve good inter-lab synchronization is therefore the variance between the output times of the PPS pulses in LNGS and at CERN. This report presents the results measurement campaigns performed to quantify the systematic offset between the two PPS pulses, used to calibrate the system, and also to study the variance of these measurements to ascertain whether the specified performance (100 ns RMS) is achievable with the basic GPSDO setup, knowing that more involved techniques, such as common view and two-way satellite time transfer [START_REF] Levine | Time and frequency distribution using satellites[END_REF] could be implemented to improve performance if necessary.

Once each laboratory has a facility with a well controlled and characterized GPSDO, it is necessary to take that time reference from there to the place of interest, i.e. the SPS extraction point at CERN and the Opera detector in LNGS. In both cases, the fibre delay incurred has to be measured and taken into account when generating the time tags. The methods used are also presented in this report, along with results from actual measurements.

Another request from the CNGS project concerns the possibility of running calibration procedures in the Opera detector during periods of no beam from CERN. In order to synchronize these calibration runs, taking into account that the SPS is a multi-cycling machine with a multitude of operating modes decided online by operators, a failsafe scheme using UDP packets though the Internet was put in place. We describe the scheme and comment on some preliminary results which confirm that this cheap, non-deterministic solution is indeed good enough for our purposes.

GPSDO CALIBRATION

A GPSDO is a special kind of GPS receiver designed to provide very accurate and stable timing signals. In order to do so, it works to find a very precise position fix during the first hours after power-up. This implies of course that the device should not move at all during normal operation. After many position averages (around 100000 in our case), the device goes into "timing" mode, i.e. it fixes its position and only solves the time equation. This means that a GPSDO can keep track of time with only one GPS satellite in sight (and not three as would be needed to find its position as well) but the more satellites the GPSDO sees the more accurate the timing will be.

The biggest source of errors in the generation of a stable PPS output in a GPSDO is the variable propagation time of GPS signals through their one-way trip from the satellites to the GPSDO, especially through the ionosphere and troposphere [START_REF] Levine | Time and frequency distribution using satellites[END_REF]. For our two relatively close locations, we assume there will be a "common mode" effect whereby the deviations from UTC of our two GPSDOs will go in the same sense and will have roughly the same values. For more local variations, typically of a faster nature, the Rubidium oscillator on each GPSDO will act as a low-pass filter, with a time constant of around 3 hours, smoothing out the short-term variations.

To validate our assumption and gain some extra confidence, we decided to calibrate CERN's GPSDO both against the official source of UTC time in Switzerland (the Swiss Federal Office of Metrology, METAS, in Bern)

THPCH125

Proceedings of EPAC 2006, Edinburgh, Scotland 3092 07 Accelerator Technology T04 Accelerator/Storage Ring Control Systems and against the LNGS GPSDO in Gran Sasso. UTC is in fact a "paper clock" built a posteriori in the "Bureau International des Poids et Mesures" (BIPM) in Paris. Measurements of GPS signals performed using Cesium clocks all around the world are taken by BIPM, properly weighted and averaged, to come up with a table of corrections each given lab should apply to their measurements. In METAS, the time metrology team also generates a "real-time" UTC, called UTC(CH.R), whose typical RMS difference with respect to the paper clock UTC(CH) is 1.3 ns and can therefore be neglected at the scale of 100 ns of interest to us. The result for the measured time offset PPS(GPSDO)-UTC(CH.R) was (+100±50) ns, with the plus sign meaning that the PPS signal coming out of the GPSDO was ahead in time of the one coming from the UTC(CH.R) realization stand. The confidence level given by METAS is such that if we correct for the 100 ns systematic offset, we can trust that the "perfect" UTC PPS pulse would fall in the ±50 ns window around our PPS with a probability of 95%.

We could have requested a similar calibration procedure to be performed on the LNGS GPSDO, but it was felt that the ultimate test would actually consist in measuring the two GPSDOs against one another. The result of these measurements, performed over a 12 day period in March 2006 in the LNGS site, can be seen in Figure 1. The average value of the data is 353 ns, meaning that the LNGS GPDO is really well ahead of UTC. The 100 ns of offset between the CERN GPSDO and UTC had already been compensated using a programmable register in the device. The excursions in the waveform are well contained within the 100 ns envelope allowed by the specifications of the CNGS project. A certain periodicity can be seen in the waveform with a period of one day. This can be explained by the fact that the GPS receiver antenna sees the same satellite configuration in the sky every 24 hours, so any systematic effects within each GPSDO will reproduce with that same period.

Time difference GPSDO(LNGS)-GPSDO(CERN)

FIBRE DELAY CALIBRATION

The GPSDO installed at CERN is located in the CERN Control Centre (CCC), while the LNGS system runs in the network routers room in the laboratory site. In both cases these systems are far away from the places were the precise timing signals are needed (point 4 of the SPS at CERN and the Opera detector in LNGS). In both labs, dedicated electronics take the PPS and 10 MHz outputs of the GPSDO and use them to encode timing information on a fibre link. At the receiving side, another set of electronic cards uses this information to generate pulses or to time-tag external events. Since the moment of arrival of the messages is as important as the message content itself, it is very important to calibrate the delay induced by the fibre in the timing transmission.

A very simple scheme requiring only one additional fibre between the two points of interest has been used at CERN to calibrate the delay between the CCC and SPS point 4 (building BB4). One optical transmitter and one optical receiver are used in a first phase to transmit the PPS pulse through the additional link from the CCC to BB4. The PPS pulses at the two ends are then time tagged using standard CERN timing receivers. Then we take the same transmitter module to BB4 and the same receiver module to the CCC and repeat the operation using the PPS produced by the timing receiver in BB4. One can prove that the time offset between the CCC and BB4 is given by: 2 4 4 4

FBB RBB FCCC RCCC CCC BB t t t t T - - + = Δ - (1)
Where t RCCC is the time tag of the "return" (BB4 to CCC) pulse as seen in the CCC, t FCCC is the time tag of the "forward" (CCC to BB4) pulse as seen in the CCC and so on. In the above formula, it has been assumed that the delay of the fibre link is the same in both directions, and also that it stayed stable during the 1 hour it took to take the material from the CCC to BB4 and vice-versa. This assumption has been validated by leaving the system running in the "return" configuration for some weeks, i.e. emitting the PPS pulse in BB4 and looking at only one half of the equation. Any fibre delay variation induced by temperature changes should be visible in this setup. The observed variations are in the 1-2 ns region and can be neglected for our purposes. The final result for the time offset between the CCC and BB4 as computed using this scheme is 10.584 V FRUUHVSRQGLQJ URXJKO\ WR NP RI fibre with a delay of 5 ns/m. More elaborate systems can be used to calibrate the fibre delay online. For example, we could have two calibrated sets consisting of a transmitter and a receiver working all the time and logging these time differences to correct the time tags offline. Online correction of the delays, which would involve active devices such as phase shifters or fibre stretchers, is not needed in our application.

SOFTWARE SYNCHRONIZATION WITH UDP PACKETS

The Opera detector in LNGS needs short calibration runs for several electronic subsystems as frequently as possible. It was agreed to have a fail-safe mechanism to inform Opera through the Internet that it can run the calibration procedure safely whenever CERN is not sending any beam to LNGS. By fail-safe we mean that the loss of packets in the network should never result in loss of useful data acquisition.

The SPS is a multi-cycling synchrotron with several users, including LNGS, the North Experimental Area and LHC. A simplified view of the multi-cycling mechanism would include different cycles meant for different destinations grouped in an entity called the super-cycle, which is run cyclically until operators decide to introduce changes. The real situation is a bit more complicated, with interlock bits driving asynchronous jumps from one normal cycle to a spare cycle online. In any event, all cycle changes occur at boundaries of the 1.2 second basic heart beat that drives all CERN accelerators. The proposed solution for the Opera calibration runs consists in sending a structure through the network at the beginning of every basic period to inform Opera of what's going on in the SPS. Typically, this structure contains the UTC time for the start of the basic period, the UTC time at which the packet was sent (for sanity checks on the receiving side), the number of basic periods for the current cycle and the current destination for the beam.

In Opera, a task can listen to the UDP packets coming from CERN and find out for example that the coming 9 basic periods (i.e. 10.8 seconds) will be used in the SPS to generate beam for the North Area, so it is safe to go ahead and run the calibration, switching automatically back to data taking after this time is elapsed.

In order to validate the concept, we tested the round trip delay of UDP packets between CERN and LNGS. As can be seen in figure 2, none of the packets makes it in less than 25 ms, and from then on the distribution presents a maximum at 26 ms and then decreases getting to negligible levels after the 50 ms region. Of the 14400 packets sent for the test, 2 were lost and 1 actually showed a round trip delay of a full second, showing that the implementation of a fail-safe mechanism over this physical layer is absolutely necessary.

Another important aspect to be taken into consideration is how to get the current UTC time within the receiving task for sanity checking purposes. The most common approach is to issue a UNIX system call that relies on the presence of a Network Time Protocol (NTP) server nearby. An NTP server typically contains a GPS receiver and an Ethernet board, and provides time-of-day services to a set of computers in the network. If a network node is very far away (in terms of network connections) from its closest NTP server, there is a risk that the network packets coming from the NTP server will take a lot of time to reach that node, and insane results will appear within our sanity check routine with no reason. The CERN site is well stocked with NTP servers, but the same is not true in LNGS. The optimal solution is to equip the GPSDO used for PPS and 10 MHz outputs with a special add-on card implementing NTP server functionality. In this case, there is no risk of incoherent results and there is the additional bonus of physical proximity between the NTP server and its client.

Figure 1 :

 1 Figure 1: time difference measurements between the PPS outputs of the LNGS and CERN GPSDOs.

FrequencyFigure 2 :

 2 Figure 2: Histogram showing typical round trip delay times for UDP packets between CERN and LNGS.

ACKNOWLEDGEMENTS

We would like to thank Laurent-Guy Bernier and Gregor Dudle of the Swiss Federal Office of Metrology (METAS) for their help in GPSDO calibration as well as for giving us general insight in time transfer issues. Many thanks as well to Konrad Elsener for his support within the CNGS project.


Catherine Ramirez
email: catherine.ramirez@ec-nantes.fr

Alban Leroyer
email: alban.leroyer@ec-nantes.fr

Michel Visonneau
email: michel.visonneau@ec-nantes.fr

Yann Roux

Ramirez V Catherine

Study of asynchronicity in the fluid-structure interaction domain

Keywords: FSI, asynchronous algorithm, numerical stabilization, strong couplings, RANS CFD, non-linear FEM, complex structures, Schwarz methods, implicit internal algorithm, Jacobian of the interface

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

When studying fluid-structure interaction (FSI), numerical instabilities are associated with the use of a partitioned approach [START_REF] Park | Stabilization of staggered solution procedures for fluid-structure interaction analysis[END_REF]. In order to guarantee robustness and accuracy, numerical methods that address applications dealing with incompressible viscous flows primarily use implicit algorithms; these being predominantly sequential with exchanges between domains synchronized.

Given the high computational cost associated with accurate Computational Fluid Dynamics (CFD) models, the existing procedures seek to reduce the global computational time, mainly by reducing the computational time devoted to the resolution of the structure. This is why it is common to simplify 3D structures with reduced-dimension models, such as beams, membranes, or multi-body systems, as is the case in [START_REF] Nayer | Interaction Fluide-Structure pour les corps élancés[END_REF], [START_REF] Durand | Interaction fluide-structure souple et legere, applications aux voiliers[END_REF], [START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF] and [START_REF] Yvin | Interaction fluide-structure pour des configurations multi-corps. Applications aux liaisons complexes, lois de commande d'actionneur et systèmes souples dans le domaine maritime[END_REF]. However, to address applications that exhibit complex geometries or complex structural behavior, these simplifications are no applicable. To avoid facing prohibitive computational times, new algorithms must be developed.

Motivated by this fact, the present work explores a strategy that allows space for complex structural models without penalizing global computational time. The proposed methodology seeks to go beyond standard exchanges governed by classical coupling schemes and to study an asynchronous approach in the FSI domain.

STATE OF THE ART

When carrying out a partitioned solution of a coupled problem, Lions [START_REF] Lions | On the schwarz alternating method. i[END_REF] and Farhat [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] show that an appropriate distribution of computational resources has to be defined for it to be efficient. However, it is not an obvious task to predict in advance the computational load of each domain. Even using techniques to dynamically balance the distribution in run-time, the difference presented in the computational loads is the reason why, unfortunately, idle times are often introduced when solving the coupled subdomains in parallel. Taking advantage of this issue has prompted investigation of an asynchronous approach in some domains.

For example, the work of Magoules et al. [START_REF] Magoules | Asynchronous optimized schwarz methods with and without overlap[END_REF] reports the achievement of a significant reduction in computational time when asynchronicity is introduced in the domain of acoustics. In this work applications described by the Helmholtz equation and the Laplacian operator are solved using an asynchronous version of the optimized Schwarz's methods, with and without overlap, and considering Robin and mixed boundary conditions at the interface. Evidence of convergence for the methods is provided, and the gain in computational time found when using an asynchronous algorithm instead of the synchronous version is in the range of 30 -50%.

Other references of successful studies using asynchronous approaches in the frame of Schwarz (classic) methods are the works of [START_REF] Frommer | Asynchronous weighted additive schwarz methods[END_REF], [START_REF] Laitinen | Asynchronous domain decomposition methods for continuous casting problem[END_REF] and [START_REF] Chau | Asynchronous schwarz methods applied to constrained mechanical structures in grid environment[END_REF]. The work of Frommer et al. [START_REF] Frommer | Asynchronous weighted additive schwarz methods[END_REF] solves a non-singular, linear, second order elliptic, boundary value problem of the form Ax = f . To do this, it uses an asynchronous version of Schwarz's methods, with and without overlapping between domains. With ten million variables in 256 processors, the results obtained in this work show improvements in the computational cost of up to 50%. Similarly, the convergence properties of the method are studied, and convergence is shown to be ensured when the coefficient matrix A is monotone and when A is an H-matrix.

On the other hand, the work of Laitinen et al. [START_REF] Laitinen | Asynchronous domain decomposition methods for continuous casting problem[END_REF] studies the non-linearities of casting modeling using finite difference methods to solve the convection-diffusion equation describing the model. Convergence of the proposed asynchronous method is studied, and an associated geometric criterion of convergence is stated. The work of Chau et al. [START_REF] Chau | Asynchronous schwarz methods applied to constrained mechanical structures in grid environment[END_REF] studies the stationary obstacle problem, using a multivalued problem type formulation. It studies the rate of convergence of the methods used in the context of convex optimization, and states that the properties obtained from discrete problems ensure the convergence of the synchronous and asynchronous algorithms studied.

Showing that asynchronous algorithms can be robust and advantageous to compensate for the disadvantages of the aforementioned load imbalance present in parallel calculations, the results observed in these works motivate the present study of a similar strategy in the field of FSI. Nonetheless, it is emphasized here that all the reviewed cases consider applications where all subdomains solve the same system of equations. Therefore, the application of an asynchronous approach in a coupled solution has not yet been validated. In addition, emphasis is placed on the fact that no reference was found in the literature to a criterion for an appropriate load balance, or a particular subdomain distribution that guarantees that asynchronous versions are always faster.

MATERIALS AND METHODS

Collectively solving all the unknown variables that characterize the fluid and structure domains and implicitly treating the interface conditions, the monolithic approach, presented by Eq. ( 1), guarantees the stability and convergence of an FSI solution.

F C-→ sf C-→ f s S x f x s = S f S s (1)
In this equation, F and S refer to the linearized fluid and structure operators, C-→ f s relates the loading of the structure given the fluid efforts and C-→ sf considers the deformation of the fluid domain due to the structural response. The external sources are represented by the blocks S f and S s , and the unknowns of each domain by x f and x s .

Exploiting the maturity of specific solvers to solve the subdomains that naturally lie in the original domain is the partitioned method's strategy. To ensure robustness and accuracy when addressing applications that deal with incompressible viscous flows, the numerical methods used in this approach primarily use implicit algorithms.

Two implicit algorithms are distinguished in this work. In the implicit external algorithm, Fig. 2, the fluid domain and structure convergence loops are independent and full convergence of each of the domain is required at each FSI iteration. The internal implicit algorithm, Fig. 1, integrates the complete resolution of the structure into the non-linear loop of the fluid's solution. This makes the convergence loop of the fluid domain solution coincide with that of the FSI solution, and the convergence of the fluid domain takes place simultaneously with the FSI convergence [START_REF] Leroyer | Etude du couplage écoulement/mouvement pour des corps solides ou à déformation imposée par résolution des équations de Navier-Stokes: contribution à la modélisation numérique de la cavitation[END_REF]. This algorithm is used within the coupling tool K-FSI to sequentially solve FSI problems. The main drawback of partitioned schemes is an inherent instability caused by the so-called added-mass effect, as shown in [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF], [START_REF] Söding | How to integrate free motions of solids in fluids[END_REF] and [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF]. To provide robust and faster convergence, the implicit external algorithm is often solved using a Quasi-Newton approach, in which a numerical Jacobian expressing the sensitivity of domains to changes in boundary conditions is evaluated. Some examples can be found in [START_REF] Degroote | Stability of a coupling technique for partitioned solvers in fsi applications[END_REF], [START_REF] Lindner | A comparison of various quasi-newton schemes for partitioned fluidstructure interaction[END_REF] and [START_REF] Gatzhammer | Efficient and Flexible Partitioned Simulation of Fluid-Structure Interactions[END_REF].

In a different way, the robustness of the implicit internal algorithm is obtained by using the Jacobian of the fluid that corrects fluid efforts taking into account the response of the structure, as shown in the works of Yvin [START_REF] Yvin | Interaction fluide-structure pour des configurations multi-corps. Applications aux liaisons complexes, lois de commande d'actionneur et systèmes souples dans le domaine maritime[END_REF] and Durand [START_REF] Durand | Interaction fluide-structure souple et legere, applications aux voiliers[END_REF]. Regardless of the application coupling level, the number of iterations required for the convergence of the FSI problem is the same as for solving the fluid domain. Furthermore, the overall computational time can be estimated as being proportional to the fluid iterations, if the time of the structural solver is limited [START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF].

To address problems where complex structural models are required and that, therefore, would significantly impact global computational time, this paper explores a parallel iterative procedure in which each process continuously performs the resolution of a subdomain without waiting for the neighbor solution to be completed. Stability issues are investigated and methods to improve the stability are provided. This work is carried out in the context of the implicit internal algorithm, and therefore, asynchronicity is posed within the non-linear loop. Nonetheless, it would be enriching to study asynchronicity with implicit external algorithms, or even explicit algorithms.

Most likely, the asynchronous approach involves context-specific numerical instabilities, but anyhow, the analysis presented in this work may represent an advance in the general understanding of parallel coupling techniques with synchronous and asynchronous behavior, which could be useful for future studies.

The K-FSI coupling tool

The present work has been implemented in K-FSI, a commercial tool coupling the structural code K-Struct and the fluid solver ISIS-CFD; The latter being part of the FINE TM /Marine suite. K-FSI resolves the Neumann-Dirichlet boundary conditions at the wet surface. The coupling between the two solvers is carried out using a block-LU decomposition of the monolithic system [START_REF] Durand | Interaction fluide-structure souple et legere, applications aux voiliers[END_REF]. This approach, along with the approximation of the first Jacobian of the interface, allows K-FSI to predict the behavior of applications ranging from weakly to strongly coupled.

Writing the lower decomposition of the system given by Eq. (1) the system of equations reads

(S -C f s F -1 C sf ) x s = S s -C f s F -1 S f F x f = S f -C sf x s (2a) (2b)
where the operators are defined in a similar way as for the Eq. (1). Using the iteration written for the second term given by Eq. (3), and manipulating Eq. (2b), as summarized by Eq. (4), the system of equations solved by K-FSI can be written as in Eq. ( 5).

C f s F -1 C sf x i+1 s = C f s F -1 C sf x i s + J(C f s F -1 C sf )(x i+1 s -x i s ) (3)
C f s F -1 S f = C f s x f + C f s F -1 C sf x s (4)    (S -J(C f s F -1 C sf )) x i+1 s = S i s -J(C f s F -1 C sf ) x i s -C f s x i f F x i+1 f = S i f -C sf x i+1 s (5a) (5b)
where the operator J F S = J(C f s F -1 C sf ), stands for the Jacobian matrix of the fluid efforts with respect to the structure variables, or first Jacobian of the interface, which is mainly characterized by inertial effects.

Meanwhile, for a non-linear FSI resolution, a fixed-point iteration of the solved structural system is given by the balance of the forces related to the structure and fluid, F s and F f respectively, and can be expressed as presented by Eq. ( 6)

S x i+1 s -S i s + C f s x i f -J(C f s F -1 C sf )(x i+1 s -x i s ) = 0 (6) F s (x i+1,j s ) F f (x i+1,j s )
or as shown by Eq. ( 7), if the internal loop of the structure is considered.

F s (x i+1,j+1 s ) + F f (x i+1,j+1 s ) = (7)
F s (x i+1,j s ) + F f (x i+1,j s ) + (J(F s ) + J(F f ))(δu j )
where J(F s ) and J(F f ) represent the Jacobian of the structural and fluid systems, and δu j = x i+1,j+1 s -x i+1,j s .

PROPOSED NUMERICAL METHOD

The implementation of the asynchronous algorithm is carried out using the ZeroMQ embeddable networking library [START_REF] Hintjens | ZeroMQ: messaging for many applications[END_REF]. With this library, instead of passing arguments on function calls, the boundary data is shared via the sockets that are available through the library.

In a first stage, a parallel algorithm was implemented. For both solvers to advance simultaneously in time, a prediction of the structural response was implemented within the fluid solver. The prediction, which takes place at each new time step, is first order during the second time step, and second order from the third time step onwards. Requiring to store some good values to be accurate and robust, this prediction was found to be not strong enough to guarantee convergence of solution.

Therefore, in order to improve the prediction a second approach was carried out. It consists in executing in a sequential way the two first FSI iterations at each time step, all along the computation. By doing this, equivalent accuracy as for the sequential algorithm is observed for all considered cases when advancing to a new time step. Tab. 1 presents the nomenclature of the algorithms used for the evaluation of this work. Three parameters are considered to define the pattern of information exchange between the two solvers: the architecture pattern: S/P → Serial/Parallel, the exchange protocol: SYN/ASY → Synchronous/Asynchronous, and the strategy to advance in time: PF/NPF → Advance in time with fluid prediction or not, and PS/NPS → Advance in time with structure prediction or not. The notation 2itSeq stands for two iteratios that are carried using a sequential scheme.

Asynchronous method

The objective of the asynchronous algorithm, depicted in Fig. 3, is to solve each of the domains using the latest available information, either local or updated by its counterpart. Considering that the number of exchanges within the nonlinear loop will depend on the computational load of each simulation, the resulting algorithm can be thought of as naturally lying between the implicit external algorithm and the implicit internal algorithm.

In order to set this type of interaction, a channel of communication is established using sockets, and blocking and non-blocking messages are used according to the stage in the loop. All along the simulation, the sockets remain open for both depositing and retrieving information. Although the exchange between the solvers is synchronized at the beginning of each new time step, all along the non-linear loop, each solver is free to send the results obtained after solving its own domain and free to retrieve the information from its counterpart in case it is already available. In the case that no new information has been placed, the solver that has finished one computation will begin a new solution using the last information received from its counterpart and its own updated information. To better understand the operation of the asynchronous algorithm, the Fig. 3 highlights the 'FSI iterations'. An 'FSI iteration' is the equivalent of an execution of the two solvers, each one using updated boundary conditions. Therefore, it will only be complete after both solvers have computed one solution using a set of information from the neighboring domain that has evolved.

In the case of synchronized operations, given that information is updated at each non-linear iteration, an FSI iteration coincides with a non-linear iteration. However, in the asynchronous algorithm the number of non-linear iterations is increased every time one solution is computed, while the number of FSI iterations will only be increased when a solution is carried out with new information from the neighboring domain.

This work added two new criteria to be used with the proposed method. The first one consists in considering simultaneously the convergence of the structure and the convergence of the fluid. The second one, depicted in Fig. 4, it consists of evaluating the criterion when the solutions are carried out using new boundary conditions, as is the case of synchronized algorithms.

Numerical instabilities

Two new sources of numerical instabilities were identified. The first consists of the additional delay introduced by the parallel operation, and the second is related to the impact of the local convergence of the fluid domain on the global convergence. Both sources are described below.

Additional delay when operating in parallel

To explain the instability caused by the parallel operation, the following analysis is performed. When an FSI solution is computed using a monolithic approach, it can be thought as if both systems were solved simultaneously for the same non-linear iteration i, as depicted in Fig. 5a. By considering the action and the reaction of both domains together, the solution is reached in a single iteration. However, when the solution is split, the same non-linear iteration i is divided in two stages, a lag in the response of the neighboring domain is introduced, and a loop of convergence becomes necessary. On one hand, as shown in Fig. 5b, for the non-linear iteration i of the sequential algorithm the structure is solved to find the unknowns x i using, as boundary conditions, the fluid efforts σ i-1 from the previous non-linear iteration. However, the fluid solver solves for the unknowns σ i using the structural variables of the current iteration x i . On the other hand, for the parallel algorithm, the structure solves for x i using the fluid efforts from the previous non-linear iteration σ i-1 , and, likewise, the fluid solves for σ i using the structural variables from the previous non-linear iteration x i-1 , as well. See Fig. 5c. If the use of the Jacobian of the fluid is considered as the way in which the structural system recovers from the delay in the boundary conditions imposed by the fluid solution, it can be said that, for the parallel algorithm, the delay in the boundary conditions imposed by the structure has yet to be overcome.

it i (itnl i) it i+1 (itnl i) (a). Monolithic approach.

Impact of the local convergence

Contrary to what was highlighted for other domains in section 2, the present work found that when an asynchronous approach is used in the FSI domain, local convergence increases numerical instabilities. One of the reasons is that during the instances that local convergence cycles take place, the nonlinearities of the coupling describing the FSI context are not resolved. The other is that, the behavior of the fluid solution within local convergence cycles is marked by the use of an under-relaxation factor used to stabilize the coupling between pressure and velocity within the non-linear loop.

To illustrate impact of the local convergence, Fig. 6 presents the convergence of the asynchronous algorithm 210 (P-ASY-PF-PS), within a time step, when solving a test case with a medium coupling level (m F /m S = 0.1). Gray circles along with numbers indicate the instances in which new boundary conditions are imposed to the fluid domain. In any of the intermediate instances between these circles the fluid domain is solved using only local updated information.

Figure 6 shows that the cycles of local convergence carried out within the fluid domain increase the amplitude of the oscillations in the fluid response. For stronger couplings there is a larger impact of the cycles of local convergence performed by the fluid solver.

Proposed stabilization method

To overcome possible sources of instability within an asynchronous approach, a three-step stabilization strategy is proposed. To begin, an adjustment is made to the formulation of the operator that corrects the structural system. Secondly, the use of an additional operator is proposed to simultaneously correct the fluid system. And finally, to accompany the isolated fluid convergence cycles, the use of a simplified auxiliary model is proposed and studied.

Fluid jacobian formulation adjusted for parallel operation

To correct the structural system, and thus stabilize the FSI solution, for a given FSI iteration the Jacobian of the fluid, or first jacobian of the interface, must be computed considering the kinematics of the structure that produced the fluid efforts that are being used, as a boundary condition, on the specific FSI iteration.

This implies that for the case of the parallel algorithm the correction should be computed with a kinematics that is one iteration further behind than in the sequential case. For the parallel operation case, the explicit term of the correction becomes

J(C f s F -1 C sf ) x i-1
s , and equations 5a and 5b becomes 8a and 8b, respectively.

   (S -J(C f s F -1 C sf )) x i+1 s = S i s -J(C f s F -1 C sf ) x i-1 s -C f s x i f F x i+1 f = S i f -C sf x i s (8a) (8b)

Second or structural Jacobian of the interface

The work of Durand et al. [START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF] presents and uses an expression for the fluid Jacobian interface operator J F S , that is internally computed within the structural solver. To overcome the instabilities arising from the additional delay introduced by the parallel operation, the present work proposes to simultaneously account for the impact of the fluid response on the structural variables through the use of the structural Jacobian interface J SF , or second Jacobian operator defined at the interface.

To do so, the 'upper decomposition', of the LU decomposition of the monolithic system is used to obtain a new formulation of the Eq. (8b), as presented by Eq. ( 9), that should be used to solve the fluid domain if a parallel operation takes place.

(F -C sf S -1 C f s ) x f = S f -C sf S -1 S s (9)
Similarly, as was done in the work of Veldman et al. [START_REF] Veldman | Accelerated free-surface flow simulations with interactively moving bodies[END_REF] and Veldman et al. [START_REF] Veldman | Preventing the added-mass instability in fluid-solid interaction for offshore applications[END_REF], the second Jacobian J SF = -C sf S -1 C f s , corresponding to the variations in the velocity and position of the nodes of the mesh at the interface with respect to the fluid variables, is here approximated using the modal analysis approach that is described in section 4.3.3. This modal approach is used both to correct the fluid system and as an auxiliary model to accompany the isolated convergence of the fluid in case of asynchronous communication within the solvers.

Modal based correction

To avoid carrying two FSI couplings when operating in parallel, the present work formulates a corrective model for the structure δx, that uses the difference in the fluid efforts between the last two nonlinear iterations δf , projected onto the modal shape φ i , for the held modes n.

The use of the 2nd order Backward Difference time scheme and the manipulation of the equations defining the existing modal model within ISIS-CFD2 , lead to an approximation of the change in the structural behavior given certain change in the fluid efforts, or the second Jacobian of the interface, as it is presented by Eq. [START_REF] Frommer | Asynchronous weighted additive schwarz methods[END_REF].

δq i k+1 = φ T i δf k e 2 c + ω 2 i 1 1 + Ca where δx = n i=1 δq i k+1 φ i (x) (10)
In this equation, ω stands for the natural frequencies and e 2 c corresponds to a coefficient linked to the time scheme. Similar to what is done in the work of Yvin et al. [START_REF] Yvin | Added mass evaluation with a finite-volume solver for applications in fluid-structure interaction problems solved with co-simulation[END_REF], the coefficient 1 /1+ Ca is used to stabilize the pseudo-coupling taking place between the fluid solver and the modal based model. Here Ca = M add/M , and M add corresponds to the added-mass coefficient related to the modal shape φ i which can be computed internally within the ISIS-CFD solver. As shown in Fig. 7, this approximation of the second Jacobian of the interface does not modify the system of equations of the fluid domain. Instead, it produces a correction to the boundary conditions calculated by the structural solver, and can be seen as a sequential coupling between the fluid solver and a modal based model of the structure. The implemented correction can be used both to correct the additional delay present in parallel operation and to work as an auxiliary solver during instances of local fluid convergence.

In Fig. 7, δx i+1 stands for the structural response found by means of the modal approach, and the hat accent ˆidentifies the modified procedures in both, structural and fluid solutions, to stabilize the FSI computation.

Local convergence with auxiliary coupling

Non-linear iterations Fig. 8 illustrates the effect of coupling the fluid solution to the proposed modal based correction during the instances when there is no communication with the main solver. The figure presents the evolution of the convergence within the non-linear loop for the case that performs exchanges with the main structural solver every 12 non-linear iterations, approximately.

In this figure it is possible to recognize a shape of convergence similar to that obtained when using the main solvers and the sequential algorithm. When comparing with the main solver, it can be seen that the auxiliary coupling has much less impact on the evolution of the solution. Given that the proposed modal formulation does not have a history of the deformations over time, but instead is based on the difference of two non-linear iterations, the resulting convergence finds a state of equilibrium different from that of the reference solution.

SIMULATION SETUP

The test case used for validation is a benchmark defined in the work of Turek and Hron [START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF]. It consists of a two-dimensional geometry composed of a fix cylinder and an elastic cantilever behind it, where an incompressible laminar flow around the geometry induces oscillations in the structure. In this work three different FSI configurations whose solution difficulty evolves incrementally are defined. Among the three presented scenarios for FSI, the socalled FSI2 and FSI3 test cases are selected here for the assessment.

A third case FSI2x100 is introduced here. It consists of a modification of the FSI2 test case, where the density of the structure is increased by a factor of 100. In this way by the modification of a single parameter, the physical coupling between the domains is reduced and the numerical solution is simpler. To illustrate the numerical complexity of each case, Fig. 9 compares the displacement of the tip of the structure for the three cases studied.

When using the implicit internal sequential algorithm, the test cases FSI2x100 and FSI2 do not require a strong numerical stabilization, carried out in this work through the Jacobian of the fluid. However, the solution of the FSI2 test case requires the use of an under-relaxation factor to achieve convergence.

In all cases, the fluid domain is discretized using an unstructured hexahedral mesh. While for the FSI2x100 test case it consists of 11206 cells, for the FSI2 and FSI3 test cases it is composed of 107672 cells. The flexible flap is discretized using 21 Euler-Bernoulli beam elements for the two initial cases, and 34 for the FSI3 test case. For all cases, the node connecting the beam and the cylinder is fully constrained. A no-slip boundary condition is used along the entire surface, and the input flow velocity is prescribed using parabolic velocity profile. A time step of 0.005s is used to solve the case FSI2x100, and the cases FSI2 and FSI3 use a time step of 0.001s. The Tab. 2 presents the structural and fluid parameters used to set up the three test cases. Given the simplicity of the structure, the computational time of the fluid is much more important than that of the structure. To simulate the hypothetical case of an application in which the solution of the structure will take longer than that of the fluid, a dummy matrix operation is introduced in the structural solution. To study the impact of additional solutions, with and without the proposed stabilization method, different computational times are defined for the structure using a dummy matrix operation. A frequency of exchanges, FSI frequency, is defined by changing the size of the matrix The FSI iterations are the criteria used to compare the results among these hypothetical scenarios.

NUMERICAL RESULTS

2D cylinder flap: FSI2x100

The test cases established through the fictitious operation are presented in Tab. 3. The frequency of FSI iterations, or FSI frequency is defined as the average number of fluid solutions for each FSI solution, within a non-linear loop. Figure 10 shows that even for a weakly coupled application the behavior caused by the use of the pressurevelocity under-relaxation factor has a negative impact on the local convergence. It is also observed that the accumulation of small errors leads to the difference in the converged values.

Figure 11 compares the fluid effort obtained with the implicit, synchronous and asynchronous, algorithms. In general, regardless of the FSI frequency, a global agreement is found for all the asynchronous tests both the fluid efforts and the structural displacements. The result obtained with an explicit algorithm is included in Fig. 11a to make clear that, although it is a weakly coupled application, it requires a minimum of interactions between the domains within the non-linear loop to obtain correct results.

Figure 11b allows noticing that slight oscillations are introduced for some of the asynchronous cases. Two reasons explain this, the first is that FSI iterations are reduced, and the second is the impact of local convergence.

2D cylinder flap: FSI2

The test cases established through the fictitious operation are presented in Tab. [START_REF] Debrabandere | A staggered method using a modal approach for fluid-structure interaction computation[END_REF]. In all the results presented, the modal based correction is used to anticipate the structural position. For an initial time step, Fig. 12 shows the convergence of the fluid efforts of the different cases, obtained without using the modal correction as an auxiliary solver. From Figure 12, it can be seen that the magnitude of the efforts increases mainly during the initial iterations. After this the amplitude of the oscillations remains of the same order. In none of the cases the local convergence of the fluid contributes to the global convergence.

Modal-based auxiliary coupling for stabilization

Using the modal based correction during the periods of the local convergence of the fluid domain has an impact on the computational time of the fluid domain. Table 5 presents the average of the FSI frequency and of the FSI iterations reached at each time step for each of the studied cases. Knowing that good results could be obtained for this case using 20 nonlinear iterations with the 201 parallel algorithm, the FSI iterations for these tests were kept close to a value of 20 so that the results in this section are comparable.

The impact of the numerical parameters of the correction when using the synchronous parallel algorithms is presented below. 100 with m=NA, r=0.5, it fsi=30(f=1), or=4 211 with m=NA, r=0.5, it fsi=8.5(f=5.9), or=4 211 with m=NA, r=0.5, it fsi=7(f=7.1), or=4

(c). FSI frequencies of 5.9 and 7.1 6.2.1.1. Amount of the correction through the prediction coefficient a.

Fy[N]

Non-linear iterations According to what was presented in section 4.3.3, the corrected fluid effort σ can be defined as the linear combination between the boundary conditions given by the structural solver x i and those obtained with the modal analysis approach δx. Therefore, the amount of correction considered can be controlled by means of a coefficient a, as σi+1 (x i + a • δx i+1 ), where the coefficient a defines the amount of computed modal correction that is applied. Figure 13 presents the results obtained with the parallel algorithm 200, using a coefficient a of 0.5, 0.7 and 1.0. Only the first natural mode of the structure is used to compute the prediction. For different instances of the computation it is seen that the larger the parameter a the lower is the amplitude of the oscillations in the response.

Number of modes.

To study the influence of the number of modes computations using 1, 2, 3 and 5 natural modes were run. To isolate the effect of the modes, the under-relaxation of forces is not used and the coefficient a is set to 1.

With this setup, it was found that none of the parallel algorithms was able to converge using only one mode. Nevertheless, they converged when two or more modes were used. Figure 14 compares the convergence, in a middle instance, for all three synchronous algorithms. It is observed that the more modes are considered the lower the amplitudes of the oscillations. However, the contribution after two modes, at least for this case, is not significant. Although the same convergence criterion is used for all algorithms, the sequential algorithm twice as fast. For this reason, for the same number of non-linear iterations, a difference in the values is observed.

Under-relaxation of forces combined with the modal correction.

Although it was seen that the problem FSI2 can be solved without using the relaxation of forces when using more than 1 mode for the modal correction, it is also possible to achieve convergence by combining a single mode and under-relaxation of forces.

Fy[N]

Non-linear iterations Nonetheless, comparing Fig. 15a and Fig. 15b, it can be noticed that if more than 1 mode is used, the relaxation of the forces generates an inverse effect decreasing the quality of the response. This is to be expected if it is considered that the relaxation of forces is carried out within the structure code, which means that the effort used to calculate the modal correction would be different from that used to calculate the structure. 6.2.1.4. Results on the asynchronous approach.

For an early instance, the Figure 16 summarizes the convergence obtained with the asynchronous algorithm 211 when using the auxiliary coupling, for the different scenarios presented in the Tab. [START_REF] Degroote | Stability of a coupling technique for partitioned solvers in fsi applications[END_REF].

Figure 16a shows that the use of a complementary structural model during the local convergence of the fluid improves the convergence of the asynchronous algorithm. However, the amplitude of oscillations is still too large and therefore none of the cases achieve a proper convergence. Curiously, similar as it was seen in the case FSI2X100, the scenario with the lower FSI frequency, and therefore performing exchanges more frequently, does not present the best performance. If the results are also compared with those of the parallel synchronous algorithm 201, which uses the modal correction to predict the structural response as shown in Fig. 16b, it can be observed that even the best results obtained with the asynchronous algorithm that uses the auxiliary coupling, do not present a better convergence than the synchronous one.

Figure 17 shows that although the global efforts oscillate, throughout the calculation these oscillations occur around the reference value. The response of the structure, Fig. 17b, on the other hand, presents a very good agreement.

2D cylinder flap: FSI3

This test case is particularly interesting since the three axes of the stabilization strategy for the asynchronous approach, 4.3.1, 4.3.3 and 4.3.2, can take place.

Regarding the correction of the dynamics, it was observed that if the added mass operator is not calculated with the adequate structural dynamics, large oscillations appear that lead to the divergence, even in the synchronous case. With respect to the correction used to overcome the additional delay, it was observed that the configurations that showed the best performance for the FSI2 case, when using the parallel synchronous algorithm, do not obtain good results for a strongly coupled application, such as the FSI3 case. Similar to what was done for the FSI2 case, the coefficient a and number of modes were studied to see their influence on the convergence for this case.

While an equivalent influence to that found for the FSI2 case was observed for the number of modes and the under-relaxation of forces, for the amount of correction used, a much more important impact on the response was observed for this case. When setting the parameter a = 1, which means fully applying the correction, without under-relaxation, the solution always led to divergence. In order to improve the convergence for the asynchronous algorithm the test are carried using a relaxation factor ω = 0.53 . Both the displacement correction and the auxiliary coupling are restricted to use only 1 mode.

Fy[N]

Non-linear iterations it is observed that the larger the amount of correction a, the greater the oscillations in the local convergence, both in the local convergence and in the global fluid efforts, figures 19 and 18a. Second, since the amount of correction has an impact on the stability, for a given number of iterations a better accuracy of the global results is observed. According to Fig. 18 the larger the amount of correction, the larger the displacement of the tip.

Still, a positive observation can be made. According to Fig. 20, the fluid efforts obtained with the asynchronous algorithm, and a prediction coefficient a = 0.5, with an average of 13 FSI iterations and an FSI frequency of 7 iterations, present smaller oscillations than those obtained with the synchronous parallel algorithm, with 30 FSI iterations.

DISCUSSION

The present work contributed to the general understanding of the impact of an asynchronous approach in the solution of coupled problems, and in particular for the solution of FSI applications. Being an issue common to the solution of any coupled problem, this work highlighted the fact of a double lag in boundary conditions when using a parallel approach. In fact, knowing that the stronger the influence between the domains the greater the impact of this delay, the need to predict the response in each of the domains can be generalized.

Likewise, particular problems link to the current FSI approach and the models used within this work are pointed out. This is, for example, about the correction in the formulation operator given by the Jacobians at the interface and the impact of the relaxation factors used in the convergence of the non-linearities of the fluid. In all cases, a solution to overcome the issue was proposed, and positive effects were found out of these proposals.

It was sought to make a generalized and qualitative formulation of the problems found so that they could be addressed using the proposed numerical methods or could be improved by different numerical approaches.

With the objective of achieving greater accuracy, greater stability in the solution, and therefore a reduction in the number of non-linear iterations, the challenge of a formulation and a model to calculate the second Jacobian of the interface is raised.

CONCLUSIONS AND FUTURE WORK

It is confirmed that, to stabilize the coupling algorithm by means of the Jacobian of the fluid at the interface, it is necessary to calculate this operator using the structural dynamics that generated the fluid efforts that are being corrected.

A positive impact is observed in the use of an additional predictor/corrector model in order to account for the additional delay presented when using a parallel algorithm. In this particular case, a model within the fluid solver predicts the structural response in advance in order to correct the fluid efforts.

If no modification of the current coupling is carried out, intermediate states of local fluid convergence do not contribute to global convergence. However, using a structural model as an auxiliary coupling model improves the global convergence of the asynchronous algorithm for an FSI problem.

The implemented asynchronous approach was shown to achieve results equivalent to those of the parallel algorithm while reducing the number of FSI iterations required for convergence. Using a different amount of correction, this was observed in particular for the weakly coupled and the strongly coupled test cases.

The findings in this work represent the possibility of hosting more complex structural models within the FSI coupling. They also represent an advance towards the use of an algorithm that is naturally found between the implicit external and the implicit internal algorithm in the FSI domain. However, there are still several subjects that need to be deepened and others that need to be addressed.

To cope with strong coupling applications, the need for a more precise formulation of the operator second Jacobian of the interface, or Jacobian of the structure, is evident. One way to improve the precision of both proposed stabilization methods, and offer greater robustness, could be to transfer a complete or simplified expression of the Jacobian of the structure, from the fluid solver to the structural solver. The use of this operator can be either externally, as is the case in this work, or internally, if the system of equations of the fluid domain is modified. In either case, a criterion must be established to establish the frequency of exchanges of this operator, so that a trade-off is found between precision and computational cost.

The evaluation of the asynchronous methodology and the stabilization methods proposed in this work was carried out within a fictitious context in which computational times were artificially manipulated. To complete the understanding of the proposed methodology, an evaluation process in light of FSI applications for which the inclusion of complex structural models is required should be carried out.

Figure 1 .Figure 2 .

 12 Figure 1. Implicit external algorithm.

Figure 3 .

 3 Figure 3. Illustration of the asynchronous operation of the parallel algorithm. Faster operation of the fluid solver is considered.

Figure 4 .

 4 Figure 4. Convergence criteria for the asynchronous approach.

  . Serial algorithm.

Figure 5 .

 5 Figure 5. Delay of information introduced by the partitioned approach.

Figure 6 .

 6 Figure 6.Effect of the pressure-velocity underrelaxation factor on the convergence of the FSI coupling of the test FSI2. 100 (S-SYN-PF), 210 (P-ASY-PF-PS).

Figure 7 .

 7 Figure 7. Algorithm predicting the deformation of the elastic body through a modal based correction. Resolution internal to the fluid code.

Figure 8 .

 8 Figure 8. Impact of the auxiliary coupling within the local convergence of the fluid. In this figure m stands for the number of modes used for the modal based correction. 211 (P-ASY-PF-NPS, 2itSeq).

Figure 9 .

 9 Figure 9. Comparison of the tip displacement for the three different FSI studied cases.

FyFigure 10 .

 10 Figure 10. FSI2x100 test case convergence. The circles indicate some cases in which local convergence take place. 100 (S-SYN-PF), 210 (P-ASY-PF-PS).

  . Detailed view.

Figure 11 .

 11 Figure 11. Comparison of the global effort obtained for the FSI2x100 test case while using the asynchronous algorithm. 100 (S-SYN-PF), 210 (P-ASY-PF-PS), 211 (P-ASY-PF-NPS, 2itSeq).

  100 with m=NA, r=0.5, it fsi=30(f=1), or=4 211 with m=NA, r=0.5, it fsi=46(f=1.1), or=4 211 with m=NA, r=0.5, it fsi=20(f=2.4), or=4 (a). FSI frequencies of 1.1 and 2m=NA, r=0.5, it fsi=30(f=1), or=4 211 with m=NA, r=0.5, it fsi=15(f=3.3), or=4 211 with m=NA, r=0.5, it fsi=6.5(f=7.7), or=4(b). FSI frequencies of 3.3 and 7

Figure 12 .

 12 Figure 12. Influence of the local convergence of the fluid domain on the global convergence within a non linear loop, FSI2 test case. 100 (S-SYN-PF), 211 (P-ASY-PF-NPS, 2itSeq).

Figure 13 .

 13 Figure 13. Influence of the a factor on the FSI2 test case. 100 (S-SYN-PF), 200 (P-SYN-PF-PS).

FyFyFigure 14 .

 14 Figure 14. Influence of the number of modes in the modal correction. Only the sequential algorithm uses under-relaxation. Results of the three algorithms are fitted to begin at the same time step. The number of indicated nonlinear iterations is used only as a reference of number of nonlinear iterations performed by each algorithm. 100 (S-SYN-PF), 200 (P-SYN-PF-PS), 201 (P-SYN-PF-NPS, 2itSeq).

FyFigure 15 .

 15 Figure 15. Convergence for the 201 parallel algorithm. Relaxation combined with the modal correction. 201 (P-SYN-PF-NPS, 2itSeq).

  m=NA, r=0.5, it fsi=20(f=2.4), or=4 211 with m=1, a=1.0, r=0.5, it fsi=28(f=1.7), or=4 211 with m=1, a=1.0, r=0.5, it fsi=20(f=2.5), or=4 211 with m=1, a=1.0, r=0.5, it fsi=18(f=2.7), or=4 211 with m=1, a=1.0, r=0.5, it fsi=19(f=2.7), or=4 . Asynchronous algorithm perform in all cases better when using the auxiliary coupling.

  100 with m=NA, r=0.5, it fsi=30(f=1), or=4 201 with m=NA, r=0.5, it fsi=30(f=1), or=4 201 with m=1, a=1.0, r=0.5, it fsi=30(f=1), or=4 211 with m=1, a=1.0, r=0.5, it fsi=28(f=1.7), or=4 211 with m=1, a=1.0, r=0.5, it fsi=20(f=2.5), or=4 (b). Synchronous and asynchronous algorithms, with and without modal correction.

Figure 16 .

 16 Convergence comparison for the FSI2 test case. In this figure m stands for the number of modes used for the modal based correction. 100 (S-SYN-PF), 201 (P-SYN-PF-NPS, 2itSeq), 211 (P-ASY-PF-NPS, 2itSeq).

  m=NA, r=0.5, it fsi=30(f=1), or=4 201 with m=NA, r=0.5, it fsi=30(f=1), or=4 211 with m=1, a=1.0, r=0.5, it fsi=28(f=1.7), or=4 211 with m=1, a=1.0, r=0.5, it fsi=20(f=2.5), or=4 211 with m=1, a=1.0, r=0.5, it fsi=18(f=2.7m=NA, r=0.5, it fsi=30(f=1), or=4 201 with m=NA, r=0.5, it fsi=30(f=1), or=4 211 with m=1, a=1.0, r=0.5, it fsi=28(f=1.7), or=4 211 with m=1, a=1.0, r=0.5, it fsi=20(f=2.5), or=4 211 with m=1, a=1.0, r=0.5, it fsi=18(f=2.7), or=4 (b). Tip displacement Y Figure 17. Global results for the FSI2 case when using the proposed modal based auxiliary coupling. In this figure m stands for the number of modes used for the modal based correction. 100 (S-SYN-PF), 201 (P-SYN-PF-NPS, 2itSeq), 211 (P-ASY-PF-NPS, 2itSeq).

Figure 18 .

 18 Figure 18. Global results for the FSI3 case when using the proposed modal based auxiliary coupling. 201 (P-SYN-PF-NPS, 2itSeq), 211 (P-ASY-PF-NPS, 2itSeq).

Figure 19 .

 19 Figure 19. Influence of the amount of correction a in the convergence of fluid efforts. In this figure m stands for the number of modes, r for the relaxation coefficient, and a for the amount of correction used. 201 (P-SYN-PF-NPS, 2itSeq), 211 (P-ASY-PF-NPS, 2itSeq).

  m=NA, r=0.5, it fsi=30(f=1), or=4 211 with m=1, a=0.5, r=0.5, it fsi=13(f=7.5), or=4

Figure 20 .

 20 Figure 20. Global fluid force Fy. 201 (P-SYN-PF-NPS, 2it-Seq), 211 (P-ASY-PF-NPS, 2itSeq).

Table 1 .

 1 Internal implicit algorithms implemented within K-FSI.

Table 2 .

 2 Fluid and structural parameters

	Test case	ρ S [ 10 3 kg m 3 ]	ν S	µ S [ 10 6 kg m s 2 ] [10 3 kg ρ F m 3 ] [10 -3 m 2 ν F s ] [ m U s ]	m F /m S	Stabilization *
	FSI2x100	1000	0.4	0.5	1	1	1	0.001	no
	FSI2	10	0.4	0.5	1	1	1	0.1	ω = 0.5
	FSI3	1	0.4	2.0	1	1	2	1	M pos add
	* Stabilization method used to solve the case using the implicit internal sequential algorithm

Table 3 .

 3 Frequency and number of FSI iterations per time step (FSI2x100).

	Algorithm	100			210				211
	Matrix size	-	0	2000 3000 4000	0	2000 3000 4000 6000
	FSI frequency	1	1	1.2	2	4	1	1.2	1.9	2.7	8
	It FSI	46	45	40	14	12	46	39	26	18	6

Table 4 .

 4 Frequency and number of FSI iterations per time step (FSI2).

	Algorithm				211
	Matrix size	4000 6000 8000 10000 11000 12000 15000
	FSI frequency	1	1.1	2.4	3.3	7.7	5.9	7.1
	It FSI	50	46	20	15	6.5	8.5	7

Table 5 .

 5 Frequency and number of FSI iterations per time step (FSI2). Modal correction using the 1st mode.

	Algorithm		211
	Matrix size	20000 30000 40000 80000
	FSI iteration frequency	1.7	2.5	2.7	12.7
	It FSI	28	20	18	19

The modal based correction used in this work consists of a modification of a modal approach that was available within the ISIS-CFD fluid solver. For a deeper insight the reader can refer to the work of Debrabandere et al.[START_REF] Debrabandere | A staggered method using a modal approach for fluid-structure interaction computation[END_REF].

With the under-relaxation being defined as f i+1 = (1 -ω)f i + ω f i+1 , where f refers to the fluid effort computed by the fluid solver and f to the relaxed fluid effort. The previous and current iteration are represented by i and i + 1, respectively