Adrien Martin
email: adrien.martin@doctorant.ec-lyon.fr

Fabrice Thouverez

Dynamic Analysis and Reduction of a Cyclic Symmetric System Subjected to Geometric Nonlinearities

Keywords:

The search for ever lighter weight has become a major goal in the aeronautical industry as it has a direct impact on fuel consumption. It also implies the design of increasingly thin structures made of sophisticated and flexible materials. This may result in nonlinear behaviours due to large structural displacements. Stator vanes can be affected by such phenomena, and as they are a critical part of turbojets, it is crucial to predict these behaviours during the design process in order to eliminate them. This paper presents a reduced order modelling process suited for the study of geometric nonlinearities. The method is derived from a classical Component Mode Synthesis with fixed interfaces, in which the reduced nonlinear terms are obtained through a STEP procedure using an adapted basis composed of linear modes completed by modal derivatives. The whole system is solved using a harmonic balance procedure and a classic iterative nonlinear solver. The application is implemented on a schematic stator vane model composed of nonlinear Euler-Bernoulli beams under von K àrm àn assumptions.

Introduction

The aim of this paper is to describe the analysis of periodic structures. They have already been studied for a long time [START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF] by focusing on their symmetry. The context of non-conservative nonlinearity, such as friction, is well known with the use of several vibration mitigation devices on the bladed disks. An overview of the methods used for the study of such phenomena has been proposed recently by Krack [START_REF] Krack | Vibration prediction of bladed disks coupled by friction joints[END_REF]. The scope of polynomial nonlinearity was considered in the 1990s by Samaranayake [START_REF] Samaranayake | Amplitude modulated dynamics and bifurcations in the resonant response of a structure with cyclic symmetry[END_REF], Vakakis [START_REF] Vakakis | Non-linear normal modes (NNMs) and their applications in vibration theory: an overview[END_REF], and more recently by Georgiades [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF] and Grolet [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF], in the context of nonlinear normal modes (NNMs) [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF], as they exhibited bifurcations in the case of simplified models composed of Duffing oscillators.

It is necessary to reduce such structures due to the size of the finite element (FE) model. Symmetry can be used to considerably reduce size by considering a reference sector and using it to determine the response of the other sectors [START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF]. The study proposed herein deals with a system subjected to nonlinear effects which means that the response of the structure can no longer be expressed with independent diameters, even in the case of a system with natural cyclic symmetry. However, it is possible to use classical Component Mode Synthesis (CMS) such as the Craig & Bampton [START_REF] Craig | Coupling of Substructures for Dynamic Analyses[END_REF] method and consider each sector as a substructure. Paired with NNMs, this method has been used for the study of friction nonlinearity [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF][START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF], and in the case of cubic nonlinearity [START_REF] Apiwattanalunggarn | Component mode synthesis using nonlinear normal modes[END_REF].

This paper proposes an original coupling of several methods used for the study of nonlinear systems. The nonlinear effects of each substructure are reduced by using linear normal modes (LNMs) completed with modal derivatives (MD) [START_REF] Idelsohn | A reduction method for nonlinear structural dynamic analysis[END_REF][START_REF] Slaats | Model reduction tools for nonlinear structural dynamics[END_REF]. The reduced nonlinear stiffnesses associated with the basis are determined using a STiffness Evaluation Procedure (STEP) [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF]. With this basis we are able to derive an original nonlinear CMS which can be used for our simluations. The originality of this work also lies in its application to a periodic structure. Indeed, the innovative approach proposed allows determining the nonlinear response of such specific structures within a reasonable computation time, demonstrating the improvement it provides. The model studied is a periodic stator vane modeled with nonlinear Euler-Bernoulli beams in conformity with the von Kàrmàn simplification of the Green-Lagrange tensor.

Reduction Procedure

The study is carried out on a n dimensional system of nonlinear equations

M ẍ(t) + C ẋ(t) + K x(t) + f nl (x(t)) = f ext (t) (1)
where x(t) is the vector of unknowns, i.e. the degrees of freedom (DOFs); M, C and K refer to the mass, damping and linear stiffness matrices; f nl (x(t)) describes the nonlinear forces due to the geometric effects and thus depends only on the displacements of the structure; and f ext (t) is the vector of external forces applied on the structure.

Fixed-interface Component Mode Synthesis: Craig-Bampton reduction

Let us consider the underlying linear system of Eqn. [START_REF] Apiwattanalunggarn | Component mode synthesis using nonlinear normal modes[END_REF], written as

M ẍ(t) + C ẋ(t) + K x(t) = f ext (t) (2)
where the nonlinear forces have been omitted. The reduction method presented herein lies on a classical Craig-Bampton procedure [START_REF] Craig | Coupling of Substructures for Dynamic Analyses[END_REF]. Thus, the matrices are partitioned into internal (x i ) and boundary (x b ) DOFs, and the Eqn. ( 2) can be re-written with respect to the ordered DOFs y

M ii M ib M bi M bb ÿ(t) + C ii C ib C bi C bb ẏ(t) + K ii K ib K bi K bb y(t) = f ext i f ext b ( 3
)
Each substructure is then reduced with the fixed-interface eigenvectors Φ j i associated with the internal DOFs (Eqn. ( 4)) and the static modes Ψ j b , with the coupling between the substructures obtained by successively applying a unitary displacement on one of the boundary DOF while the others are clamped (Eqn. ( 5)).

K j ii -ω 2 M j ii Φ j i = 0 (4)
Ψ j b = -K j ii -1 K j ib ( 5
)
The transformation from the physical coordinates y(t) to the reduced coordinates z(t) is done using matrix T CB

y(t) =      x 1 i (t) . . . x N S i (t) x b (t)      =      Φ 1 i Ψ 1 b . . . . . . Φ N S i Ψ N S b I b           q 1 i (t) . . . q N S i (t) x b (t)      = T CB z(t) (6)
Thus, the projection of Eqn. (3) onto the subspaces spanned by T CB leads to the reduced system

M CB z(t) + C CB ż(t) + K CB z(t) = f ext CB (t) (7)
where A CB = T CB A T CB for the matrices, and a CB = T CB a for the vectors.

Computation of the reduced nonlinear terms 1.2.1 Formulation and hypothesises

Two hypothesises are made to extend the CMS procedure to the nonlinear forces. The first one, already used in the context of polynomial nonlinearity [START_REF] Apiwattanalunggarn | Component mode synthesis using nonlinear normal modes[END_REF][START_REF] Nguyen | Non linear dynamic of coupled mechanical systems: model reduction and identification[END_REF], assumes that the static modes, and thus the DOFs associated with them, remain linear. It has shown its reliability in the framework of Duffing oscillators [START_REF] Apiwattanalunggarn | Component mode synthesis using nonlinear normal modes[END_REF][START_REF] Nguyen | Non linear dynamic of coupled mechanical systems: model reduction and identification[END_REF]. The second hypothesis assumes that the nonlinear terms depend only on the reduced coordinates associated with the internal modes. This results in expressing the reduced nonlinear forces as

f nl CB (z(t)) =       Φ 1 i f 1,nl Φ 1 i q 1 i (t) . . . Φ N S i f N S ,nl Φ N S i q N S i (t) 0       (8)

Completing the internal modes: modal derivatives

The basis used for the projection of the nonlinear terms is critical for the precision of the reduction. Some of them have already been studied and used. LNMs have shown their capacities to take into account the membrane-flexion coupling [START_REF] Lazarus | Finite element reduced order models for nonlinear vibrations of piezoelectric layered beams with applications to nems[END_REF] induced by the nonlinearities, but they require numerous vectors to do so. Regarding the procedure for projecting the nonlinear terms onto the basis presented further on, it must be kept as small as possible. Thus, the choice was made to use modal derivatives [START_REF] Idelsohn | A reduction method for nonlinear structural dynamic analysis[END_REF][START_REF] Slaats | Model reduction tools for nonlinear structural dynamics[END_REF], due to their capacities to take into account coupling effects along with LNMs. Their mathematical origin lies in the fact that NNMs depend on the amplitude of the DOFs. Using a Galerkin procedure, the solution can be developed by applying Taylor's expansion

x(q) = Φ(q) q = x(q = 0) + r ∑ k=1 ∂ x ∂ q k q=0 q k + 1 2 r ∑ k,l=1 ∂ 2 x ∂ q k ∂ q l q=0 q k q l (9)
Identifying each term in Eqn. ( 9) [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF] allows writing the solution naturally as a sum of a set of LNMs and modal derivatives.

∂ x ∂ q k q=0 = Φ k + r ∑ p=1 ∂ Φ p ∂ q k q=0 q p = Φ k ∂ 2 x ∂ q k ∂ q l q=0 = ∂ Φ k ∂ q l + ∂ Φ l ∂ q k + r ∑ p=1 ∂ 2 Φ p ∂ q k ∂ q l q=0 q p = ∂ Φ k ∂ q l + ∂ Φ l ∂ q k = ∂ Φ kl
x(q) = Φ ∂ Φ q p ( 11
)
where p refers to the unified variables q k q l associated with the MDs.

The computation of the MDs can be performed in several ways [START_REF] Slaats | Model reduction tools for nonlinear structural dynamics[END_REF]. Analytic formulas are proposed which take the mass into account or not and they are well-suited for simple problems where the jacobian of the nonlinear terms is known. A numerical solution is also proposed that relies on a finite difference method between two eigenvalue problems (Eqn. ( 12)): the LNMs of the structure around an equilibrium position Φ, x = 0 in our case, which is not necessarily a stable equilibrium position; and the LNMs of the structure deformed along one of the previous LNMs Λ(Φ k ) and which take into account the nonlinear effects, i.e. computed using the tangent stiffness.

∂ Φ k ∂ q l = Φ k -Λ k (δ h l Φ l ) δ h l ( 12
)
This method is particularly well-adapted to the computation of MDs using a commercial FE solver.

Projection of the nonlinear forces onto the basis : STEP procedure

The projection of the nonlinear forces can be time-consuming as it is necessary to re-build the solution vector, construct the nonlinear terms and then project them onto the reduction basis. One solution is to estimate the reduced stiffnesses associated with the projection basis a priori, and compute the nonlinear forces directly as proposed by the STiffness Evaluation Procedure (STEP) [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF].

Integrating the Green-Lagrange tensor in the FE method [START_REF] Zienkiewicz | 5 -geometrically non-linear problems -finite deformation[END_REF] leads to formulating the nonlinear internal forces as a third degree polynomial

f tot p (x) = n ∑ k=1 K pk x k + n ∑ k=1 n ∑ l=k K q pkl x k x l + n ∑ k=1 n ∑ l=k n ∑ m=l K c pklm x k x l x m (13)
where K q and K c indicate the quadratic and cubic stiffness matrices. Projecting Eqn. ( 13) onto the basis Φ leads to the reduced internal forces

f tot red,p (x) = Φ p f tot (Φ q) = r ∑ k=1 β pk q k + r ∑ k=1 r ∑ l=k β q pkl q k q l + r ∑ k=1 r ∑ l=k r ∑ m=l β c pklm q k q l q m (14)
where β , β q and β c are the reduced linear, quadratic and cubic stiffness matrices. The method proposed by Muravyov [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF] allows determining the reduced stiffnesses through several static loads along linear combinations of vectors of the basis. For example, to compute the terms of the form β pi , β q pii and β c piii , at least three loadings are required along the vector Φ i at different amplitudes, providing the necessary equations (Eqn. [START_REF] Samaranayake | Amplitude modulated dynamics and bifurcations in the resonant response of a structure with cyclic symmetry[END_REF]). Note that β can be obtained simply by projecting the linear stiffness K onto the basis Φ. Similarly, β q pij , β c pijj and β c pii j are obtained through loadings along h i Φ i + h j Φ j with differents h i and h j (Eqn. ( 16)), and

β c pi jk with h i Φ i + h j Φ j + h k Φ k (Eqn. ( 17
)
). In all the cases i j k.

f tot red,p (h i Φ i ) = Φ p f tot (h i Φ i ) = β pi h i + β q pii h 2 i + β c piii h 3 i (15) Φ p f tot (h i Φ i + h j Φ j ) = β pi h i + β q pii h 2 i + β c piii h 3 i + β p j h j + β q p j j h 2 j + β c p j j j h 3 j + β q pi j h i h j + β c pii j h 2 i h j + β c pi j j h i h 2 j ( 16
)
Φ p f tot (h i Φ i + h j Φ j + h k Φ k ) = β pi h i + β q pii h 2 i + β c piii h 3 i + β p j h j + β q p j j h 2 j + β c p j j j h 3 j + β pk h k + β q pkk h 2 k + β c pkkk h 3 k + β q pi j h i h j + β c pii j h 2 i h j + β c pi j j h i h 2 j + β q pik h i h k + β c piik h 2 i h k + β c pikk h i h 2 k + β q p jk h j h k + β c p j jk h 2 j h k + β c p jkk h j h 2 k + β c pi jk h i h j h k (17)
Prediction from previous point

z = z init r = Z CB (ω) z + f nl CB ( z) -f ext CB z r < ε Prediction to next point z(t) = F z f nl CB ( z) = F f nl CB (z(t)) yes no z = z iter+1 Evaluation of f nl CB ( z) : AFT Procedure Nonlinear Solver f nl CB
Eq. 8

Figure 1. HBM procedure of the reduced system

Frequency Analysis: the Harmonic Balance Method (HBM)

The complete and reduced solutions of Eqn. ( 1) are studied in the framework of a frequency analysis. Their limit cycles are searched by assuming a solution of the form

x(t) = a 0 + N h ∑ k=1 (a k cos (kωt) + b k sin (kωt)) (18)
Substituting Eqn. [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF] in the general equation, complete Eqn. (1) or reduced Eqn. [START_REF] Idelsohn | A reduction method for nonlinear structural dynamic analysis[END_REF], and projecting it onto the basis of the harmonic function with respect to the adapted scalar product

f (t), g(t) = ω π 2π ω 0 f (t) g(t) dt (19)
leads to the general equation of the HBM

Z(ω) x + f nl ( x) = f ext ( 20
)
where x, f nl and f ext are the Fourier coefficients of x, f nl and f ext , while Z is the dynamic stiffness expressed as

Z =      2K Z 1 . . . Z N h      , Z k = K -(kω) 2 M kωC -kωC K -(kω) 2 M (21)
The Fourier coefficients of the nonlinear terms are not easily computable analytically, except for very simple systems. Thus, an Alternating Frequency-Time (AFT) method [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF] is used to determine them as it allows computing nonlinear terms regardless of the nonlinearity involved. The method consists in passing forwards and backwards in the frequency and time domains using direct (F) and inverse (F) discrete Fourier transformation matrices. A flow chart of the method is presented in (Fig. 1).

Dynamic Analysis and

Application

Presentation of the model

The structure considered for the application of the method proposed is a schematic stator vane composed of 6 blades (• bl ) and 2 rings (• c ) providing the coupling between the blades (Fig. 2a). A beam model is proposed Fig. 2b, and the repeating sector is shown in Fig. 2c. The dimensions of the beams representing the blades and the rings are given in (Tab. 1) (Length, width, thickness), and the finite element discretization in (Fig. 2d, Tab. 1). The springs can be assimilated with locating pins used to fix the vane on the carter. Thus, the complete model is composed of 726 DOFs.

Nonlinear Euler-Bernoulli theory

The beams are modelled using a classical Euler-Bernoulli theory, with u being the longitudinal displacement, v the transversal displacement and θ the rotation of the sections. The general displacements are thus expressed as

u(x, y) = u(x) -y θ (x) v(x) (22)
The nonlinearities are a consequence of the large strains of the beam, which leads to considering the nonlinear part of the Green tensor.

ε = 1 2 (∇u + ∇ u) + 1 2 ∇u ∇ u (23)
The constraint along the longitudinal displacement of the model allows applying the von Kàrmàn hypothesis [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF] and simplifying Eqn. ( 23)

ε = ε xx = u ,x -y v ,xx + 1 2 v 2 ,x (24)
Using Eqn. (24) in the strain energy

E p = 1 2 L 0 ES u 2 ,x + u ,x v 2 ,x + 1 4 v 4 ,x + EI v 2 ,xx dx (25)
where S = w × t and I = w × t 3 12 , and applying Lagrange's equations gives the linear and nonlinear stiffness matrices. The linear mass matrices are determined thanks to the kinetic energy. The material used is a classical steel (E = 210 GPa, ρ = 7800 kg/m 3 ).

Study of the underlying linear system

Complete system

Let us consider the system described by Eqn. (2). In the case of cyclic symmetry, the system will be subject to simple and double eigenvalues, characterized by nodal diameter mode shapes [START_REF] Samaranayake | Amplitude modulated dynamics and bifurcations in the resonant response of a structure with cyclic symmetry[END_REF]. This study focuses on the first flexural mode family of the blades. The mode shapes are represented in Fig. 3, where the amplitude is taken at the middle of each beam transversally. Classical results are obtained, with 0 and 3 diameters (as there are 6 sectors) as simple modes and 1 and 2 diameters as double modes.

The forced response of the system is computed in order to account for all the diameters. Thus, the excitation is applied at the middle of one of the beams transversally. The excited beam is numbered 1, the other following beams are numbered consecutively. A Rayleigh damping is applied that takes into account the mass matrix and a damping ratio of 0.1% on the 3 diameter modes. The result is plotted in Fig. 4. The graph represents the transversal amplitude at the middle of each beam.

The 0 diameter mode is clearly identified with all the beams having the same amplitude (ω ≈ 6900 rad/s). The other diameter modes, especially the 3 diameters (1 st peak), do not appear so clearly, certainly due to the proximity of those 3 modes.

Linear Craig & Bampton reduction

As the simulations are conducted for the first flexural mode family, the reduction of each substructure is performed with the first flexural mode Φ 1F . In order to process the nonlinearities in the followings parts, a modal derivative is added to the LNM. To simulate a simple system like a beam, the modal derivative ∂ Φ 11 is sufficient, in addition to the LNM, to take the nonlinear effects into account [START_REF] Grolet | Model Reduction for nonlinear vibration analysis of structural systems using modal derivatives and stiffness evaluation procedure[END_REF]. Therefore the basis for the internal modes of the substructures is Φ i = Φ 1F ∂ Φ 11 . The reduction allows reducing the number of unknowns from 726 to 48 (6 × 2 internal modes plus 2 × 3 statics modes for each of the 6 junctions between substructures), equivalent to 93.4 %. The modal derivative ∂ Φ 11 was computed with δ h = 1e -3 (Eqn. ( 12)).

The results of the linear reduction (Fig. 5 and Tab. 3) show good correlation with the complete model with less than 7e -3 % error on the frequencies. Thus, as this basis shows its reliability in the linear context, it will be used again in the nonlinear reduction procedure.

Simulation of the full nonlinear system

Let us assume that only the thin blades are subjected to nonlinearities. As the rings are thick compared to the blades, they do not undergo such effects. The solutions presented in Fig. 6 are computed on the full system with respect to the HBM. The solutions for 3 and 5 harmonics are compared and they are perfectly superposed, the 3 harmonics solution can thus be considered as harmonically converged. In the following, the solutions presented will be computed with 3 harmonics.

In order to present the global results, the following outcomes are plotted in a pulsation-energy diagram [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF] such that

E(x) = a 0 2 + N h ∑ k=0 a k 2 + b k 2 ( 26
)
An example is given in Fig. 7 with responses of the system to several excitations. As the amplitude of excitation increases, the nonlinearity becomes increasingly visible (Fig. 7). From F=0.6N, the system seems to follow a bifurcated branch, like those preesented by Grolet [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF]. The time series in Fig. 8b-8c, representing the transversal amplitude at the middle of each beam as a function of time, clearly highlights that the amplitude is mostly localized on the excited beam, number 1, in the nonlinear case (Fig. 8b), whereas the linear time response is similar to the 1 diameter linear mode (Fig. 8c). In addition, still referring to the time series (Fig. 8b), beam 1 is in phase opposition with all the other beams, which totally breaks the nodal diameter characteristic of the mode (Fig. 3c). In the following, the time series will be taken at the resonance peak of the 1 diameter mode, for both linear and nonlinear signals.

The behaviours in play demonstrate the usefulness of the study and the need to perform it within an acceptable time. Although the analysis is performed along a simplified structure, the phenomena exhibited appear at moderate amplitudes and highlight a strong localisation on the beam 1, clearly showing the limitations of the linear approach and the risks subjected by the structure.

Reduction of the model

The basis presented in the linear study is re-utilized for the reduction of the nonlinear system. The whole nonlinear procedure presented previously is applied, and the nonlinear terms of each substructure are projected using the STEP procedure on the basis and then assembled along the boundary DOFs. The reduced stiffnesses are computed with h i = 1e -3 and h j = 5e -3 (Eqns. [START_REF] Samaranayake | Amplitude modulated dynamics and bifurcations in the resonant response of a structure with cyclic symmetry[END_REF][START_REF] Slaats | Model reduction tools for nonlinear structural dynamics[END_REF]) (2 coefficients are needed as there are 2 vectors in the basis). The results of the simulations can be seen in Fig. 9 and the correlation between the full and the reduced models is quite good, whether for the bifurcated diameter mode or the unbifurcated ones (Figs. 9a). The behaviour of the structure around the bifurcation is also well reproduced (Figs. 9b-9c) compared to the full nonlinear response. The computation times are presented in Tab. ( 4) and demonstrate a significant gain thanks to the reduction, more than 99 %. It is noteworthy that the whole offline computation, including the computation of the linear mode, the modal derivative and the STEP coefficients, lasts less than 0.2 s.

Parametric study: impact of the rings

Thanks to the capacities of the method, a parametric study is performed on the thickness of the rings. The thickness was chosen in the range -40% +40% of the initial model. The results in terms of forced responses and time series of the 1 diameter mode are visible in Fig. 10. The first impact of the rings' thickness on the structure is observable on the linear frequencies of the diameter modes. Indeed, the thinner the rings are, the farther the frequencies are from each other due to the strength of the coupling performed by the rings.

This also has a clear impact on the nonlinear response of the system. For a weak coupling, i.e. t c = 0.6 cm, the system appears to show no bifurcation, as highlighted in Fig. 10a. This is supported by the time series in Fig. 10b where the 1 diameter mode appears clearly. For stronger coupling, t c > 0.6 cm, the 1 diameter mode reaches the bifurcated branch, as shown in Figs. 10c-10e-10g and the phase opposition of the beam 1 with the others (Figs. 10d-10f-10h). For approximately the same energy level, the bifurcation of the 1 diameter mode is thus expected to occur for a coupling thickness between 0.6 cm and 0.8 cm.

The time series Figs. 10b-10d-10f-10h also leads to a configuration where the amplitude localization is particularly marked for t c = 0.8 cm. In this case, the risk of reaching the fatigue limit is more significant than the others, highlighting the impact that such nonlinear localization phenomena can have on the design process.

Conclusion

This work was devoted to the presentation of a reduction procedure suited to the study of structures subjected to geometric nonlinearities. The method is based on a component mode synthesis, completed by the use of modal derivatives and a stiffness evaluation procedure in order to reduce the nonlinear terms. The application was tested on a schematic stator vane composed of nonlinear beams. The results showed that the procedure is efficient in taking into account the nonlinear behaviours of the structure, especially the specific dynamical effects associated with cyclic symmetry such as energy localization. What is more it provides an interesting gain in time.

Besides the reduction procedure, this paper illustrated the extreme behaviours that a cyclic system can be exposed to. The nonlinear effects exhibited strong localizations on the structure at low amplitudes which must be avoided to ensure the integrity of the components. The method proposed here allowed performing this study in accordance with the imperative of reducing design time, for example, with respect to the thickness of the rings, and highlighting the presence or not of localization and its significance.

In further works, it would be interesting to perform such a study with NNMs and compare it with the results of the forced responses presented in this paper. This would allow predicting the bifurcation points more easily in a design procedure. Also, the extension of the method to a more realistic and complex system, such as a 3D bladed disk, is planned to assess its real potential. Lastly, coupling this method with that proposed by Joannin [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF] to incorporate friction non-linearity in the model will also be considered.

Figure 2 .

 2 Figure 2. Stator vane model

Figure 3 .

 3 Figure 3. Mode shapes of the first flexural mode family

Figure 4 . 11 Figure 5 .Figure 6 .

 41156 Figure 4. Linear forced response of the system for F=1N

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Pulsation-Energy diagram of the full nonlinear system for several forcing amplitudes (Linear: --, Nonlinear: 0.1N, 0.3N, 0.6N, 0.8N, 1.1N)

Figure 10 .

 10 Figure 10. Influence of the thickness of the rings (Linear: --, Nonlinear: 0.1N, 0.5N, 0.8N, 1.1N)

Table 1 .

 1 Properties of the reference sector

	L c	w c	t c	N c	L bl	w bl	t bl	N bl	k
	0.1 m 0.02 m 0.01 m	6	0.1 m 0.02 m 0.002 m	30	1e5 N.m -1

Table 2 .

 2 Pulsations of the first flexural mode family

	Nodal diameters	3	2	1	0
	Pulsation [rad/s] 6638 6662.3 6676.7 6907.4

Table 3 .

 3 Error on the linear reduction

	Nodal diameters	3	2	1	0
	Reference puls. [rad/s]	6638	6662.3 6676.7 6907.4
	Reduced puls. [rad/s]	6638.1 6662.3 6676.7 6907.9
	Relative error [%]	5.4e -4 1.4e -4 1.7e -5 6.6e -3

Table 4 .

 4 Gain in time of the reduction

	0.6	11927	89.97	99.246
	0.8	9264	44.4	99.5
	1.1	10243	73.2	99.286

F [N] Full [s] Reduced [s] Gain in time [%]


Guy Archambault

Rock Flamand

Sylvie Gentier

Joëlle Riss

Colette Sirieix

Rock joint shear mechanical behavior with 3D surfaces morphology and degradation during shear displacement

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


est destinée au dépôt

The behavior of induced pore fluid pressure in undrained triaxial shear tests on fractured porous analog rock material specimens.

Mechanics of Jointed and Faulted Rock (MJFR-3) 1998, Rossmanith ed .. A.A. Balkema, pp.555-560.

  The behavior of induced pore fluid pressure in undrained triaxial shear tests on fractured porous analog rock material specimens Guy Archambault, Stephane Poirier, Alain Rouleau, Sylvie Gentier, Joëlle Riss To cite this version: Guy Archambault, Stephane Poirier, Alain Rouleau, Sylvie Gentier, Joëlle Riss. The behavior of induced pore fluid pressure in undrained triaxial shear tests on fractured porous analog rock material specimens. Mechanics of Jointed and Faulted Rock (MJFR-3) 1998, pp.1-6, 1998. hal-03740177


Corentin Lothode
email: corentin.lothode@univ-rouen.fr

Jules Poncin

Didier Lemosse
email: didier.lemosse@insa-rouen.fr

David Gross

Eduardo Souza De Cursi
email: eduardo.souza@insa-rouen.fr

E Souza

Investigation of blade-mast fluid-structure interaction of a tidal turbine

In this article, we investigate the movement and vibrations of a blade due to the presence of the mast. When the blade passes in front of the mast, a sudden pressure spike induces vibrations in the blade.

. In a third section, we present cases of fluid-structure interaction. Several structures are setup by modifying the stiffness of the material. Their steady open-water (without a mast) behaviors are compared. And finally, two dynamic fluid-structure computations are performed to compare the behavior of an elastic blade passing next to a mast. For all the cases, we use K-FSI developed by K-Epsilon to solve the fluid-structure interaction (FSI).

Introduction

Tidal currents as a potential source of renewable energy has been under study for more than 10 years now. Mainly horizontal tidal turbines have been investigated both numerically and experimentally [START_REF] Bahaj | Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank[END_REF][START_REF] Mycek | Experimental study of the turbulence intensity effects on marine current turbines behaviour. part i: One single turbine[END_REF][START_REF] Tedds | Near-wake characteristics of a model horizontal axis tidal stream turbine[END_REF][START_REF] De | The effects of wave-current interaction on the performance of a model horizontal axis tidal turbine[END_REF][START_REF] Fernandez-Rodriguez | Experimental study of extreme thrust on a tidal stream rotor due to turbulent flow and with opposing waves[END_REF]. Vertical tidal turbines have also being studied [START_REF] Stallard | Interactions between tidal turbine wakes: experimental study of a group of three-bladed rotors[END_REF][START_REF] Derakhshan | Experimental and numerical study of a vertical axis tidal turbine performance[END_REF]. Two designs were also compared extensively in the literature, see [START_REF] Khan | Hydrokinetic energy conversion systems and assessment of horizontal and vertical axis turbines for river and tidal applications: A technology status review[END_REF] for example. Current research focuses on optimizing array of tidal turbines and their interaction [START_REF] Funke | Tidal turbine array optimisation using the adjoint approach[END_REF], erosion of the seabed [START_REF] Verbeek | Estimating the stability of a bed protection of a weir-mounted tidal turbine[END_REF], performance and optimization of the shape or structure (Nicholls-Lee, 2011), fatigue and reliability [START_REF] Finnegan | Operational fatigue loading on tidal turbine blades using computational fluid dynamics[END_REF].

To optimize the performance of a tidal turbine, there is a need for accurate and efficient tools to perform simulations. In this article, we focus on the development and evaluation of a new tool to better understand the physics and calculate the loads acting on a single flexible tidal turbine including mast effects. This tool could be used to improve estimations of fatigue for example.

The first simulation framework used to perform computations was the Blade Element Momentum Theory (BEMT), for example in [START_REF] Batten | Experimentally validated numerical method for the hydrodynamic design of horizontal axis tidal turbines[END_REF]). BEMT is a good approach to assess the open water performance of one turbine, but it fails to accurately simulate the performance in the pres-ence of obstacles or other turbines. To avoid this problem, other approaches have been developed such as the Vortex Lattice Method (VLM) in [START_REF] Pinon | Numerical sim-ulation of the wake of marine current turbines with a particle method[END_REF]. Their focus is the wake of the turbine to study the interaction between two or more turbines (for example in [START_REF] Mycek | Numerical and experimental characterisation of interactions between two marine current turbines[END_REF]). Their results are accurate until stall which is a known limitation of the model: their method force the flow to be attached until the trailing edge, hence overestimating hydrodynamical forces.

Attempts to use Computational Fluid Dynamics (CFD) on wind or tidal turbine has been made in the past. To avoid too much computational efforts, many authors modelled the behavior of the turbine instead of solving the complete geometry. For instance, [START_REF] Jimenez | Advances in largeeddy simulation of a wind turbine wake[END_REF] used Large Eddy Simulation (LES) with the turbine modeled by an approximate model of a concentrated drag force to study the wake development. Also using an approximate model for the turbine, [START_REF] Calaf | Large eddy simulation study of fully developed wind-turbine array boundary layers[END_REF] performed a LES computation using an actuator disk.

Fully resolved blade geometry CFD computations are computationally expensive but can yield much more information about the flow behavior and force distribution along the blade. [START_REF] Afgan | Turbulent flow and loading on a tidal stream turbine by LES and RANS[END_REF] modeled the fluid using different turbulence models (𝑘 -𝜔 SST, Launder-Reece-Rodi and LES) on the 20°pitch angle case of [START_REF] Bahaj | Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank[END_REF]. They performed unsteady simulations including the mast and a simplified geometry of the cavitation tunnel. [START_REF] Yan | Free-surface flow modeling and of horizontal-axis tidal-stream turbines[END_REF] was the first work to perform a full-scale free-surface flow simulation of a tidal turbine. All of those methods need a lot of computing power, and are only slightly more accurate than BEMT if the goal is only to obtain information on performance curves. Their main advantage is their ability to add geometry details (mast, cables, proximity to the seabed or duct for example) that would be difficult or even impossible to model with current BEMT models. To account for the fluid when simulating an unsteady simulation, a linear model is often used. This is, for example, what is used by the industry standard software Bladed and the open-source software FAST [START_REF] Jonkman | Nwtc design codes (fast)[END_REF]. It is with the work of [START_REF] Bazilevs | 3d simulation of wind turbine rotors at full scale. part ii: fluid-structure interaction modeling with composite blades[END_REF] and [START_REF] Takizawa | Numerical-performance studies for the stabilized space-time computation of wind-turbine rotor aerodynamics[END_REF] that the first dynamic FSI computations were performed on wind turbines using CFD on fully resolved geometries.

Tidal and wind turbines are different when considering fluid-structure interaction in the sense that water has a density of three order of magnitude higher than air. The fluidstructure coupling is therefore much greater and comparable to what can be observed for hydrofoils [START_REF] Lothodé | Dynamic fluid structure interaction of a foil[END_REF]. [START_REF] Murray | Added-mass effects on a horizontal-axis tidal turbine using fast v8[END_REF] investigated the influence of the added mass on the response behavior of tidal turbines in order to model it in a BEMT code. They showed that the natural frequency is decreased especially when using light blades. Turbulence influx events could not induce vibration because their frequency was higher than the natural frequency of the blade. Nicholls-Lee (2011) coupled a structure solver with a surface panel code in order to obtain a steady-state solution and to optimize a tidal turbine structure. Zilic De [START_REF] De Arcos | Hydrodynamic modelling of flexible tidal turbine blades[END_REF] studied the performance and deflection of a turbine with flexible blades using two-way coupling, but no dynamic effects were investigated. However, no study of the interaction between the various elements of the tidal turbine has been carried out. Therefore, understanding what happens for a blade-mast interaction including the effects due to flexible blades involves the use of strongly coupled fluidstructure interaction solver. For example, these results could be used to understand the dynamics of passively adaptive blades [START_REF] Nicholls-Lee | Application of bend-twist coupled blades for horizontal axis tidal turbines[END_REF].

We first present our numerical schemes. The fluid solver is described in Section 2.1 and the structure solver is detailed in Section 2.2. Our fluid-structure coupling scheme is shown in Section 2.3 and finally mesh morphing is quickly introduced in Section 2.4. Next, we validate our simulation procedure in Section 3.3, first ignoring blade deformation using [START_REF] Bahaj | Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank[END_REF] as a baseline validation for rigid geometry. Later, we introduce the different structures used in the article and reproduce the performance curve taking into account the flexibility of the blades in Section 3.4. Finally, simulations of two different structures focusing on the blade-mast interaction are discussed in Section 4.

Numerical method

In this section, we describe the different solvers used in this article and the coupling scheme that allows to solve the steady and unsteady simulations. In Figure 1, we show a schematic of the domain and we describe the different notations used in this article in Table 1.

Fluid solver

The fluid solver is ISIS-CFD [START_REF] Deng | Ship flow simulations with the isis cfd code[END_REF]. It is developed by the METHRIC team of LHEEA laboratory and commercialized by Cadence. It is included in FINE™/Marine. ISIS-CFD solves the incompressible unsteady Reynolds Averaged Navier-Stokes equations. The code is fully parallel using the MPI (Message Passing Interface) protocol.

It is formulated through a fully unstructured finite-volume method to obtain the spatial discretization of the conserva- tion equations. Arbitrarily shaped polyhedrons are assembled through a list of faces. All unknown variables are cellcentered.

The temporal integration scheme is the implicit, second order Backward Differentiation Formula (BDF2) scheme when dealing with unsteady configurations. For each time step, the residual of the errors is minimized through an inner loop in order to solve the non-linearities of the system. Pressurevelocity coupling is enforced by a SIMPLE like algorithm: at each time step, the velocity updates come from the momentum equation and the pressure is given by the mass conservation law.

Turbulence is accounted using additional transport equations for modeled variables. They are solved in a form similar to the momentum equations and they can be discretized and solved using the same principles [START_REF] Duvigneau | On the role played by turbulence closures in hull shape optimization at model and full scale[END_REF]. In this article, we use the 𝑘 -𝜔 SST model, a two-equations eddy-viscosity model.

An Arbitrary Lagrangian Eulerian (ALE) formulation is used to account for fluid mesh changes due to the body deformation [START_REF] Hughes | Lagrangian-eulerian finite element formulation for incompressible viscous flows[END_REF][START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF]. Basically, the flow velocity 𝒖 fluid is decomposed in two velocities: a velocity of the mesh 𝒖 mesh , and an additional velocity relative to the mesh 𝒖 relative .

𝒖 fluid = 𝒖 mesh + 𝒖 relative (1)
It is possible to have a non-matching sliding interface [START_REF] Queutey | Sliding grids and adaptive grid refinement for RANS simulation of ship-propeller interaction[END_REF]. It allows to have a rotating mesh inside a fixed mesh for example, and to avoid using more costly procedures involved in overlapping meshes. This is very useful for modeling a propeller or a tidal turbine.

A procedure called Rotating Frame Method (RFM) is available. It adds a mesh velocity corresponding to a rotation in a domain whilst the mesh stays fixed. It allows to compute performance curves with larger time steps than typical rotating mesh methods. A comparison of the two methods is made in Section 3.4.

Structure solver

We use the solver K-Struct which is developed by the company K-Epsilon (K-Epsilon, 2020). The solver is based on a non-linear finite element formulation derived through the use of the virtual work principle. The temporal discretization is driven by a Newmark-Bossak scheme [START_REF] Wood | An alpha modification of newmark's method[END_REF] and the resolution is ensured by a Newton method through the computation of the tangent matrix. Non-linear iterations are relaxed with an Aitken acceleration.

The code initially aimed to simulate the dynamic behavior of sailboat rigging: sails, mast and cables [START_REF] Durand | Fsi investigation on stability of downwind sails with an automatic dynamic trimming[END_REF]. A non-linear finite element method with a large displacement formulation is implemented. While numerous element types have been implemented in the structural code, in the present study, only Euler-Bernoulli beam elements are used, as shear deformation was considered negligible. This element accounts for torsion, bending and axial deformations. It offers a good compromise between accuracy and solving speed compared to more complex shell elements. This kind of element is used in the industry standard software FAST [START_REF] Wang | BeamDyn: A high-fidelity wind turbine blade solver in the fast modular framework[END_REF]. The topology and properties of the different structures used in this work will be further detailed in Section 4.

Coupling scheme

We use K-FSI, a FSI software developed by K-Epsilon. K-FSI implements a strong coupling scheme, more details are given in [START_REF] Durand | Interaction fluide-structure souple et legere[END_REF]. As stated previously, the coupling was first developed to tackle sails and rigging issues. Even though air is not a dense material, sails have a large wetted area and are lightweight structures. As shown in [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF], large wetted areas and light structures induce strongly coupled fluid-structure problems. In the case of a tidal turbine, the material used can be dense (typically metals), but the wetted area is significant and the surrounding fluid is denser than air (about 10 3 times more). It follows that the problem is also strongly coupled.

To tackle strongly coupled problems, numerous solutions exist [START_REF] Fernández | Coupling schemes for incompressible fluidstructure interaction: implicit, semi-implicit and explicit[END_REF]. Most of them are partitioned algorithms where the fluid problem and the structure problem are solved separately. The most naive algorithm is often referred as explicit: during one time step, the fluid problem and the structure problem are solved only once, one after each other. It is the fastest and simplest algorithm for weakly coupled problems, but it is unstable in the case of strongly coupled cases. Monolithic schemes, where all the equations are solved together, can be very stable. Since all equations must be developed in a similar framework, and all unknowns are shared between systems, they make it more difficult to use already existing software.

In the case of partitioned schemes used for strongly coupled cases, an implicit scheme with a relaxation (applied to forces given by the fluid or position of the nodes given by the structure) can be stable if the relaxation rate is low enough [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF]. During one time step, the fluid equations and the structure equations are solved one after another until a convergence is reached. It can be very costly if the convergence is not fast enough. The choice of the relaxation rate is a matter of trade-off between greater stability and convergence rate. An Aitken procedure is useful for dynamically adjusting the relaxation rate and accelerating convergence, as in [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF].

Instead of solving a fixed point problem to reduce the residual of errors between the structure and fluid forces, one can use a Newton algorithm including an exact Jacobian operator as in [START_REF] Fernández | A newton method using exact jacobians for solving fluid-structure coupling[END_REF]. It turned out to be very efficient, even though the fluid-structure problem was greatly simplified in order to be able to develop such an operator. K-FSI takes advantage of this idea to accelerate convergence. Such operators have been developed in [START_REF] Lothodé | Modélisation des pales d'éoliennes ou d'hydroliennes en environnement naturel à l'aide d'un code fluide-structure[END_REF]. For example, it allowed to solve problems on membranes [START_REF] Gross | New quasi-monolithic method to solve dynamic fluidstructure interaction problems on membranes[END_REF] or vortex induced vibrations on offshore risers [START_REF] Lothodé | Numerical study of viv over a flexible riser[END_REF].

Forces (stresses on the wetted surface) are given through an arbitrary polygonal surface mesh by the fluid solver. They are interpolated on the structure by projecting the force vector onto the neutral axis of the beam.

𝒇 (𝑖,𝑛) structure = 𝑃 𝑖↦𝑠 ( 𝒇 (𝑖,𝑛) interface ) (2)
In between two beam nodes, a linear interpolation is done for the torsion and deflection depending on the position of the projected forces. The displacements of the structure (at a time step 𝑖 and sub-iteration 𝑛: 𝛿𝒙 (𝑖,𝑛) structure ) are propagated to the interface surface mesh through a cubic spline representing the neutral axis of the beam with 𝑃 𝑠↦𝑖 the operator mapping the structure position to the interface position:

𝒙 (𝑖,𝑛) interface = 𝑃 𝑠↦𝑖 ( 𝒙 (𝑖,𝑛-1) structure + 𝛿𝒙 (𝑖,𝑛) structure ) (3)
Once the structure has moved the interface surface mesh, a comparison is done between the old position of the surface nodes and the new ones:

𝛿𝒙 (𝑖,𝑛) interface = 𝒙 (𝑖,𝑛) interface -𝒙 (𝑖,𝑛-1) interface (4)
If the maximum observed displacement of the interface nodes is non-zero:

‖ ‖ ‖ 𝛿𝒙 (𝑖,𝑛) interface ‖ ‖ ‖∞ > 𝜀 displacement (5)
a mesh-morphing algorithm is performed to modify the position of all the nodes in the volume in order to maintain a good quality mesh. The algorithm is described in Algorithm 1.

A mesh movement velocity 𝒖 𝑖,𝑛) mesh is computed from the new position of the mesh and used through the ALE capability of the fluid solver. update structure 𝛿𝒙 (𝑖,𝑛) structure ; deform mesh according to 𝛿𝒙 (𝑖,𝑛) interface ; compute 𝒖 (𝑖,𝑛) mesh ; while ‖ ‖ ‖

𝑟 (𝑖,𝑛) fluid ‖ ‖ ‖
< 𝜀 do 𝑛 ← 𝑛 + 1; update linearly 𝒖 (𝑖,𝑛) fluid and 𝒑 (𝑖,𝑛) fluid using 𝒖 (𝑖,𝑛-1) mesh ; compute

‖ ‖ ‖ 𝑟 (𝑖,𝑛) fluid ‖ ‖ ‖ ;
compute 𝒇 (𝑖,𝑛) fluid ; distribution of 𝒇 (𝑖,𝑛) fluid on the structure ; update structure 𝛿𝒙 (𝑖,𝑛) structure ; deform mesh according to 𝛿𝒙 (𝑖,𝑛) interface ; compute 𝒖 (𝑖,𝑛) mesh ; end end

Mesh morphing

After the interface deformation, the whole of the fluid mesh needs to be moved. This deformation occurs at each coupling iteration. The number of calls to this procedure being non-negligible, the mesh deformation needs to be fast. For this purpose, a method has been developed which propagates the deformation state to the fluid mesh. The basic principle is to propagate, from the interface to the boundaries, the information about displacements. A smoothing step can be needed to improve the quality of the cells. The details of the method are given in [START_REF] Durand | Interaction fluide-structure souple et legere[END_REF].

Fluid Validations

Description of the reference case

In [START_REF] Bahaj | Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank[END_REF], the authors describe an experiment of a small-scale tidal turbine. The tests were carried out in a cavitation tunnel at the University of Southampton (see Figure 2). The rotor radius of the turbine is 𝑟 = 400 mm. It was chosen as a compromise between maximizing the Reynolds number and reducing the tunnel blockage correction.

In the original work, the blockage correction is based on an actuator disk model of the flow through the turbine in which the flow is presumed to be uniform across any crosssection of the stream tube enclosing the turbine disc. For example, with a single rotor and a thrust coefficient of 0.8, the corrections amounted up to an 18% decrease in power coefficient and an 11% decrease in thrust coefficient.

All the curves presented in this article take into account this blockage correction. It means that all results are presented as open water curves. These numerical computations are open water computations and no correction are necessary.

The blades are made out of the NACA 63-8xx series. The distribution of pitch, thickness and chord can be found in Table 2. The pitch angle reference is the angle of the section 𝑟 = 80𝑚𝑚. It means that a pitch of 15°is simply taking the blade untouched and 20°means adding a rotation of 5°to the blade. To compare the performance of the tidal turbine, we introduce three dimensionless quantities. The thrust coefficient 𝐶 𝑇 is equal to the force in the flow direction divided by 1 2 𝜌𝑉 2 𝑆. The power coefficient 𝐶 𝑃 is equal to the torque times the rotation speed divided by 1 2 𝜌𝑉 3 𝑆. The tip speed ratio, TSR, is the ratio of blade tip speed to flow speed. These notations and definitions are reported in Table 3. 𝜌 is the density of the fluid (here 998.3 kg/m 3 ). 𝑟 is the radius of the turbine (here 0.4 m).

𝐶 𝑇 = 𝑇 1 2 𝜌𝑉 2 𝑆 thrust coefficient 𝐶 𝑃 = 𝑄 × Ω 1 2 𝜌𝑉 3 𝑆 power coefficient

Description of the simulations

The The boundary layer is captured with a 𝑦 + ≈ 30, thanks to the insertion of cell layers. Depending on the surface, 10 to 20 layers were added in order to reach the needed size, with a geometric expansion law with a factor of 1.2 increase in thickness between each pair of layers. Refinement boxes were placed to capture the tip vortices and the close wake correctly.

Finally, an orthogonality criterion was setup for the mesh optimization step (smoothing) of 20°. The mesh refinement size was chosen in order to satisfy both best practices and constraints on the allocated CPU time of the project. As shown in the next section, the mesh yields satisfactory results even if a mesh independence should be pursued in future research.

The sides and outlet boundary conditions use a zero-gradient boundary condition for pressure (𝛁𝑝 ⋅ 𝒏 = 0). The inlet boundary condition is an imposed velocity of 𝒖 fluid = 𝑉 .

With the Rotating Frame Method, a rotation mesh velocity is added to the cylinder through the ALE formulation (𝒖 mesh = Ω ∥ 𝒅 ∥ 2 with 𝒅 the distance perpendicular to the rotation axis) while the mesh if physically fixed:

𝜕𝒙 mesh 𝜕𝑡 = 0 (6)
In other words, the mesh velocity is not computed from the mesh movement, because it is fixed, but is imposed as if the mesh was rotating. With this method, the wake shows a steady behavior. The time step only plays a role for the convergence toward the steady solution and therefore can be very large.

With the dynamically rotating mesh, the mesh is rotating at an imposed velocity. The fluid velocity is not imposed, instead the mesh velocity is computed from the mesh movement :

𝒖 mesh = 𝜕𝒙 mesh 𝜕𝑡 (7)
The derivative is approximated through a Backward Differentiation Formula of order one or two. In this case, the wake is unsteady from the point of view of the observer. Both methods exploit the ALE capability of the fluid solver, but in a different manner.

When the mast is included in the domain, the Rotating Frame Method cannot be used accurately anymore, and only the dynamically rotating mesh can be used.

Validating all pitches

To validate our fluid model, we perform the five runs listed in Table 4 using the Rotating Frame Method. A time step of Δ𝑡 = 1∕20 𝑡ℎ of a rotation

( Δ𝑡 = 1 20 𝑉 ×𝑇 𝑆𝑅 𝑟×2𝜋
) is used, and iterations continue until convergence is reached. For each pitch angle, seven TSR are simulated (4, 4.5, 5, 5.5, 6, 7 and 8). Results are shown in Figure 4.

The blade pitch angle of 15°corresponds to the highest angle of attack. The flow around the blade is fully separated at a TSR of 3 and still partially detached for a TSR from 4 to 5. For TSR between 4.5 and 5.5, the error observed between the experimental data and the simulation is about 5% for the 𝐶 𝑃 and 2% for the 𝐶 𝑇 . Outside this range, the error is larger, especially for high TSR where the experimental 𝐶 𝑃 drops quicker than what is observed with the simulations. The observed 𝐶 𝑃 peak is observed for a TSR of 5, which matches the experimental peak, with a slightly larger value (0.4626 for the simulation, 0.44 for the experiment).

The blade pitch angle of 20°corresponds to the design (optimum) angle. The flow is only partially detached for a TSR of 3 on the upper part of the blade, and fully separated for the other half of the blade. Some detached flow can still be seen for TSR up to 4.5, but is limited to a small area of the blade. The error made between the experimental and the simulation results is very small (about 1% for both 𝐶 𝑃 and 𝐶 𝑇 ). It is probably due to an accurate blockage correction, and an attached flow behavior around the blades. The 𝐶 𝑃 peak is not as clear as the one with an angle of 15°, and occurs at a TSR of 5.5 (𝐶 𝑃 = 0.4533), although the value obtained for a TSR of 5 and 6 are close (0.4450 and 0.4531 respectively).

The blade pitch angle of 25°is starting to show a significant loss in term of power peak, and it is even worse for the blade pitch angles of 27°and 30°. With the angle of 25°, the flow behavior is similar to the blade pitch angle of 20°with only the lower third part of the blade showing separation. For the blade pitch angles of 27°and 30°, this separation is even smaller. For the blade pitch angle of 30°, the flow starts to separate at the tip, on the front side. The difference observed between the simulation and the experiment is only significant for the blade pitch angle of 30°, otherwise the agreement is good (less than 5% difference). For the blade pitch angles of 25°, 27°and 30°, the 𝐶 𝑃 peaks at a TSR of 5, 4 and 4 respectively, with a 𝐶 𝑃 of 0.3491, 0.3012 and 0.2391.

Comparison between RFM and dynamic rotating mesh

The goal of this section is to introduce a moving rotating mesh and to study the time step that should be used in this case. We use the same domain and mesh as in Section 3.3. In the case of a RFM computation, the best practice Numeca (2018) is to use a time step of a 1∕20 𝑡ℎ of a rotation. In the case of a dynamic rotating mesh (DRM), we tried two time steps: 1∕100 𝑡ℎ and 1∕200 𝑡ℎ of a rotation.

Performing a first set of computations using the best practices, the first results obtained were different and are shown in full and dashed lines in Figure 5. A second set of DRM computations was launched using a time step of 1∕200 𝑡ℎ of a rotation, and the results were in agreement with the RFM and with the experiment, see dotted line in Figure 5. The remaining difference is due to the unsteady behavior of the flow that is described with a better accuracy in the case of the DRM computation. In addition, the wake of the rotor outside the cylindrical inner domain is different when using DRM computations because the position of the turbine change with time whereas in the case of RFM it does not. In conclusion, the use of a time step below a 1∕200 𝑡ℎ of a rotation should be used in the case of a DRM computation. Reducing the time step size further would lead to a greater computational effort.

Fluid-Structure Computations

For all fluid-structure computations, the pitch angle used is 20°and the flow speed is 𝑉 = 1.73 m/s, as it was the design angle for this tidal turbine.

Description of the structure

In the case of the experiment from [START_REF] Bahaj | Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank[END_REF], the blades are made of solid aluminum and appears to be rigid. In order to have more flexible structures, we modeled the structure as a restricted section of the blade, using a beam inside the blade. The structure box is located at the largest part of the profile (see Figure 6).

The root of the blade is rendered rigid in our computations at a radius smaller than 0.08 m, in grey in Figure 7. For the structure labeled as 𝐸 4 , we used a Young modulus of 69 GPa, corresponding to aluminum. The Poisson's ratio is 0.35. We introduce two other Young moduli:𝐸 2 = 𝐸 4 2 and 𝐸 1 = 𝐸 4 4 . The results using the module Young labeled 𝐸 4 are referred to as 𝐸 4 . We do the same for 𝐸 2 and 𝐸 1 .

In the local frame of the beam, 𝑥 is the direction of the neutral axis, 𝑦 is orthogonal of 𝑥 and corresponds to the flapwise direction, 𝑧 is orthogonal to 𝑥 and 𝑦 and corresponds to the turbine rotation axis. The properties of the beam are given with the notations in Table 5. The obtained properties are shown in the Figure 7. As 𝐸 is constant, and the other quantities vary for each cross section, we chose to show 𝐸×𝑆, the axial stiffness, 𝐸×𝐼 𝑦 the flap wise bending stiffness, 𝐸 × 𝐼 𝑧 the rotation wise bending stiffness, and 𝐺 × 𝐽 the torsional stiffness.

E × I y E 1 = E 4 4 E 2 = E 4 2 E 4 0.1 0.2 0.3 0.4 Span (m)
E × I z E 1 = E 4 4 E 2 = E 4 2 E 4 0.1 0.2 0.3 0.4 Span (m) 0.0 0.2 0.4 0.6 0.8 1.0 E × S ×10 7 E 1 = E 4 4 E 2 = E 4 2 E 4 0.1 0.2 0.3 0.4 Span (m) 0 20 40 60 80 G × J E 1 = E 4 4 E 2 = E 4 2 E 4

Quasi-static computations

The performance curves were obtained for different structures using quasi-static computations. A quasi-static computation consists of running a long steady fluid computation, then after each convergence of the fluid forces, those forces are applied to the structure and the new structural position is computed using a steady state scheme. After a few iterations, the forces between the structure and the fluid are in equilibrium. The obtained result is a steady state.

In Figure 8, we can see that the results between the fixed computation and 𝐸 4 are close. The 𝐶 𝑇 is a little bit superior, but the 𝐶 𝑃 is similar. 𝐸 2 and 𝐸 4 are also showing an increasing 𝐶 𝑇 , but there is a decay for high TSR, because the blade aligns itself with the flow. 𝐶 𝑃 performances deteriorates significantly with decreasing stiffness, for the same reason.

Figure 9 shows the deformed blade for the different Young moduli at TSR 5. 𝐸 1 has the biggest deflection, and is quite unrealistic. The deflection occurs both in the direction of flow and in the direction of rotation.

Blade-mast dynamic interaction with FSI

A mast is added to the computation domain. The mast is considered rigid, and the hub split into two parts. The first part is fixed and linked to the mast, and the second part is rotating with the blades. As previously, the domain is split into two domains, but the fixed domain now includes the mast and a portion of the hub now. An outline of the mesh is shown in Figure 10. This mesh consists of 7.9 million cells.

We chose to consider 𝐸 4 and 𝐸 2 only, because 𝐸 1 showed poor performance results and was therefore considered unrealistic (see Section 4.2). We used a time step of 10 -3 s, which is lower than a 1 200 𝑡ℎ of a rotation period. The water turbine is rotating at 21.625 rad/s which corresponds to TSR 5. The frequency of the blade passage in front of the mast is 3.44 Hz. We record the position and forces on all structural nodes of the blade. Acceleration at the tip of the blade is shown in Figure 11. The start of the graph shows the recovery after passing in front of the mast. At 𝑡 = 3.78 s, the blade arrives near the mast. The acceleration reaches 30 m/s 2 for 𝐸 2 and 15 m/s 2 for 𝐸 4 . The recovery starts with an acceleration peaking at 8 m/s 2 for 𝐸 2 and 6 m/s 2 for 𝐸 4 at 𝑡 = 3.82 s for 𝐸 2 and 3.81 s for 𝐸 4 . The time for recovery (acceleration of 0 m/s 2 ) is approximately the same and take 0.17 s. A cycle takes 0.29 s at TSR 5. The first acceleration and its recovery takes ≈ 60% of a cycle. For larger TSR, this would be even more noticeable.

Figure 12 shows the linear force compared to the average force on different sections of the blade. Upper and middle sections show a similar behavior: there is a spike in relative negative force when going past the mast and a recovery afterward. The bottom portion is much more chaotic which is due to vortex shedding. The upper blade portion shows some difference between 𝐸 2 and 𝐸 4 while for the middle portion this difference is less noticeable. For the upper blade portion, the difference between 𝐸 2 and 𝐸 4 is 2.5 N/m at 𝑡 = 3.75 s and is almost 0 N/m during recovery.

Figure 13 shows the absolute forces acting on each structural node of the blade. The loading is the greatest at the center of the beam (23 N at 𝑟 = 0.25 m) and the lowest at the Figure 14 shows the difference between the fluid force and overall time averaged force for each structural node of the blade over time. With this post-processing, it is easier to see the transient behavior. For blade radii greater than 𝑟 = 0.1 m, there is a negative spike in forces just at the passage of the mast. It is then recovered with a more continuous positive pressure. The negative spike occurs on the whole blade, especially for 𝐸 4 , the upper part for 𝐸 2 being less spiky. There are very fluctuating forces at the root of the blade.

Figure 15 shows the Fourier transform of the forces in the frequency space. for 𝐸 2 and 𝐸 4 the frequency corresponding to the rotation frequency are highly excited as expected. The second mode (≈ 7 Hz) is higher for 𝐸 2 than for 𝐸 4 . It is the contrary for the third mode (≈ 10.3 Hz). The fourth mode (≈ 14 Hz) is existent in 𝐸 2 but barely present for 𝐸 4 . At the root of the blade, the frequency we can see are related to the vortices advected in the flow.

Figure 16 shows the acceleration of each structural node of the blade. The acceleration is the greatest at the tip of the blade, as one should expect. As seen before, the acceleration peak is half the magnitude for 𝐸 4 compared to 𝐸 2 . The most interesting feature is the positive acceleration after the peak, called here recovery, which is shorter in duration when the stiffness is higher. For 𝐸 2 , the acceleration is non-negligible for a large portion of the cycle, while for 𝐸 4 it concerns less than half of the cycle.

Figure 17 shows the Fourier transform of the acceleration for each structural node of the blade. It is scaled radius by radius to show which mode is excited for each configuration. For 𝐸 2 , the first mode -which is the frequency of rotation 3.44 Hz -is highly excited, while for 𝐸 4 it is less the case. The second mode (≈ 7 Hz) is the principal mode to be excited for 𝐸 4 and is also excited for 𝐸 2 . The third mode (≈ 10 Hz) is highly excited for 𝐸 2 . The fourth mode (≈ 14 Hz) is much more excited for 𝐸 4 than 𝐸 2 . There is not a huge difference in spectrum between the tip and the root except for very high frequencies.

The Figure 18 is showing the vortices in the domain using the Q criterion (TSR 5 and 𝐸 2 ). A zoom at the root of the blade shows a vortex due to the cylindrical shape at this section. As stated before, the vortex oscillation at the root of the blade has a different frequency than the rotation frequency.

Figure 19 shows a cut of the mesh showing the deformed mesh. The sliding interface is treated as non-deformable mesh morphing algorithm.

Conclusion

A dynamic methodology has been developed to perform fluid computations and was validated against results from [START_REF] Batten | Experimentally validated numerical method for the hydrodynamic design of horizontal axis tidal turbines[END_REF]. The results show a good agreement with the experiment especially for the design pitch angle of 20°w here the error is minimal. In order to use a dynamic rotating mesh, a criterion for the time stepping rule was found. Three different structures were setup and simulated using quasi-static fluid-structure computations. The difference between a flexible and a rigid blade was most evident at a TSR greater than 5. These cases excluded the most flexible blade (𝐸 1 ) for the rest of the study because it produced a poor power curve and thus did not seem suitable for production.

𝐸 2 and 𝐸 4 were then setup with a mast. The rotation of the blades in front of the mast induced vibrations. The vibrations and forces were recorded for comparison. Dif-ferent behaviors were highlighted. The most flexible of the two structures (𝐸 2 ) showed the greatest displacements and accelerations but the lowest vibrational frequencies.

For future work, wave induced vibrations could be introduced. There are already some experimental work done in this area, for example [START_REF] Luznik | The effect of surface waves on the performance characteristics of a model tidal turbine[END_REF]; [START_REF] Gaurier | Flume tank characterization of marine current turbine blade behaviour under current and wave loading[END_REF]. An optimization based on fatigue analysis could help to maximize the durability of such a turbine. For example, the recent work of [START_REF] Finnegan | Operational fatigue loading on tidal turbine blades using computational fluid dynamics[END_REF] shows the interest of fatigue loading for tidal turbines. It would also be very interesting to compare against experimental data.

Figure 1 :

 1 Figure 1: Notations used for the fluid-structure domain.

Algorithm 1 :

 1 Quasi-monolithic algorithm. while 𝑡 𝑖 < 𝑡 final do 𝑖 ← 𝑖 + 1; 𝑡 𝑖 ← 𝑡 𝑖-1 + Δ𝑡; 𝑛 ← 0; distribution of 𝒇 (𝑖-1,𝑛 final ) fluid on the structure ;

Figure 2 :

 2 Figure 2: Photo of the experiment (from Bahaj et al. (2007)).

  fluid domain is decomposed in two parts. The outer domain is made of a box of width 𝐿 𝑥 = 4 m, length 𝐿 𝑦 = 6 m and height 𝐿 𝑧 = 4 m. A smaller domain, cylindrical, represents the domain in rotation and has a length of 𝐿 ′ 𝑦 = 0.6 m and a radius of 𝑟 ′ = 0.6 m. It includes the hub part that is rotating and the blades. The two domains are shown in Figure 3. The inner domain can rotate inside the outer domain. The two domains are linked through a sliding interface. The mesh is generated with Hexpress™, part of FINE™/Marine, with an octree method and consists of 4.1 million hexagonal cells.

Figure 3 :

 3 Figure 3: Different views of the fluid mesh. From left to right, top to bottom: computational domain, internal cylindrical domain (used as a sliding interface), tidal surface mesh, volume mesh cut.

Figure 5 :

 5 Figure 5: Comparison between RFM and DRM. Full line, dashed line and dotted line represents results obtained respectively with the RFM method, with a rotating mesh with a time step of Δ𝑡 = 𝑐𝑦𝑐𝑙𝑒 100 and with Δ𝑡 = 𝑐𝑦𝑐𝑙𝑒 200 . Markers are experimental results from Bahaj et al. (2007).

Figure 6 :

 6 Figure 6: The box reinforced in blue for section at 𝑟 = 80𝑚𝑚.

Figure 7 :

 7 Figure 7: Structural properties of the beam along the span of the blade.

Figure 8 :

 8 Figure 8: Performance curves for different structural properties.

Figure 9 :

 9 Figure 9: Different views of deformed blades (grey: fixed, blue:𝐸 4 , green:𝐸 2 , red:𝐸 1 ). The left figure is from the turbine side, the middle figure shows the rear view, and the right figure shows the top view.

Figure 10 :

 10 Figure 10: Different views of the initial fluid mesh for the blade-mast dynamic interaction. The figure on the left shows the sliding interface with the mesh rotating in the inner cylinder and the rest is remaining fixed. The figure on the right shows a slice of the mesh viewed from the top.

Figure 11 :

 11 Figure 11: Acceleration of the tip of the blade (in m/s 2 ).

Figure 12 :

 12 Figure 12: Linear force compared to average linear force on different sections of the blade (in N/m).

Figure 13 :

 13 Figure 13: Linear force on the blade according to time (in N/m).

Figure 14 :

 14 Figure 14: Linear force compared to average forces acting on the blade according to time (in N/m).

Figure 15 :

 15 Figure 15: Frequency analysis of forces acting on the blade.

Figure 16 :

 16 Figure 16: Acceleration of each structural node of the blade according to time (in m/s 2 ).

Figure 17 :

 17 Figure 17: Frequency analysis of acceleration of each blade element blade.

Figure 18 :

 18 Figure 18: Instant view of the simulation showing the Q criterion > 100.

Figure 19 :

 19 Figure 19: View of the deformed mesh along the blade.

Table 1

 1 Notations for the fluid domain, structure domain and interface (see Figure1).

	Symbol	Description
	Ω	complete domain Ω = Ω 𝑓 ∪ Ω 𝑠
	Ω 𝑓	fluid domain
	Ω 𝑠	structure domain
	Σ	fluid-structure interface, Σ = Ω 𝑓 ∩ Ω 𝑠
	Γ 𝑓	fluid boundaries
	𝒏	normal outer vector

Table 2

 2 Blade geometry description.

	r (𝑚𝑚)	c/R	Pitch (°)	t/c (%)
	80	0.125	15	24
	100	0.1203	12.1	22.5
	120	0.1156	9.5	20.7
	140	0.1109	7.6	19.5
	160	0.1063	6.1	18.7
	180	0.1016	4.9	18.1
	200	0.0969	3.9	17.6
	220	0.0922	3.1	17.1
	240	0.0875	2.4	16.6
	260	0.0828	1.9	16.1
	280	0.0781	1.5	15.6
	300	0.0734	1.2	15.1
	320	0.0688	0.9	14
	340	0.0641	0.6	14.6
	360	0.0594	0.4	13.6
	380	0.0547	0.2	13.1
	400	0.05	0	12.6

Table 3

 3 Notations used for the dimension and dimensionless numbers.

	Name		Property
	𝑟		radius of the turbine (m)
	𝑆 = 𝜋𝑟 2		swept area (m 2 )
	𝑉		inlet velocity (m/s)
	𝜌		density of the fluid (kg/m 3 )
	𝑄		turbine torque (Nm)
	𝑇		turbine drag (N)
	Ω		rotation speed (rad/s)
	𝑇 𝑆𝑅 =	Ω𝑟 𝑉	tip speed ratio

Table 4

 4 Pitch angle and corresponding flow speed. Power and thrust coefficients for all cases in Table4with the Rotating Frame Method. Simulations are displayed with solid circular markers connected by a dotted line and experiment data are displayed with empty circular markers.

	Pitch angle (°)			Flow speed
	15°1.40 m/s
	20°1.73 m/s
	25°1.54 m/s
	27°1.30 m/s
	30°1.54 m/s
		0.5
		0.4
		0.3
	CP
		0.2
		0.1					exp 15°s
							im 15°e
							xp 20°s
		0.0					im 20°e
							xp 25°s
		1.4					im 25°e xp 27°s
		1.2					im 27°e xp 30°s
		1.0					im 30°F
	CT	0.8
		0.6
		0.4
		0.2
		0.0	4	5	6	7	8
					TSR
	igure 4:

Table 5

 5 Notation for the structure properties.

	Symbol	Description
	𝐸	Young modulus (Pa)
	𝜈	Poisson's ratio
	𝐺 = 𝐸∕(2(1+𝜈) shear modulus (Pa)
	𝑆	section surface of the structure (m 2 )
	𝐼 𝑦	second moment of inertia for bending
		flap wise (around 𝑦, m 4 )
	𝐼 𝑧	second moment of inertia for bending
		rotation wise (around 𝑧, m 4 )
	𝐽	polar moment of inertia (around 𝑥, m 4 )

Acknowledgment

This work was initially funded by K-Epsilon and the ANRT via a CIFRE doctoral grant and a partnership with the LMN of INSA Rouen. K-Epsilon supported the project through a partnership with LMRS -CNRS. We are thankful for the computing support provided by CRIANN. We would like to thank NUMECA International, now Cadence, for the R&D license.

Marine Renewable Energy (MRE) lmn (D. Lemosse); / (D. Gross);


Guy Archambault

Stéphane Poirier

Alain Rouleau

Sylvie Gentier

Joëlle Riss

The behavior of induced pore fluid pressure in undrained triaxial shear tests on fractured porous analog rock material specimens

à la diffusion de documents scientifiques 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.

The behavior of induced pore fluid pressure in undrained triaxial shear tests on fractured porous analog rock material specimens.

Mechanics of Jointed and Faulted Rock (MJFR-3) 1998, Rossmanith ed .. A.A. Balkema, pp.555-560.