Towards Unsteady Approach For Future Flutter Calculations

: In prior work to define an improved hydrodynamic approach to flutter calculations, Centrale Nantes, Bureau Veritas Marine & Offshore and Farr Yacht design investigated the possibility of defining a linearized unsteady hydrodynamic model using a fluid response database coming from a series of 2D unsteady RANSE computations. The approach is compared to the Theodorsen theory. The linearized hydrodynamic model was used in a strip theory model for frequency domain flutter analysis. In this latest work, the IMOCA 2006 keel which has been used previously in frequential domain flutter calculation is also analyzed using an alternative and more accurate solution, featuring a fully coupled FSI modal approach with CFD. As described in the literature, the results present a large impact of the unsteadiness on the phase and module for both lift and moment with a fairly good match compared to Theodorsen theory. The implementation of the results on a frequency domain flutter analysis tool reduces the critical speed for the studied model. Thus the results are closer to the 3D modal CFD approach which gave an lower critical speed.


NOMENCLATURE
Dimensionless

INTRODUCTION
At the 2015 HPYD conference, Burns Fallow from North Sails noted how quickly we become accustomed to extraordinary. Indeed five years ago, it was extremely rare for boats to actually 'fly'; only Hydroptere and a few foiling dinghies had this potential. After the 2013 America's Cup the number of fully foiling designs exploded. We are now witnessing the first fully foiling ocean going yachts. It is relevant to consider whether flutter will be an issue for larger and larger foils. Are foils not concerned by the hydro-structural instability, or is it simply that the design envelope is still below the dangerous area? And what is the margin between the actual designs and the flutter limit?
As a reminder, an unstable hydro-structural vibration or flutter caused the failure of the keel of the IMOCA 60 SILL in 2002 (Balze and Devaux 2013). Since then several IMOCA 60's were subject to the phenomenon which resulted in expensive and time consuming redesign efforts. Mouton and Finkelstein (2015) presented the promising results of a modal hydro-structural flutter analysis. This approach used a very basic hydrodynamic model based on the heave and pitch motion of 2D hydrofoil strips. For structures such as keels, the model gave fairly accurate results and improved understanding of the phenomenon. The previously developed model can be used successfully to predict flutter in the sailing domain of a canting keel. For example this model can be used to produce optimized designs for yachts such as mini Transat 6.50 and maxis where the class rules do not require one design keels. In contrast other canting keel classes, like the VOR65 or the present IMOCA 60's have avoided flutter by requiring one-design steel keels.
Unfortunately foils and keels are very different structures in terms of their dynamic response. Hydro-structural effect induced by lead carriers (keels) is very different from the response observed on a high modal frequency foil. Hence, the first step was to study the difference between the two appendages and the available models for unsteady flows.
The present paper relates the effort produced by the LHEEA lab and Bureau Veritas Marine & Offshore to enhance the flutter analysis possibilities from a quasi-static approach to an unsteady model. The first part of the work deals with the unsteady strip model of the hydrostructural modal calculation. The second part describes preliminary results using a fully coupled modal approach with CFD.

Flutter
The phenomenon of flutter is well known in the aviation industry. It was discovered at the beginning of the twentieth century on aircraft wings and was shown to be linked to a vibration problem. During flutter, a coupling of torsional and bending vibration modes of the structure with the aerodynamic forces leads to a transfer of energy from the fluid to the vibrating structure.
Below the flutter speed, and at low vibration frequencies, the hydrodynamic lift induced by a heave motion (bending) is in opposition of phase with the motions speed. It damps the vibration by transferring energy from the vibrating structure to the fluid. The lift induced by the torsional motion is in phase with the motion. It transfers energy from the structure to the fluid and reversely at each period. It is close to be energy conservative. In such condition, the damping is positive and any vibration of the structure reduces.
Approaching flutter speed, with low vibration frequencies, the coupling of the bending and torsional modes allows a mode to be combination of non-phased torsional and bending motions. During a vibration the lift of the torsional component may have a phase different from the bending motion. Thus it can bring energy to the bending component. The damping then becomes null or even negative. This transfer of energy, even for a few seconds, can be enough to induce vibration amplitudes that can lead to complete failure of the structure. The phenomenon is said to be explosive.

Lead bulbs and foils, quasi-static and unsteadiness
The main objective of a modern keel is to carry a heavy lead bulb that will produce a large part of the overturning moment necessary to keep the boat upright under the massive force of the sails. If we consider a typical 2006 IMOCA keel design that experienced the onset of flutter, the bulb inertia lowers the natural frequencies of the keel to: -1.0 Hz for the bending mode, -2.15 Hz for the torsion mode.
The appendage modal vibration creates a harmonic oscillation of the fin's section. A vibration of the torsional mode will see a full period of fin section oscillation in 0.46 second. During this time, a particle of water travelling at 30 knots will travel a distance of 7.2 m. This means that after one oscillation, the particle has travelled 18 chord lengths downstream. In other words, the vortex wake induced by an oscillation of the profile will be very far aft of the profile at the end of the initiating oscillation. As the influence of the vortex is inversely proportional to the distance, it will have a small impact on the profile. The hydrodynamic response can be considered quasi-static. This is characterised by a low reduced frequency equal to 0.3 in this case.
On the other hand, a light and stiff foil will have natural frequencies classically around 20 Hz. For the same chord and speed, a particle of water will travel only 1.1 chord lengths downstream during an oscillation. The vortex wake induced by the oscillation will exist directly downstream of the section of the foil at the end of the period. This proximity creates a potential influence of the vortex wake on the hydrodynamic forces of the next oscillation. The flow is said to be unsteady. The reduced frequency of this example is 5.3. Figure   The number of chord lengths travelled by a particle of water (or a vortex wake) during a period of the oscillation is 2 . This explains why the quasi-static hydrodynamic model can give good results in the flutter calculation of a ballasted keel. However, to enable the faithful modelling of hydro-structural effect on a foil, it is necessary to take into account unsteady hydrodynamic effects.
Although turbulence models mainly designed from steady cases still need validation for such FSI configurations, the unsteady effects are intrinsically considered in a URANSE code. However, the use of complex CFD tools for flutter evaluation is cumbersome as the creation of such a model is complex and time consuming. Therefore, the frequency domain model remains interesting but requires enhancement of its hydrodynamic models to consider high reduced frequency problems.
In order to better take into consideration the unsteady effects, two options are considered and compared. The first is the analytical Theodorsen model. This model considers a 2D wing section harmonically oscillating in a flow. As a function of the reduced frequency, it evaluates the lift amplitude and phase lag. Momentum can also be expressed, but not drag.
The second option is to build a 2D unsteady model representing a foil profile pitching or heaving in a uniform flow produced by CFD. By simulating the fluid forces produced by different sets of harmonic motion in a database, the latter can be used as input to the frequency domain tool in order to incorporate unsteady effects into the fluid model.
The methodology is then based on two tools: the first one, developed by Bureau Veritas takes its background from the aeronautical industry. It is based on the frequency domain analysis of a hydro-structure model, the results of the 2D unsteady calculations are integrated via a strip theory model as hydrodynamic models of the code (Mouton & Finkelstein 2015). The second one, ISIS-CFD, is a Navier-Stokes solver developed at the LHEEA Lab., which is used here both to build a database as input to the previous tool and also to carry out a fully coupled time-marching resolution of this fluid-structure interaction problem.

Modal hydro-structural flutter analysis
The structural model is based on a truncated series of the n most energetic modes, without limitation on the number of modes. Accounting for two modes in a keel analysis is sufficient, but would be insufficient to analyze properly light structures such as a hydrofoil. The inputs to the flutter analysis software is a chosen number of n dry modes natural frequencies and shapes (normalized using the diagonalized mass matrix) provided as an output from Finite Element Analysis (FEA) of the structure.
Laws of dynamics give the following equation for the dry system with or without heel: As a reminder, and as detailed in Mouton and Finkelstein (2015), the strip model consists of several slices of the fin profile along the span. During any vibration of the keel or foil, the motion of each slice can be described with 6 degrees of freedom (see figure (2)) among which two are identified as having a sensible impact on the hydrodynamic flow.
In the local axis of a fin section (span wise, chord wise and perpendicular to the mean plane) the motion of the section profile can be named heave and pitch as described in figure (3). The hydrodynamic loads in the respective direction are the drag, lift and pitching moment. For the sake of simplicity, only the lift will be shown here, but it should be noted that what follows can be reproduced for the other hydrodynamic loads.

2D flat foil section harmonically oscillating in pitch and heave
Let's consider a 2D profile section oscillating of small variations ̃ around an equilibrium position 0 : If we consider harmonic motions around the equilibrium position, the motion can be written as equation (3).
Assuming the uncoupled hydrodynamic response we can write the unsteady lift per unit span through equation (4).
In the hypothesis of small motions around the equilibrium position, we make the assumption that the hydrodynamic response is harmonic and its spectral content is purely in the first harmonic of the initial motion, leading to equation (5): Where 0 pitch ; 0 heave ; ; ℎ are real numbers.
Finally, by introducing α ℎ = − ℎ∞ and, 0 = − ℎ 0 ∞ dimensional analysis leads to the following equation (6): where, Note that for very slow dynamic ( ≪ 1) the result should be in accordance with the potential theory. The phase of the pitch motion lift shall tend to 0 for slow motions (lift in phase with pitch). Then we should have for the derivative of the lift coefficient: It should be noted that the heave motion produces a lift which will tend to be in phase with the heave induce angle ℎ , therefore ℎ   0 k 0. This leads to − 2 ⁄ of phase lag with the heave motion itself.
The unsteady effect will act on the amplitude and the phase lag of the lift response. It is important to note that any steady and unsteady effect (wake, added mass) can be modelled through this approach. It is therefore a practical way to compute a very complete hydrodynamic model for harmonic motions. Following the same approach we get for the pitching moment:

2D CFD Calculations
The objective of the 2D CFD analysis is to compute the amplitude 0 pitch ( ,Re) and the phase ( ,Re) (or the equivalent for the heaving motion) of the hydrodynamic response.
Each motion (pitch or heave) is imposed to a profile in a constant infinite flow defining a couple ( ,Re). The unsteady calculation is run and the lift computed for each time step. Thus, the Fast Fourier Transformation of the resulting signal is achieved. Only the first harmonic is kept here assuming that it contains the major part of the response. This hypothesis is validated in the next section.
To build the database, the profile shape from the 2006 Farr Yacht Design IMOCA 60 fin is used (see figure (4)). Computations were carried out using the ISIS-CFD solver described in section Navier-Stokes solver. The 2D fluid domain of 51610 cells is 12 chord lengths and 10 chord heights. Far field conditions are used for all external boundaries, except the outlet where a constant pressure is imposed. A wall-function approach associated with the − SST Menter turbulence model is used (Menter 1994). Grid and time step independency were checked before running the complete database.  The domain of calculation was 9.6e4 to 1.0e7 for the Reynolds number, and from 0.1 to 6 for the reduced frequency.

Theodorsen solution
The model of Theodorsen (Theodorsen 1935;Hodge and Pierce 2002) is derived from a theory of unsteady aerodynamics for a thin 2D airfoil undergoing small harmonic oscillations in incompressible and inviscid flow.
It contains added-mass and quasi-steady contributions as well as the effect of the wake, included in the Theodorsen's function ( ), for both lift force and momentum. No information concerning the drag can be deduced from this model.
Since the publication of Theodorsen (1935) Hunsaker, Phillips (2015) also used it to analyze the efficiency of flapping foils and compare the results with CFD.
With the notations used in section Modal hydro-structural flutter analysis and in figure (3), the expressions of lift and momentum (expressed at the center of rotation) are given by equation (9).
( ) being the n th Hankel function of second species.

Brief description of the ISIS-CFD code
The ISIS-CFD code, developed by the METHRIC team of the LHEEA Lab. of Ecole Centrale Nantes, UMR-CNRS 6598, solves the incompressible Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations in a strongly conservative way. It is based on a fully unstructured, cell-centered finite-volume method to build the spatial discretization of the conservation equations. Pressure-velocity coupling is obtained through a Rhie & Chow SIMPLE-type method: in each time step, the velocity update comes from the momentum equations and the pressure is given by the mass conservation, transformed into a pressure equation. The temporal discretization scheme is the Backward Difference Formula of order 2 (BDF2) when dealing with unsteady configurations. For each time step, an inner loop (denoted by non-linear loop) associated to a Picard linearization is used to solve the non-linearities of the system and converge all the sequential coupled equations.
In the case of turbulent flows, additional transport equations for modeled variables are solved in a form similar to the momentum equations and they can be discretized and solved using the same principles. An Arbitrary Lagrangian Eulerian (ALE) formulation is used to take into account modification of the fluid spatial domain due to body motion and deformation (Leroyer et al. 2008). Only monofluid configurations are investigated here. However, free surface flowcan be addressed with an interface capturing method, by solving a conservation equation for the volume fraction of water, with specific compressive discretization schemes (Queutey and Visonneau 2007). The code is fully parallel using the MPI (Message Passing Interface) protocol.

Fluid-Structure Interaction coupling with modal approach
The modal approach in the ISIS-CFD code uses the n first natural modes of vibration of the dry structure as input. Since the eigen modal matrix makes the mass matrix and stiffness matrix orthogonal, the temporal resolution of the structure is then reduced to a set of uncoupled degree of freedom (DOF) governed by the n modal equations given by equation (10).
where is the i th modal vector normalized by the mass, the amplitude of the mode, = 2 the pulsation, the frequency, and is a possible damping coefficient assuming the Rayleigh damping hypothesis (damping matrix proportional to the mass and stiffness matrix). ( ) refers to the fluid force acting on the structure. The total shape deformation is then given by equation (11).
These equations are fully coupled with the resolution of the fluid flow through an internal implicit coupling. At each non-linear iteration of the fluid solver (within the same time step), all the DOF of the structure are solved and updated. A Radial Basis Function (RBF) approach is used here to compute the source term of equation (10) and as a mesh deformation technique to recover a body-fitted mesh after each resolution of the structure.
This approach was previously validated by extracting some configurations from Hollowel (1981) and Verma et al. (2015), based on aerodynamic flutter.

Verification of the harmonic content
Using the FFT, it was verified that all disregarded harmonics are less than 5% of amplitude of the fundamental mode.

Verification of the linearity
The first step was to validate the assumption of linearity. Additive and multiplicative linearities were checked. This was performed for various reduced frequencies. Results are very good as shown in figure (6). In this figure, we observe good agreement between the hydrodynamic force generated by a combined motion, and the reconstructed force 2L1+2L2 generated with separated pitch L1 and heave L2 motions at medium to high dynamic motion (k=2.51 et Re=3.96e5).
The slight difference can be attributed to the amplitude of the motion which becomes larger with the combination of both pitch and heave associated with the factor 2 for the amplitude.

Comparison 2D CFD database versus Theodorsen
The results of 2D CFD are compared to the Theodorsen theory for lift and pitching moment in terms of module and phase. The figures presented here show the results as a function of the reduced frequency and for all computed Reynolds numbers. It seems that the Reynolds number has a low influence on the results as it is not possible to differentiate the URANSE results at different Reynolds Number except for . However these results need to be consolidated as the database was calculated on the path (in k and Re) of the modes of the present study. Therefore it should be completed to fully describe the unsteadiness at various reduced frequencies and Reynolds. The computed points are presented in figure (7). Two main reduced frequencies have been assessed for each Reynolds number.

Figure 7. Reduced frequency versus Reynolds points computed in 2D unsteady CFD
Comparisons for lift are given in figures (8,9). Even though ℎ ; are identical for the quasi-static regime, large differences appear for high dynamic regime. The unsteady effect has a large impact on the lift curve coefficient as the minimum is more than 30% below quasi-static value, and for large value of it increases continuously, (figure (8)). Regarding the Theodorsen model, the body motion must produce small variations of velocity compared to the inflow velocity to fulfill the hypothesis of small disturbances. As a consequence, this approach becomes less and less valid when the reduced frequency increases. This explains why we can observe larger differences when the reduced frequency increases.
We observe that both approaches follow very similar trends. For = 0 both curves tend to the potential lift curve slope coefficient 2 of a thin plate assumed by the Theodorsen theory. It should be noted that the CFD calculation presents the actual foil profile and therefore tends to the actual value which is close but not equal to 2 . In figure (9), the phases present a good agreement between the two methods. For low values of phases are in line with the quasi-static, , ℎ tend to 0, which traduce no phase delay with ̃ the pitch motion, and α ℎ the heave induced angle. Concerning the pitching moment, a good agreement can also be noticed between CFD and Theodorsen theory, figures (10, 11). Again the quasi-static value is matched as the pitching moments tend to zero on this profile where the rotation center is defined at one quarter of the chord from the leading edge. This is verified for both pitch and heave. Similar to the lift, the quasi-static pitching moment is not a proper model to define a configuration involving high dynamic phenomena. Regarding the phase, figure (11), for large values of reduced frequency, the models agree fairly well. However, there are large differences at low frequencies. Even though the impact is probably limited as it happens as the amplitude tends to zero, the reason for these discrepancies is not known at this point. It may be related to the viscous friction that becomes more important than the pressure loads. The thickness of the profile may play a role since the Theodorsen model only considers a zero-thickness profile. Also, for high values of the presented results can be influenced by the lower Reynolds number used compared to the ones used for low values of , see figure (7).
It can be noticed the present CFD results which uses a quite thick profile (see Fig.4) match well with the trends obtained by Motta (2015) which study the influence of the thickness of the profile. Even if the study was restricted to a range of reduced frequency from 0 to 0.75, the results of Motta (2015) for quite thick NACA profile (NACA 0018 and NACA 0024) are very similar to the one presented here concerning the position of the curves relative to the Theodorsen ones for the amplitudes and for the phase. It can be noticed that the influence of the thickness of the profile is not negligible.
It seems that the hydrodynamic focal point at a quarter of the chord could be reconsidered at high reduced frequencies as the moment increases largely. However, this results should be further translated in terms of position of the focal point.
2D unsteady hydrodynamic calculations, by nature induce 2D unsteady wakes that tend to produce larger influence of the unsteady effects. In real 3D structure, the wake and vortices will fold in 3D and their influence tend to be reduced. Thus the 2D hydrodynamic model integrated along a fin to represent a 3D flow represents rather overestimated effect of the unsteadiness. The 2006 Farr Yacht Design Imoca 60 keel was analyzed and compared to the previous version of the model which considered quasi-static hydrodynamic with pitch and heave induced angles.
In figure (12) the evolution of the mode is given as a function of boat speed. On the left hand figure the evolution of the frequency is shown (top). The bottom left graph presents the evolution of the reduced frequency of mode 1 and mode 2 and the right hand side graph presents the evolution of the damping rates. The critical speed corresponds the abscissa where the mode 1 curve cut the 0 axis.
Steady and unsteady results are very close for this keel which presents low reduced frequencies as shown at the bottom of the left figure (12) . As shown in the right figure (12), the slight differences change the flutter speed from 20.9 down to 18.9 knots.
As expected no tremendous differences are found on this keel. This is due to the low natural frequencies that induce an almost quasi-static hydrodynamic regime. Here the quasi-static model is sufficient to model correctly the physics of the flow. However, even for such a low reduced frequency the flutter speed is reduced by 10% which is not negligible.
It is different for a lighter structure, and the unsteady effects will have a larger influence. The method is now ready for this future step. Next steps are to include pitching moment into the code, and to allow the definition of curved foils. FSI coupled modal approach with CFD For this configuration, simulations have been done with a full hexahedron mesh of 3.8 million cells designed for a wall-function approach of the turbulence (y+=30). Figure (13) show a global and a closer view of the mesh. The time step is set at 0.002 s. The Reynolds stress tensor is computed using the k-w (SST-Menter) model, (Menter 1994). Concerning the structural part, the information of two first dry modes used for the frequential hydro-elastic analysis was simply converted into an input data file formatted for the CFD code. At the beginning of the simulation, the flow is at rest. The body is sped up to reach its nominal forward velocity within a short time (1s). FSI resolution is activated from the start. Perturbations generated during this violent transient state excite the structure modes. Thereafter, we can follow the evolution of their amplitudes in time and then the evolution of the displacement at some locations of the keel.
Results show that the approach succeeds in catching the flutter occurrence for a velocity around 8 m/s (figure (14)).

Figure 14. Evolution the transversal velocity at the rear bulb tip for different velocities
For velocity of 7 m/s, it can be shown that the time derivative of the rear bulb tip displacement tends to slowly decrease, whereas at 9 m/s, the flutter velocity is exceeded, which results in an increase of the velocity in time. At 8 m/s, close to the flutter critical velocity found by this approach, the response tends to a limit periodic cycle. The range of the flutter velocity found here is a slightly lower than the value found by the hydro-elastic tool. However, additional simulations should be run to assess numerical convergence (grid, time step,...).
The fully coupled time marching approach between modal approach and RANSE solver is capable of detecting flutter phenomenon even if it requires large CPU time to accurately determine the range of the flutter velocity.
It is interesting to consider the differences between the methods presented in this work. As the dichotomic approach of the CFD 3D simply furnished upper and lower bound of the flutter critical speed we will consider the average value for this technique and compare to the frequency domain approach. Figure 15 summarizes the results. It can be seen that there is quite a large difference (around 26%) between the quasi-static hydrodynamic frequency approach and the 3D CFD modal calculation. This is reduced to 14% when considering the 2D unsteady linearized hydrodynamic in the frequency tool based on the CFD matrix response of the 2D profile.
The actual flutter speed on this boat was estimated around 18 to 20 knots. This suggests that the frequency domain approach is more precise than the modal CFD. However, in both analysis, the structural damping was not modelled. As explained in (Mouton & Finkelstein 2015) the structural damping pushes the flutter speed higher as the negative damping from the hydrodynamics need to counter the positive structural damping. The impact on the critical speed was not assessed in the present work. Not considering the structural damping is conservative. Thus the quasi-static frequency domain approach is not conservative while CFD is probably the most accurate.

CONCLUSION
Two main approaches of flutter evaluation have been investigated in this article. The first one corresponds to the frequency domain analysis which was developed by Bureau Veritas M&O. It requires a linearized 2D hydrodynamic model. The work presented in this paper uses CFD URANSE computations to set up a linearized unsteady hydrodynamic model. The linearity of the model requires verification, and is only valid for small motions around the equilibrium position. However, it can include unsteady effects, such as phase delay, amplitude modifications, and added mass in a very simple linear approach. The comparison with Theodorsen's theory allows good confidence in the results. The presented database is sparse, and further work is required. However, the method allows the use of the frequency flutter analysis tool for foils with high dynamic regime even if only driven by the Theodorsen model maybe corrected for thickness as it is suggested in Motta (2015). In the present test case the 2D unsteady model presents a lower flutter speed. However it is not sufficient to generalize.
The methodology is efficient to evaluate the onset of flutter which should be avoided in a design approach. In some cases limit cycles stabilize the flutter effect for large vibration amplitude. This constitutes a limit case and many further investigations will have to be performed prior to accept such effect in the sailing domain. On the contrary, a comfortable speed margin should be used.
The modal approach that is fully coupled with CFD, contains fewer assumptions and captures more of the physics, such as 3D effects. This approach requires running several unsteady simulations at different speeds via a dichotomic search to determine the critical flutter velocity range. More accurate estimations of the flutter phenomenon can be expected, at the expense of far larger computational power.
The tools are complementary; on the one hand the frequency analysis is able to investigate the onset of flutter for a large number of designs in a limited amount of time, on the other hand the FSI computation can provide a precise verification including a broader range of physical phenomenon. It can be noted that the use of linearized unsteady 2D hydrodynamic model in the frequency domain provides a critical speed closer to the URANSE calculated value.