Mihai Arghir
email: mihai.arghir@univ-poitiers.fr

Omar Benchekroun

A simplified structural model of bump-type foil bearings based on contact mechanics including gaps and friction

Keywords: foil bearing structure, contact mechanics, friction

) was implemented for solving the structural problem with friction and close/loose gaps.

The procedure uses the augmented Lagrangian algorithm for the normal contact forces and the penalty method for the tangential friction force. The procedure is very robust in dealing with stickslip situations and close/loose contacts. The model was validated by comparisons with results from the literature ([31], [32]).

By considering gaps in the foil structure and close/loose contacts, the model is not only closer to reality but can also handle manufacturing errors of the foil structure. This is shown by considering bump foil height manufacturing errors.

Introduction

Aerodynamic foil journal bearings represent a technology under scrutiny since decades. They represent a viable, proven technology for small and medium size, high speed rotating machinery operating on air or with gases. The key point of foil bearings is the compliant foil structure that brings a considerable damping compared to the very thin air film during operation. There are many types of foil bearings or, more generally, bearings with a compliant structure. The first bearing of this type, still in use, is the foil leaf bearing. It consists of many beam-like, thin foils that slightly overlap. These kind of bearings have a low ratio of load capacity versus damping. The bump-type foil bearings are more elaborate and use two foils (Figure 1): the smooth top foil and the underneath corrugated foil provided with bumps that act as cross-coupling springs. This kind of foil bearings have better characteristics in terms of load capacity and damping. Another kind of foil bearings use beam-like springs instead of bumps. However, these solutions are less popular. Among these types, the bump-type foil bearings are the most widely used and analyzed. The difficulties in accurately describing the characteristics of the foil structure are the reason behind systematic research efforts developed for years. An overview of the large number of publications dealing with this subject can be found in the review article [START_REF] Branagan | Compliant Gas Foil Bearings and Analysis Tools[END_REF].

The first model [START_REF] Heshmat | Analysis of Gas-Lubricated Compliant Thrust Bearings[END_REF], [START_REF] Heshmat | Analysis of Gas-Lubricated Foil Journal Bearings[END_REF] ignored the stiffness of the top foil and considered the bumps as noninteracting springs. By considering the stiffness as being uniformly distributed, the foil structure appeared as a Winkler foundation. Friction was subsidiary taken into account [START_REF] Iordanoff | Analysis of an Aerodynamic Compliant Foil Thrust Bearing: Method for a Rapid Design[END_REF] but the model could not predict the damping of the compliant structure and therefore additional viscous damping [START_REF] Peng | Calculation of Stiffness and Damping Coefficients for Elastically Supported Gas Foil Bearings[END_REF] or structural loss [START_REF] Rubio | Structural Stiffness, Dry Friction Coefficient, and Equivalent Viscous Damping in a Bump-Type Foil Gas Bearing[END_REF] parameters had to be added. This structural model became widespread used due to its easy implementation ( [START_REF] Kim | Effect of Side Feed Pressurization on the Dynamic Performance of Gas Foil Bearings: A Model Anchored to Test Data[END_REF], [START_REF] Ryu | Prediction of Axial and Circumferential Flow Conditions in a High Temperature Foil Bearing With Axial Cooling Flow[END_REF], [START_REF] Rubio | Structural Stiffness, Dry Friction Coefficient, and Equivalent Viscous Damping in a Bump-Type Foil Gas Bearing[END_REF], [START_REF] Bou-Said | On Nonlinear Rotor Dynamic Effects of Aerodynamic Bearings With Simple Flexible Rotors[END_REF], [START_REF] Bonello | The efficient computation of the nonlinear dynamic response of a foilair bearing rotor system[END_REF]). The model was extended by taking into account the top foil.

Carpino et al. [START_REF] Carpino | Effects of Membrane Stresses in the Prediction of Foil Bearing Performance[END_REF], [START_REF] Peng | Finite Element Approach to the Prediction of Foil Bearing Rotor Dynamic Coefficients[END_REF] introduced a membrane and bending effect of the top foil while the bump foil was still modeled as an elastic foundation. The elastic foundation model was extended in [START_REF] Carpino | A Fully Coupled Finite Element Formulation for Elastically Supported Foil Journal Bearings[END_REF] for taking into account radial and circumferential displacements of bumps. Other references modeled the top foil either as a beam or a shell ( [START_REF] Bruckner | Simulation and Modeling of the Hydrodynamic, Thermal and Structural Behavior of Foil Thrust Bearings[END_REF], [START_REF] San Andrés | Analysis of gas foil bearings integrating FE top foil models[END_REF], [START_REF] Kim | Hydrostatic Air Foil Bearings: Analytical and Experimental Investigations[END_REF], [START_REF] Nielsen | Transient and steady state behavior of elasto-aerodynamic air foil bearings, considering bump foil compliance and top foil inertia and flexibility: A numerical investigation[END_REF], [START_REF] Leister | Computational Analysis of Foil Air Journal Bearings Using a Runtime-efficient Segmented Foil Model[END_REF]) in contact with isolated bumps. These latter are described by the stiffness analytic formulas from [START_REF] Iordanoff | Analysis of an Aerodynamic Compliant Foil Thrust Bearing: Method for a Rapid Design[END_REF]. The models have the advantage of predicting the sagging effects of the top foil. Moreover, the model in [START_REF] Nielsen | Transient and steady state behavior of elasto-aerodynamic air foil bearings, considering bump foil compliance and top foil inertia and flexibility: A numerical investigation[END_REF] considers also the inertia of the top foil. However, their accuracy can be questioned because they still discard any direct interaction between bumps (a.e. the bending moment that links two successive bumps). Ku and Heshmat [START_REF] Ku | Compliant Foil Bearing Structural Stiffness Analysis-Part 1: Theoretical Model Including Strip and Variable Bump Foil Geometry[END_REF] have partially remedied this limitation by adding plate elements between successive bumps and by bookkeeping the different stick-slip state of bumps under load and friction forces. However, this algorithm is quite difficult to implement.

A different category of models discretizes both the top and the bump foil with structural finite elements. Lee et al. [START_REF] Lee | The static performance analysis of foil journal bearings considering three-dimensional shape of the foil structure[END_REF], [START_REF] Lee | The Dynamic Performance Analysis of Foil Journal Bearings Considering Coulomb Friction: Rotating Unbalance Response[END_REF], [START_REF] Lee | The Effect of Coulomb Friction on the Static Performance of Foil Journal Bearings[END_REF] modeled the foils with bar and beam elements and took friction into account by using a stick-slip algorithm similar to the one employed in [START_REF] Ku | Compliant Foil Bearing Structural Stiffness Analysis-Part 1: Theoretical Model Including Strip and Variable Bump Foil Geometry[END_REF]. Lehn et al. [START_REF] Lehn | Elasto-gasdynamic modeling of air foil thrust bearings with a two-dimensional shell model for top and bump foil[END_REF] used two dimensional shell elements and modeled Coulomb forces with a penalty contact model and a regularized friction law.

A finite element model with a Lagrange multiplier contact model was employed in [START_REF] Li | Modeling of a multi-layer foil gas thrust bearing and its load carrying mechanism study[END_REF] for foil leaf type bearings but friction was discarded.

Barzem et al. [START_REF] Barzem | Aero-elastic bearing effects on rotor dynamics: a numerical analysis[END_REF] introduced a non-linear model using shell elements for the bumps and for the top foil. The model took into account large (non-linear) deformations of the bumps including buckling.

Friction was taken into account following an updated Lagrangian algorithm from Kalker [START_REF] Kalker | Three Dimensional Elastic Bodies in Rolling Contact[END_REF].

Larsen et al. [START_REF] Larsen | Numerical and experimental investigation of bump foil mechanical behavior[END_REF] adopted a similar approach based on shell elements and large displacements for the bump and the top foil. However, friction was taken into account by using a very specific approach related to penalty method.

Other researchers employed non-linear elasticity commercial codes for modeling the foil structure ( [START_REF] Zywica | The static performance analysis of the foil bearing structure[END_REF], [START_REF] Temis | Gas-dynamics Foil Bearing Model[END_REF], [START_REF] Liu | Simulation of Compliant Bump Foil Journal Bearing Using Coupled Reynolds Equation and Finite Element Model Method[END_REF], [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF]). These approaches benefited from the implementation of modern, full nonlinear contact and friction algorithms but the models were somewhat difficult to handle due to the large computational effort. Le Lez and Arghir introduced a different kind of model in 2007 [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF]. The model neglected the stiffness of the top foil and modeled the bump foil as a network of interacting trusses behaving like springs. This lead to a linear stiffness matrix of the bump foil defined in a manner similar to the finite element method where each bump is an element. The model was validated by comparison against full non-linear simulations performed with a commercial code [START_REF] Lez | Static and Dynamic Characterization of a Bump-Type Foil Bearing Structure[END_REF] and was also coupled to the gas film [START_REF] Lez | A Dynamic Model for Dissipative Structures Used in Bump-Type Foil Bearings[END_REF] in journal bearings. Coulomb forces were taken into account by using a dynamic model from [START_REF] Petrov | Generic Friction Models for Time-Domain Vibration Analysis of Bladed Disks[END_REF] where regularized friction forces were estimated from first order differential equations.

Structural models of the bump foil similar to [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] were successfully developed by Hryniewicz et al. [START_REF] Hryniewicz | Structural Properties of Foil Bearings : A Closed Form Solution Validated with Finite Element Analysis[END_REF], Feng and Kaneko [START_REF] Feng | Analytical Model of Bump-Type Foil Bearings Using a Link-Spring Structure and a Finite Element Shell Model[END_REF], Gad and Kaneko [START_REF] Gad | A New Structural Stiffness Model for Bump-Type Foil Bearings: Application to Generation II Gas Lubricated Foil Thrust Bearing[END_REF], von Osmanski [START_REF] Osmanski | A fully coupled air foil bearing model considering friction -Theory & experiment[END_REF] or by Hoffman et al. [START_REF] Hoffmann | A valid method of gas foil bearing parameter estimation: A model anchored on experimental data[END_REF]. Some of these models took the top foil into account. The top foil was modeled either as a shell [START_REF] Feng | Analytical Model of Bump-Type Foil Bearings Using a Link-Spring Structure and a Finite Element Shell Model[END_REF], [START_REF] Gad | A New Structural Stiffness Model for Bump-Type Foil Bearings: Application to Generation II Gas Lubricated Foil Thrust Bearing[END_REF] or as a beam [START_REF] Osmanski | A fully coupled air foil bearing model considering friction -Theory & experiment[END_REF]. Friction was taken into account using either stick-slip bookkeeping algorithms ( [START_REF] Hryniewicz | Structural Properties of Foil Bearings : A Closed Form Solution Validated with Finite Element Analysis[END_REF], [START_REF] Feng | Analytical Model of Bump-Type Foil Bearings Using a Link-Spring Structure and a Finite Element Shell Model[END_REF], [START_REF] Gad | A New Structural Stiffness Model for Bump-Type Foil Bearings: Application to Generation II Gas Lubricated Foil Thrust Bearing[END_REF]) or dynamic, regularized friction forces ( [START_REF] Osmanski | A fully coupled air foil bearing model considering friction -Theory & experiment[END_REF], [START_REF] Hoffmann | A valid method of gas foil bearing parameter estimation: A model anchored on experimental data[END_REF]). In [START_REF] Osmanski | A fully coupled air foil bearing model considering friction -Theory & experiment[END_REF] the inertia of the top and bump foil was also considered in the dynamic model.

The model introduced in [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] is considered as a compromise between accuracy and computational effort because the structural model remains linear while the bumps are fully coupled. However, the model supposes that the contacts between the top and the bump foils as well as the contacts between the bump foil and the sleeve are always closed. This is an ad-hoc assumption aimed to simplify the numerical algorithm but calculations performed in [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF] with a commercial, non-linear elasticity code showed that this not true. Even for very simple loading cases, the contact state might change between loose or close. The model introduced in [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] could not handle these situations because it neglected the top foil stiffness and any possible gaps in the foil structure.

The purpose of the present work is to extend the model introduced in [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] for taking into account close/loose contacts between the top and the bump foil, between the bump foil and the sleeve and between the top foil and the rotor. The top foil stiffness is now taken into account and an additional degree of freedom is added to each bump.

Normal and friction forces between the foils and between the bump foil and the sleeve are solved with the augmented Lagrangian method (Uzawa algorithm) [START_REF] Wriggers | Computational Contact Mechanics[END_REF], [START_REF] Zavarise | A Guide for Engineers to Computational Contact Mechanics[END_REF]. The contact between the rotor and the top foil is solved with the penalty method.

Comparisons with the models from [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] (linear model) and from [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF] (full non-linear) for simple static loadings show hysteresis loops in good agreement. Manufacturing errors of the bump foil are then added to the new model and the results are again compared with [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF] showing the same level of agreement. The contact algorithm proved to be very robust and efficient in terms of computational time even for irregular foil structures.

Figure 1. Foil bearing (first generation design)

The foil bearing structure

The gaps

The structural model presented in [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] was based on a linear approximation of the bump foil with structure made of trusses. Curvature effects were neglected and the foil structure was unwrapped.

Figure 2 depicts two typical bumps and the associated simplified model. Each bump had three degrees of freedom (DOF), a vertical and a horizontal displacement for the top point of the bump and a horizontal displacement for the point connecting two successive bumps. It was supposed that the bump foil is continuously in contact with the top foil and with the bearing casing. Moreover, it was supposed that the stiffness of the top foil was negligible compared to the bump foil. Therefore, the structural model neglected the top foil. Closed and open contacts between the bump foil and the top foil and between the bump foil the bearing casing were not taken into account. All contacts were closed and the friction forces were calculated using Coulomb law. If the contacts may become loose then additional DOFs must be taken into account as shown in Figure 3. For the contacts between the bump and the top foil, an additional DOF, the radial displacement of the top foil, must be added per bump. The gap between the top foil and the bump foil is:

, , , , (1)
with i the bump number and , the initial gaps between the bump and the top foil. The noninterference condition between the bump and the top foil is, , 0, 1 … .

The gap between the rotor and the top foil can also be described with the aid of the radial displacement of the top foil:

, , 0 (2)
were is the distance between the rotor and the undeformed top foil.

1 6 !" # ! $" (3)
The non-interference condition is, , 0, 1 … .

An additional DOF must be added to the bump foil if its contacts with the sleeve may become loose.

The additional DOF corresponds to the radial displacement of the bottom node of the bump (Figure 3). The gap between the bump foil and the sleeve is:

, , , (4)
where , are the initial gaps between the bump foil and the sleeve. The non-interference condition between the bump and sleeve is, , 0, 1 … .

Due to the additional DOFs a stiffness matrix for the top foil must be introduced and the stiffness matrix of the corrugated bump foil must be reviewed.

The stiffness of the top foil

The top foil is modeled as a circular (curved) beam with the appropriate plate correction. The beam end corresponding to the welding is fixed while the other end is free (Figure 1). The discretization nodes correspond to the top of the bumps and each node has a one DOF represented on the unwrapped foil on Figure 3 as a vertical displacement, , (downwards positive). For the circular beam , corresponds to a radial displacement (outwards positive). The elasticity matrix is obtained with the aid of Bresse formulas and yields:

! ,% & ' () ' 12 ⁄ 1 , - . / 0 / 1 !" ! $" % 2 " % !2" " % 3 2 , 4 " % 5 " !" % ! $" 2 " !2" % " 3 2 , 67) 8 !)
" , " % 9 0, , :

The stiffness matrix of the top foil is then:

; < = (6)

Figure 4. Discretization of the bump foil with 4DOF per bump

The stiffness of the bump foil

The stiffness of the bump foil with 3 DOF per bump was presented in [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF]. For the bump foil with 4DOF per bump, the parameters k1…k4 and θd depicted in Figure 4 remain the same as in [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] but an additional DOF, namely the vertical displacement of the bump bottom, , must be taken into account. For example, the potential energy of the corrugated foil with two bumps depicted in Figure 4 is:

Π 1 2 ?@ 2ΔB ΔB ' ΔB C ΔB D 3 @ 2ΔB E ΔB F 3 @ ' ΔB G @ E ΔB H (7) ΔB ΔB G ' ΔB I ! ΔB C ' I ' ! I ! ΔB ' I ! I ! ΔB D E ' I ' ! I E ! ΔB E ΔB F E
where !" J , ! ! $" J and its DOF are the elements of the following vector:

K LM N , , ,' ,E ,G ,C ,D ,F O P (8)
The stiffness matrix is obtained by deriving the potential energy Π following each DOF. This yields:

(9) ; LM Q R R R R R R R S 2@ @ ' @ E 0 @ @ ! @ ' 0 0 0 0 2@ ! @ ! @ ! 0 0 0 0 @ @ ! 2@ 2@ 0 @ @ ! @ 0 @ ! @ ! 0 2@ ! @ ! @ ! 0 0 @ ' 0 @ @ ! @ ' 2@ 0 @ @ ! 0 0 @ ! @ ! 0 2@ ! @ ! @ ! 0 0 @ 0 @ @ ! @ @ @ ! 0 0 0 0 @ ! @ ! @ ! @ ! T U U U U U U U V
This matrix can be cast as the sum of individual stiffness matrixes of each bump plus a coupling matrix (Appendix 1).

The contact model

The mathematical model of the aerodynamic foil bearing structure must include the rotor, the top, the bump foil and the bearing casing plus the non-interference constraints that link these four elements. The non-interference constraints are expressed with the aid of the gap functions described by ( 1)-( 4). If the gap functions are positive, then contacts are loose and the normal and the tangential contact forces are zero. On the other hand, if the gap functions are nil, then the contacts are closed and the normal forces are non-zero. This is mathematically expressed by the following equations:

0, W X, 5 0, W X, 0 (10)
where stands for any of the gap functions and W X, for the corresponding normal contact forces. In contact mechanics, these are known as the Hertz-Signorini-Moreau conditions [START_REF] Wriggers | Computational Contact Mechanics[END_REF].

The elastic problem of the foil structure is solved by minimizing the potential energy of the elastic structure (i.e. a functional) Π,

Π 1 2 K M P ; M K M 1 2 K < P ; < K < (11)
with constraints issued from contact conditions, 0.

The methods usually adopted in contact mechanics include the effect of constraints into the potential energy of the elastic structure. The penalty method, the Lagrange multiplier method or the augmented Lagrange method can be employed. Each method has its advantages and drawbacks. For the penalty method the extended functional of the potential energy is:

Π Y Π 1 2 Z Y [ P [ (12)
where g in eq. ( 12) are the constraints for the contacts supposed to be closed (contacts supposed loose are not included) and Z Y is a parameter similar to a spring stiffness. For high accuracy Z must be of the order 10 12 N/m or larger. The method has the advantage of being simple because the stiffness matrix obtained when minimizing Π Y has the same size as the linear elastic stiffness matrix.

However, it includes the parameters Z that might damage its conditioning.

For the Lagrange multipliers method the extended functional of the potential energy is:

Π \] Π ^P[ (13)
where the Lagrange multiplier _ is part of the solution and represents the normal contact force. As in the previous method, the Lagrange multipliers are added only for the contacts supposed to be closed. The stiffness matrix obtained by minimizing Π \] is well conditioned but its size is not a priori known (it must include the Lagrange multipliers and therefore depends on the number of contacts supposed to be closed). Compared to the linear stiffness matrix of the elastic structure without constraints, the extended matrix it is not positive definite any more.

The augmented Lagrange method avoids the ill conditioning problems of the penalty method and the problems related to the variable size of then extended matrix of the Lagrange multipliers method.

The extended functional of the potential energy is :

Π `\] Π ^P[ 1 2 Z `\] [ P [ (14)
and the Lagrange multipliers _ are iteratively updated following Uzawa algorithm [START_REF] Zavarise | A Guide for Engineers to Computational Contact Mechanics[END_REF]:

^Xab ^cdJ Z `\] [ Xab (15)
The iterative algorithm can start with _ 0 and with a limited value of the parameter Z `\] for preserving the conditioning of the resulting stiffness matrix. The bar notation _ denotes that the Lagrange multipliers are kept constant and are updated only within what is addressed as "outer" iterations loop. Their significance is the same, i.e. normal contact forces.

If friction forces are discarded, the advantages and drawbacks of these three methods were mentioned above. However, when friction forces are present the advantage of the augmented Lagrange method becomes clear.

The Coulomb friction force The friction forces are described by the Coulomb law. In its natural form this law is written in terms of the sliding velocity, e f . gW h g 5 4|W X | so that j 4 gW h g k 4|W X | 67)$ e f 0, !6 @ 67) 8 !) gW h g 4|W X | l$m ∃o 9 0 ! 67l6 e f oW h , !B p

The friction cone graphically depicts Coulomb law (Figure 5). The law is not conservative; therefore, the potential energy of the contact problem with friction has no minimum anymore. The functional (i.e. the potential energy) associated to Signorini's problem with Coulomb friction is nonconvex, nonquadratic and nondifferentiable. Thus, standard uniqueness and existence results from convex analysis cannot be applied. However, the solution has been shown to exist for sufficiently small friction coefficients [START_REF] Wriggers | Computational Contact Mechanics[END_REF].

A fixed point iteration method in the normal force is used for the numerical solution of this problem.

A contact problem with the so called Tresca friction will be solved at every iteration step. The Tresca friction law is a reduced frictional problem with a prescribed normal force. It is expressed in the same way as the Coulomb law gW h g 5 q so that j 4 gW h g k q 67)$ e f 0, !6 @ 67) 8 !) gW h g q l$m ∃o 9 0 ! 67l6 e f oq, !B p

where q 4|W r X | and W r X is the prescribed normal force. The bar in the notation of the normal force indicates it remains constant during the iterative algorithm for the calculation of the friction force W h .

This will be addressed as the "inner" iterations loop. The bar notation also indicates the link between the prescribed nodal force and the Lagrange multipliers in the augmented Lagrange algorithm. In fact W r X will be exactly the Lagrange multipliers _ updated in an « outer » iterations loop. 2In the "inner" iterations loop, the friction forces are calculated depending on the stick or slip state. It is supposed that closed contacts are initially in stick state and the friction force is described by a penalty approach.

W h sd Z h t f u f v (18)
where f u is the tangential displacement3 at iteration (k) and f is the tangential displacement at the previous time or loading step. The corresponding potential energy is:

1 2 Z h t f u f v (19)
and is added to the total potential energy of the structure. The displacements of the structure are calculated by minimizing the total potential energy in a step that will be further detailed. After updating the displacements of the structure, u , the trial value of the friction force ( 18) can be estimated and the algorithm checks for stick or slip state:

w 4 xW h sd x k 4|W r X | 67)$ W h u W h sd !6 @ !6l6) 67) 8 !) W h u 4W r X ! $ t e f u v !B p !6l6) (20)
where it should be reminded that following the Hertz-Signorini-Moreau conditions, W r X k 0. It is to be underlined that for the stick state, f u and f are very close but not exactly equal. Following the penalty approach, the largest the value of the "stiffness" Z h , the lowest the error y f u f y.

In the next iteration, nodes that are in a stick state are continued to be dealt with the penalty method and their contribution to the total potential energy is given by eq. ( 19) but the contribution of nodes that are in slip state is described by the mechanical work of the friction force:

W h u t f u f v 4W r X ! $ t e f u v t f u f v (21)
The displacements of the structure are calculated again by minimizing the total potential energy that includes the mechanical work of friction forces in slip state [START_REF] Lee | The Dynamic Performance Analysis of Foil Journal Bearings Considering Coulomb Friction: Rotating Unbalance Response[END_REF].

The state of the friction force may change again between stick and slip. Therefore, the "inner" iteration loop is repeated until there is no change of the state of the friction forces. The convergence is possible because the normal contact forces (i.e. the Lagrange multipliers) are kept constant during the "inner" iteration loop. The advantage of the augmented Lagrange method is now clear.

The total potential energy and the "inner" iterations loop

Following the described algorithm, the total potential energy of the bump foil structure is:

Π `\] 1 2 K M P ; M K M 1 2 K < P ; < K < z {_ ̅ 1 2 Z `\] } ~•€••• ' 1 2 z Z h t , , v ~ƒ•"•… ' z 4_ ̅ ! $2 e , 3 † ‡ ‡ ‡ˆ ‡ ‡ ‡‰ Š ‹," t , , v ~ƒOE"• ' 1 2 z Z , ~Ž•••• ' (22)
where 'X ' are all nodes where the contacts between bumps/top foil and bumps/sleeve are assumed to be closed ( 'X ' 5 2 ), ' 'u and 'd are the nodes with friction state stick and slip, respectively 4 and "c c are the nodes where the rotor interferes with the top foil "c c 5

. The third term on the right hand side of eq. ( 22) corresponds to the potential energy due to bump normal contacts (thus the gap functions are either , or , ), the fourth term is the potential energy due to nodes that are in stick state, the fifth term corresponds to the mechanical work of slip friction forces and the sixth term is the potential energy due to the interference between 4 $6 6 !6 @ !B p , !6 @ ∩ !B p Ø the rotor and the top foil. The penalty parameters acting like springs Z `\] , Z h , Z have different values and generally, Z `\] ≪ Z h , Z .

For example, the equations of the algorithm can be detailed for a foil structure with one bump following the notations used in Figure 3. The extension to an arbitrary number of bumps is immediate.

Π `\] 1 2 K M P ; M K M 1 2 - _ ̅ 2 , 3 1 2 Z `\] 2 , 3 _ ̅ 2 ,E 3 1 2 Z `\] 2 ,E 3 1 2 Z h t , , v ˜™š˜ -4_ ̅ ! $2 e , 3 † ‡ ‡ ‡ˆ ‡ ‡ ‡‰ Š ‹,• t , , v ˜™š˜ -1 1 2 Z h t ,' ,' v ˜™š˜ - 4_ ̅ ! $2 e ,' 3 † ‡ ‡ ‡ˆ ‡ ‡ ‡‰ Š ‹,› t ,' ,' v ˜™š˜ - 1 1 2 Z ( 23
)
where ISTICK=1 when the contact is stick and ISTICK=0 otherwise (i.e. slip).

The minimization of the total potential energy following the degrees of freedom of the problem yields:

oeΠ `\] oe , z -,% ,% E %' Z h t , , v ˜™š˜ -4_ ̅ ! $2 e , 3 ˜™š˜ -1 0 oeΠ `\] oe , z -,% ,% _ ̅ Z `\] 2 , 3 E %' 0 oeΠ `\] oe ,' z -',% ,% E %' Z h t ,' ,' v ˜™š˜ - 4_ ̅ ! $2 e ,' 3 ˜™š˜ - 1 0 (24) oeΠ `\] oe ,E z -E,% ,% _ ̅ Z `\] 2 ,E 3 E %' 0 oeΠ `\] oe - _ ̅ Z `\] 2 , 3

Z 0

In a matrix form:

• ž Ÿ ; M 0 0 - ¡ Q R R R S Z h ˜™š˜ - 0 0 0 0 0 Z `\] 0 0 Z `\] 0 0 Z h ˜™š˜ - 0 0 0 0 0 Z `\] 0 0 Z `\] 0 0 Z `\] Z T U U U V ¢ £ ¤ † ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ˆ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡‰ ¥¦ §¨© ª . / 0 / 1 , , ,' ,E « / ¬ / - . / / 0 / / 1 Z h , ˜™š˜ -4_ ̅ ! $2 e , 3 ˜™š˜ -1 _ ̅ Z `\] Z h ,' ˜™š˜ - 4_ ̅ ! $2 e ,' 3 ˜™š˜ - 1 _ ̅ Z `\] _ ̅ Z `\] Z « / / ¬ / / - (25)
The sign function is discontinuous and thus strongly non-linear. Therefore, it is replaced by a continuous approximation and a Newton-Raphson algorithm explained in Appendix 2 can now solve the non-linear system.

The "outer" iterations loop

After solving for friction forces, the Lagrange multipliers are updated following eq. ( 15) in an "outer" iterations loop. It was previously stated that the Tresca approximation of the friction force described by eq. ( 17) can be interpreted as a solution search inside of a cylinder. Figure 6 then shows how Tresca friction cylinders approximate the Coulomb friction cone. A structured pseudo-code of the "inner" and "outer" iterations loops for friction and normal forces, respectively, is given in Figure 7.

The "inner" and "outer" lops can be also identified on the flowchart depicted in Figure 17 of the Appendix 3. supposed open and the value of the gap function becomes negative. This is a third iterative loop that encompasses the "outer" one. Figure 8 displays the pseudo-code of the contact algorithm with friction. The algorithm contains three embedded loops. The first loop detects contacts, the second loop calculates the normal forces and the third loop the friction forces. The calculation of the nodal displacements given by the non-linear system (25) could be considered the fourth innermost loop.

The flow chart of the complete contact algorithm is depicted in Figure 17 of the Appendix 3. The decision algorithm describing the tests and the actions taken for updating the close/loose state of the contacts can be identified in the lower part of this flowchart.

Figure 8. Algorithm for contact states detection and updating

Numerical results

The contact model is applied for the static analysis of the bump foil bearing structure described in [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF]. It is a first generation foil bearing that was thoroughly analyzed in the literature and could now serve as a workbench. The characteristics of the foil structure are depicted in Figure 9 and listed in Table 1. The numerical tests consist of progressively pushing the rotor against the foil structure in positive and negative X and Y directions (Figure 1) by incrementing the eccentricity up to maximum radial displacement and then pulling the rotor back to the centered position. The force response of the compliant structure is calculated by summing the projected contributions of loaded bumps.

Foil structure free of manufacturing errors

The first comparisons are made with the model presented in [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] and for the bearing free of manufacturing errors. The results are depicted in Figure 10 for cases with friction (f=0.1) and without.

Figure 10 depicts only the direct response forces (parallel to the rotor displacement), cross coupling forces being of small magnitude. The maximum radial displacement of the rotor was 79.5 µm (i.e.

2.5Cr

). This means the rotor interfered with the top foil over 47.7 µm. The former contact model from [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] considered Coulomb friction but was based on an algorithm that made no distinction between the static loading and unloading. Therefore, only a single (dotted) curve can be obtained. The present model makes a clear distinction between the loading and unloading part of the cycle with friction and therefore predicts the expected hysteresis loop. When friction is absent, both model predict a single loading-unloading curve.

The influence of friction on the static load is clearly shown by the results depicted in Figure 10. A friction coefficient of 0.1 leads to an increase of 20-60% of the static load depending on the loading direction.

Figure 10 shows that for the +X and -Y directions the results obtained with the contact model from [START_REF] Lez | A New Bump-Type Foil Bearing Structure Analytical Model[END_REF] and the present model coincide while differences are present for the -X and +Y directions. These two latter loading cases correspond to radial displacements in the direction of the welding point of the foils (Figure 1). The differences can be explained when analyzing the load supported by each bump. Figure 11 depicts the load on each bump for the maximum radial displacement of the rotor (i.e. 2.5Cr or 79.5µm). It can be seen that for the -X direction, the first bump close to the welding is much more loaded than the rest of the bumps. While respecting this trend, the new model predicts a lower force for the first bump while the next bumps (from 2 to 5) are slightly more loaded. This result holds for both cases, with and without friction. This difference is due to the top foil that was absent in the former model and that enables a more uniform distribution of the bump load. The load predicted by the new model for the -X direction is therefore slightly lower.

For the +Y direction, the results in Figure 10 obtained with friction show a slightly larger force predicted by the new model while when friction is absent both models give the same results. Figure 11 shows that the forces predicted by the new model with friction are larger for bumps 2 to 4

(especially for bump 2) while being identical with the previous model for the rest of the bumps. This difference can be due to a combined effect of the top foil and the friction force. In all cases, the presence of the top foil in the new model extends the coupling between bumps while decreasing slightly the load gradient from bump to bump. were presented in [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF]. In that reference, both the top and the bump foils were discretized with conventional shell elements and the sleeve was considered rigid. The interactions between the bumps, the sleeve and the top foil were taken into count with a penalty contact algorithm. Friction was considered by using the law of Coulomb. The results of the two models are very close but not identical. The differences come from the fact that in the non-linear structural model, the top foil and the bump foils are modeled as shells while the rotor and the sleeve were considered rigid. Surface contacts may then occur between discretization nodes and elements of the top and bump foils and are dealt with in a master-slave approach [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. During loading and unloading, the master node may successively lay on the faces of different slave elements. Moreover, the contact is not reduced to a single pair of master-slave elements: many neighboring master nodes may lay on neighboring slave elements. This is a consequence of modeling the top foil and the bump foil as shells.

The model presented in this work for the bump foil is much simple. The bumps are replaced by trusses and the contacts between the elements of the foil structure may occur only in a node-to node approach: every node of a foil knows a priori the node of the other foil that it may contact.

Even if a more elaborate model would have been used for the top foil the results would not have been different because the node to node contact is mainly imposed by the simplified bump foil model.

The present approach has the advantage of being very efficient in terms computation time compared to the complete structural model. For example, a loading-unloading cycle necessitates hours in the complete non-linear structural model and only minutes in the simplified model with 100 equally spaced intermediate steps.

Foil structure with manufacturing errors

Further comparisons between the new model and the full non-linear elasticity model from [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF] can be made in the case of the foil bearing with manufacturing errors. Bump foil height manufacturing errors can be imposed in the new simplified model by using the gap functions.

It was considered that the manufacturing errors followed a normal distribution with an average value equal to the theoretical bump height (i.e. 508 µm) and a given standard deviation.

The results obtained for a rotor displacement of 2.5Cr toward +X are depicted in Figure 13 and Figure 14 for 10µm and 20µm standard deviations (i.e. 2% and 4% of the design bump height), respectively.

The heights of the 26 bumps were randomly selected from a normal distribution. They are depicted as bars for each calculation case and are compared with the design bump height and the position of the rotor surface.

Six different cases are depicted in Figure 13 and Figure 14: the five cases presented in Figs. 12-13 of reference [START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF] plus an additional one. The prediction of the new simplified model coincide well with full non-linear elasticity model. For some cases, the radial forces predicted by the new simplified model are higher but the hysteresis loop delimited by the loading-unloading cycle is always well predicted.

The bump heights depicted as bars in Figure 13 and Figure 14 give also an insight of how the radial clearance is affected by manufacturing errors. The numerical analysis always needs a radial clearance for completely describing the geometry of the foil bearing. This design clearance can be even negative when it corresponds to an interference between the rotor and the top foil. However, due to the elasticity of the foil structure, the radial clearance is quite difficult to measure accurately even for an ideal foil structure, free of manufacturing errors. Therefore, the experimental measurements give a radial clearance affected by quite high uncertainties. For example, if the results in Figure 13 and Figure 14 would have been experimental data, the clearance in the +X direction could be estimated as the radial displacement of the rotor for which the direct structural force response Fx is lower than a given threshold. Figure 15 depicts the +X clearance estimated from Figure 13 and Figure 14 for 5N and 10N thresholds for Fx. The incertitude appears quite clear when comparing these results with the 31.8µm design value of the radial clearance (dotted line).

The difficulties encountered when measuring the radial clearance of real foil bearings (i.e. affected by manufacturing errors) can be explained by analyzing the deformation of the bumps and of the top foil during the loading and unloading process. The right hand columns depict the same results in a cylindrical coordinate system. The successive snapshots in Figure 16a-e correspond to the points indicated in Figure 14. Figure 16a shows that interferences between the rotor, the top foil and the bump foil may be present even if the rotor is centered, although the bearing was designed with a 31.8µm radial clearance. This situation can occur even for lower standard deviation of the manufacturing errors.

Figure 2 .Figure 3 .

 23 Figure 2. Structural model of the bump foil used in [32] (3DOF per bump)

Figure 5 .

 5 Figure 5. The friction cone of the Coulomb law

Figure 6 .Figure 7 .

 67 Figure 6. The approximation of the Coulomb friction cone by Tresca friction cylinders

Figure 9 .

 9 Figure 9. The geometry of the bump foil structure

Figure 10 .

 10 Figure 10. Direct structural response forces for rotor push and pull; comparisons between the new model and the model from [32] (left f=0.1, right f=0)

Figure 11 .

 11 Figure 11. Bump loading for 2.5Cr radial displacement of the rotor; comparisons between the new model and the model from [32] (left f=0.1, right f=0)

Figure 12 .

 12 Figure 12. Direct structural response forces for rotor push and pull; comparisons between the new model and the non-linear results from[START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF] (f=0.[START_REF] Branagan | Compliant Gas Foil Bearings and Analysis Tools[END_REF]

Figure 16 depicts

 16 Figure 16 depicts the deformed top and bump foils for successive +X loading displacements. The left hand side and the mid column figures depict the unwrapped top and bump foils following the bearing circumference. For a better visualization of the contacts, the left hand column shows only the upper part of the bumps depicted as bars while the mid column depicts the entire bump height.

  Figure 16b-e show how the top and bump foil are deformed by the displacement of the rotor. It can be observed that not only contacts between the bumps and the top foil change state (close/loose) but also contacts between the top foil and the rotor and between the bump foil and the casing may change state.

Figure 13 .Figure 14 .Figure 15 . 13 5 5 16 _, 1 Π 3 Superscripts²Figure 17 .

 131415135161317 Figure 13. Direct structural response forces for rotor push and pull in "+X" direction; comparisons between the new model and the non-linear results from[START_REF] Fatu | Numerical Analysis of the Impact of Manufacturing Errors on the Structural Stiffness of Foil Bearings[END_REF], bump heights manufacturing errors σhb=10µm

Table 1 Geometric characteristics of the foil structure

 1

	Axial bearing length (mm), L	38.1
	Shaft radius (mm), R	19.05
	Radial clearance (μm),	31.8
	Number of bumps, Nbump	26
	Foil thickness (mm), ) h Young modulus (GPa), E	0.102 214
	Poisson's ratio, υ Bump pitch (mm), p	0.29 4.572
	Bump length (mm), 2B	3.556
	Bump height (mm), hb	0.508

The angle θ is measured from the welding point as shown in Figure1.

The Tresca approximation of the friction force described by eq. (18) can be interpreted as a solution search inside of a cylinder instead of a cone.

For the present problem, these are lateral displacements of the bump top and bottom points. For a general case, these are relative displacements between the two sliding surfaces.

Acknowledgement

The authors are grateful to Centre National d'Etudes Spatiales (CNES) and to Airbus Safran Launchers for supporting this work. εX=250% Bump number Bump number Radial distance [m]

Figure 16. Top and bump foil deformation for successive +X rotor displacements (20 µm standard deviation of the bump height error). (left and center: unwrapped, right: Cartesian coordinate system )

Summary and conclusions

The present work introduces a new simplified model for taking into account the interactions with friction and gaps in a bump-type foil bearing structure. Simplified models are used for the bump and the top foil in order to reduce the computational effort. However, the contacts are dealt with in a very rigorous manner, by taking Coulomb friction into account and the possibility of having close/loose contacts. Three gap functions were employed for taking into account the close/loose states: between the bump foil and the sleeve, between the bump foil and the top foil and between the top foil and the rotor.

The algorithm was tested by simulating a static loading of the foil structure. The rotor was progressively pushed against the foil structure up to a maximum radial displacement and then pulled back to the centered position. The results of the present model were compared with simulations performed with a previous model that could not consider loose contacts and with a full non-linear elasticity code where the foils were modeled as shells. The results compared well and dissimilarities were explained by the different contact characteristics : the present simplified model deals with point contacts while the non-linear model handles contacts between surfaces.

The present model is able to deal very easily with manufacturing errors. This was proved by considering normally distributed errors of the bump height. It was shown that although the radial clearance of aerodynamic foil journal bearings is an important design parameter, its definition in a bearing with manufacturing errors is very difficult and prone to large measurement incertitude.

The contact algorithm proved to be very robust and efficient even when dealing with irregular gaps and friction forces.

Appendix 1

The stiffness matrix of the 4DOF/bump model depicted in Figure 4 can be cast as the sum of individual stiffness matrixes of each bump plus a coupling matrix.

Equation (A1) suggests how the stiffness matrix of the corrugated foil with two bumps can be extended to an arbitrary number of bumps.

Appendix 2

For a robust solution of the non-linear system of equation ( 25), the discontinuous sign function must be approximated by a regular function. An approximation with the atan function is presently employed although other functions can be used (a.e. tanh)

The non-linear system for the displacements now yields: 30

O P .

Following the Newton-Raphson algorithm, the iterative solution writes:

where the Jacobian matrix is: 6

The inverse of this symmetric and positive definite matrix is quite rapid.

The time scale τ h was calibrated with a trial and error approach. It should be underlined that the present numerical approximation of the Coulomb friction depends on two user calibrated parameters: Z h and τ h . They can be interpreted as the stiffness of the contacts between surface roughness and the time scale needed for changing the sign of the sliding velocity.

Appendix 3

A general flowchart of the complete contact algorithm is given in Figure 17. 6 The DIAG function indicates a matrix with the first diagonal given by the argument vector.


P.-Y Mechin

V Keryvin

J.-C Grandidier

D Glehen

An experimental protocol to measure the parameters affecting the compressive strength of CFRP with a fibre micro-buckling failure criterion

Keywords: carbon fibres, thermosetting resin, strength, mechanical testing

Understanding axial compressive failure mechanisms and estimating the related strength in continuous fibre composites is of paramount importance in the design of their parts. The mechanism at stake is the micro-buckling instability of fibres, which is contained by the matrix. In experimental measurements, compressive strength in bending is consistently higher than in axial compression. Indeed, the induced strain gradient provides an additional containment, namely structural effect. The characterisation of all geometrical and materials properties is usually lacking for models describing these combined mechanisms. A comprehensive experimental protocol is proposed in this paper to measure all the input parameters involved in a design oriented failure criterion. An epoxy matrix/high-modulus carbon fibre composite material illustrates this protocol. The influence of some key parameters, including the initial misalignment of the fibre, is discussed, thanks to additional experimental results.

Introduction

Carbon Fibre Reinforced Polymers (CFRP) are widely used in many engineering fields (aeronautics, offshore, renewable energies, marine industry. . . ) and specific custom applications (racing cars, racing yachts. . . ). They bring both high specific stiffness and strength. From an engineering point of view, composite parts undergo compressive stresses (compression or bending loadings). In the case of wings (mainsail wing or the wing of an aircraft), bending loads will induce tension on one side and compression on the other side.

Masts of racing yachts undergo quasi pure compression loadings due to the base tuning (dock-tuning) to ensure a sufficient tension in cables. Based on the observation that compressive strengths are usually lower than tension ones [START_REF] R R Effendi | Failure mechanism analysis under compression loading of unidirectional carbon/epoxy composites using micromechanical modelling[END_REF], compressive failure should be considered carefully. Therefore, a proper estimation of compressive strength will allow for an efficient design of structures.

In terms of predictive modelling, several approaches have been suggested [START_REF] Rosen | Mechanics of composite strengthening[END_REF][START_REF] Argon | Treatise on Materials Science and Technology[END_REF]10,[START_REF] Budiansky | Compressive failure of fibre composites[END_REF][START_REF] Berbinau | Effect of off-axis ply orientation on 0 degrees fibre microbuckling[END_REF][START_REF] Pimenta | A micromechanical model for kinkband formation: Part II -Analytical modelling[END_REF][START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF][START_REF] Soutis | Compressive strength of composite laminates with an open hole: Effect of ply blocking[END_REF]. The widely accepted scenario leading to failure consists in the displacement of the fibre constrained by the matrix stiffness (see Figure 1). Unidirectional plies (UDs) of a laminate exhibit an initial misalignment of fibres (φ 0 ) due to the manufacturing process. When the compressive stress on the UD (σ in Figure 1) increases, the misalignment of the fibre increases. (Actually the term fibre waviness is more appropriate but thie former term will be used throughout the paper.) The matrix contains the propensity of the fibre to buckle under this compressive loading, mainly via its elastic stiffness. When the stress increases, the matrix yields and exhibits plasticity, resulting in a lower stiffness than the elastic one, and buckling is made easier. The matrix undergoes a huge shearing which reaches its maximum at the location of the maximum misalignment along the fibre. Bending of the fibre occurs, generating a deformation gradient through the fibre's section. Different scenarios of failure can then occur: namely fibre breakage; matrix failure or interface debonding. It usually occurs in the fibre (locally in bending) for CFRP with a specific failure mechanism named kinking, detailed by Argon [START_REF] Argon | Treatise on Materials Science and Technology[END_REF] or Budiansky and Fleck [START_REF] Budiansky | Compressive failure of fibre composites[END_REF].

Early models used a specific relationship between the macroscopic compressive stress on the UD and local shearing of the matrix due to the misalignment of the fibre [START_REF] Argon | Treatise on Materials Science and Technology[END_REF]10,[START_REF] Budiansky | Compressive failure of fibre composites[END_REF]. The difference between all these models lies in the way the shear behaviour of UD or its constituents is described: from pure elastic to elasto-plastic with or without hardening. Among recent contributions to this field [START_REF] Laurin | Laminated composite structures subjected to compressive loading : A material and structural buckling analysis[END_REF], Pimenta et al. [START_REF] Pimenta | A micromechanical model for kinkband formation: Part II -Analytical modelling[END_REF] proposed an analytical micro-mechanical model taking into account the misalignment of the fibre and an elasto-perfecly plastic matrix. Extension of this work was developed by Gutkin et al. [START_REF] Gutkin | Micro-mechanical modelling of shear-driven fibre compressive failure and of fibre kinking for failure envelope generation in CFRP laminates[END_REF] by establishing a failure envelop under both compression and in-plane shear. Since fibre buckling may initiate from the specimen free edge, where the fibre support is substantially reduced, Berbinau [START_REF] Berbinau | A study of compression loading of composites laminates[END_REF] considered the importance of estimating the through-thickness stresses and include their influence on the stability of the UDs, especially the inter-laminar shear stresses of neighbouring plies. He nevertheless concluded [START_REF] Berbinau | Effect of off-axis ply orientation on 0 degrees fibre microbuckling[END_REF] that the effect of resin softening and fibre waviness were the most significant parameters.

In terms of measurement of compressive strength, direct methods of com-pressive testing are usually difficult to carry out. Applying a pure compressive stress field without any stress raisers and preventing the buckling of the specimens is tricky and compression tests show a lot of variability [START_REF] Michael | The effect of specimen size on the bending strength of unidirectional carbon fibre-epoxy[END_REF][START_REF] Hart | Garteur compression behaviour of advanced CRFP[END_REF][START_REF] O Anthoine | Pure Compression Testing of Advanced Fibre Composites[END_REF]. As quoted in [20, chap. 5] by F.J. Matthews, It is clear that whichever compression method is used, great care must be taken with specimen preparation, operator training and in the execution of the test. Indirect and easier to carry out methods, such as bending, are alternatively employed [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF]. They are also close to service use loadings in composites structures. However, significant differences were observed between the direct and indirect ways of testing [START_REF] Michael R Wisnom | On the high compressive strains achieved in bending tests on unidirectional carbon-fibre / epoxy[END_REF].

Indeed higher compressive stresses may develop in bending as compared to the compressive strength measured in carefully carried out compression experiments [START_REF] Michael R Wisnom | On the high compressive strains achieved in bending tests on unidirectional carbon-fibre / epoxy[END_REF]. A structural effect accounting for the strain gradient induced by bending was then introduced by Drapier et al. [START_REF] Drapier | Structure effect and microbuckling[END_REF] and Gardin et al. [START_REF] Gardin | Homogenized models for the modelling of instability in long fibre media[END_REF].

Further evidence of this strain gradient effect was also observed by Wisnom and Atkinson in constrained buckling tests [START_REF] Wisnom | Constrained buckling tests show increasing compressive strain to failure with increasing strain gradient[END_REF]. This structural effect takes into account both the effect of the deformation gradient and the influence of off-axis neighbouring plies in laminates. It indicates that the compressive failure depends on the laminate sequence (thickness of the successive UDs [START_REF] Allix | Modelling and identification of the mechanical behaviour of composite laminates in compression[END_REF], stiffness of neighbouring plies [START_REF] Berbinau | Effect of off-axis ply orientation on 0 degrees fibre microbuckling[END_REF]) and the deformation gradient.

A model accounting for the combined contributions of fibre micro-buckling [START_REF] Budiansky | Compressive failure of fibre composites[END_REF] and the structural effect (stacking sequence, strain gradient) [START_REF] Drapier | Structure effect and microbuckling[END_REF][START_REF] Wisnom | Constrained buckling tests show increasing compressive strain to failure with increasing strain gradient[END_REF][START_REF] Gardin | Homogenized models for the modelling of instability in long fibre media[END_REF] was proposed to describe the mechanism of failure in compression under flexural loading [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF]. This model is also analytical so that it can be used for the fast design loop of composite structures. It requires a number of parameters to properly estimate the compressive strength. With some few exceptions [START_REF] Jumahat | Fracture mechanisms and failure analysis of carbon fibre / toughened epoxy composites subjected to compressive loading[END_REF][START_REF] Jelf | The failure of composite tubes due to combined compression and torsion[END_REF], the necessary parameters for the reported compressive strength es-timations are not measured or assessed in a comprehensive way. (Neither of them considers the structural effect.) This is the case for the latter model [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF]. The purpose of this study is therefore to establish a clear and comprehensive experimental protocol to measure all the parameters involved in this latter model [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF], being geometrical or material, required to use the design oriented failure criterion to estimate the compressive strength of UD.

The paper is organised as follows. We will first recall the salient features of this model. We will then establish an experimental protocol, able to measure all the necessary parameters of this latter model. A composite material composed of an epoxy matrix with high modulus carbon fibres will be used as an example. The measured parameters will then be used to estimate the compressive strength of such a CFRP. The experimental value of compressive strength on this material will be measured in parallel by performing bending tests. Eventually, the estimated and experimental values will be compared and the sensitivity of some key parameters on the estimation will be discussed.

A design oriented criterion for estimating the compressive strength

of unidirectional plies by fibre micro-buckling accounting for a structural effect

Fibre micro-buckling

Budiansky and Fleck [START_REF] Budiansky | Compressive failure of fibre composites[END_REF] modelled the non-linear shear behaviour of UD γ (shear strain) vs. τ (shearstress) by using a Ramberg-Osgood (RO) description:

γ γ y UD = τ τ y UD + 3 7 × τ τ y UD n UD (1)
It involves three material parameters (γ y UD , G UD , n UD ) or (γ y UD , τ y UD , n UD ). G UD is the shear modulus of the composite, τ y UD is a nominal shear yield stress for which the secant modulus is reduced to 70% of its initial value G UD . γ y UD is defined as the elastic strain at the shear stress τ y UD so that τ y UD = G UD ×γ y UD . The coefficient 3 7 is specific to a choice made in [START_REF] Budiansky | Compressive failure of fibre composites[END_REF] related to the definition of τ y UD .

The compressive strength, referred to as σ stab UD , and calculated using Budiansky and Fleck's model is defined as the maximum stress applied on UD before the instability of the fibre appears:

σ stab UD = G UD 1 + n UD 3 7 1 n UD φ 0 /γ y UD n UD -1 n UD -1 n UD (2)

Structural effect

The structural effect results in the increase of UD failure compressive strain and stress with respect to σ stab UD . It is due to several mechanisms that do not take place at the UD mesoscopic scale but rather at a larger, macroscopic scale (laminate). Three factors are included in the structural effect, namely the deformation gradient through the thickness due to bending loadings, the thickness of UD and the stiffness of off-axis neighbouring plies.

Experimental observations on the structural effect have been experimentally detailed and compared by Grandidier et al. [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF], ending with the proposal of an analytical model (Eq. ( 3)) taking into account both the microbuckling mechanism [START_REF] Budiansky | Compressive failure of fibre composites[END_REF] and the structural effect [START_REF] Gardin | Homogenized models for the modelling of instability in long fibre media[END_REF]. The estimated compressive strength is called critical, referred to as σ crit UD , is higher than σ stab UD and defined by:

σ crit UD = σ stab UD + 2 r g f π e b E m E f 1 -ν 2 m V f (1 -V f ) (3)
where r g f is the gyration radius of fibre (r g f = I S , with I the second moment of area of the fibre, and S the area of the fibre's cross section), E m and ν m are respectively the Young's modulus and the Poisson's ratio of the matrix, E f the longitudinal elastic modulus of the fibre, V f the volume fraction of fibres and e b the characteristic thickness of UD involved in the instability mode (a fraction of the total thickness of UD), which is discussed below.

Let us start with the definition of a critical thickness, e c . As an example, for a UD layer containing the neutral axis (see Figure 2 -a), only the part in compression is considered for critical thickness (e c ), while for a laminate with different ply orientations [0 2 , ±45, 0 2 , ±45, 0 4 ] S (see Figure 2 -b), e c is equivalent to the thickness of the two outer UD plies subjected to pure compression (the deformation gradient is neglected for small thicknesses). For a monolithic UD stacking (only 0°plies), the instability mode is located near the free edge due to the deformation gradient. As it has been demonstrated by Drapier et al. [START_REF] Drapier | A structural approach of plastic microbuckling in long fibre composites: comparison with theoretical and experimental results[END_REF], the modal shape of the plastic instability has an impact on critical values. In the case of UDs under bending only a reduced part of thickness is considered to contribute to compression failure due to the distribution of the strain. In order to easily approximate this effect, Grandidier et al. [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF] proposed the characteristic thickness to be e b = 0.4×e c which permits correlation with experimental results in similar configurations (see Ref. [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF] for details). On the other hand, in case of a laminate where UD plies located near the edge are relatively thin, these latter plies bear a high compression stress with a reduced gradient due to their low thickness. Consequently the modal shape of instability extends on all UDs thickness, and the characteristic thickness should be approximate by e b = e c as proposed in [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF].

Necessary parameters to identify

In order to use the model ( 3), a number of experiments is required. Most of them were briefly described and some parameters measured in [START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF][START_REF] Drapier | A structural approach of plastic microbuckling in long fibre composites: comparison with theoretical and experimental results[END_REF][START_REF] Jelf | The failure of composite tubes due to combined compression and torsion[END_REF],

but never in a comprehensive way. A comprehensive experimental protocol including specific tests is required and described schematically in Figure 3.

These tests involve:

• microstructural observations to measure the geometrical parameters:

optical microscopy (OM) to measure the ply thicknesses e ply as well as e c and e b and also the fibre misalignment φ 0 ;

scanning electron microscopy (SEM) to measure the fibre volume fraction V f and fibre diameter D f .

• mechanical tests on the constituents (fibre and matrix) as well as on the composite (UD) to measure the material parameters:

elasto-plastic shear behaviour of the UD ply, described by a RO behaviour (γ y UD , G UD , n UD ); -elastic properties of the matrix (E m , ν m ); -longitudinal elastic modulus of fibre E f .

A comprehensive experimental protocol to measure all the model parameters

This section has two objectives. Firstly, to introduce the composite material we used as a example of our protocol. Secondly, to present the set of experiments carried out (see Figure 3) and the results obtained.

Materials

The present study is focused on high modulus carbon fibres (HR40 -Mitsubishi) for both 0°plies (UD) and off-axis plies, which are commonly used in racing yachts design.

The associated matrix is a low curing temperature (∼ 110°C) epoxy resin system (for confidentiality reasons, the reference of the matrix is not provided). The pre-impregnated CFRP supplier was able to produce a block of matrix (100×100×20 mm 3 ). To ensure that the reticulation of the thermoset polymer is exactly the same at all locations in the block, quality control experiments were performed using Differential Scanning Calorimetry (DSC822, Mettler Toledo) to define the quality of reticulation. No differences between samples from the skin or the core of the matrix block were found. The question of a possible difference betwen in situ and ex situ properties of the matrix has been addressed by using a micromechanical numerical model (periodic boundary conditions) with the properties of both matrix and fibres (see § 3.3.2) and their volume fractions to simulate the shear response of a UD.

Comparisons with experimental results (see § 3.3.3) were found very conclusive. They are not reported here for sake of clarity and space limitations.

As for the laminate, the following stacking sequence was used: [+45, 0 9 , -45, 0 11 , -45, 0 9 , +45]. To ensure the highest quality of the composite, compaction is carried out every three plies to reduce void content in composites. Due to this compaction, plies located on the mould side are more compacted than those located on the vacuum bag side. Therefore, the properties of composites on the mould side and the vacuum bag side are different.

Coupons have been heated up in an autoclave according to standard curing cycle recommanded by manufacturer and pressure is set to 3 bars throughout the curing cycle.

Microstructure

The analysis of the microstructure allows one to measure the parameters related to the structural effect: the thickness of the UD plies, the volume fraction of fibres and the radius of these fibres. For the analysis of the fibre micro-buckling mechanism, the measurement of initial existing misalignment in the composite is also required.

Laminate or macroscopic scale

The side of the composite laminate is polished by means of standard grinding techniques. The thickness of each ply e ply of the composite laminate is then measured by optical microscopy (OM, Leica MZ16) and digital image analysis (ImageJ software [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF]). Close to the mould side, e ply = 292 ± 2 µm, and e ply = 300 ± 3 µm close to the vacuum bag side.

Sub-ply or microscopic analysis

SEM (Jeol JSM 6460 LV) allows measurement of both the volume fraction of fibres (V f ) and their associated diameters (D f ). The latter is needed to calculate the gyration radius (r g f ). By using a magnification of x500 and an energy of ∼ 20 kV, SEM precise views of fibres, as seen in Figure 4, are possible.

The measurement of V f and D f is carried out by digital image analysis (ImageJ - [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF]). However, other solutions exist to measure the volume fraction of fibres, such as the calcination method [START_REF] Yee | A TGA technique for determining graphite fiber con-tent in epoxy composites[END_REF]. This latter method consists in heating the composite in order to burn-off the resin without affecting the fibres (glass fibres or carbon fibres). The difference in mass between the initial composite and the fibres alone at the end of the measurement makes it possible to identify the mass fraction of matrix then back calculate the volume fraction of fibres. This measurement gives an average value of V f obtained on a large volume while SEM measurements are localised measurements on reduced surfaces raising the issue of the representativeness of the measurement (Figure 4). The results were equivalent between measurement by imaging and calcination method (for the volume fraction of fibres considered), giving V f = 54 ± 2 %.

Initial misalignment of fibres

Measuring the initial misalignment of fibres φ 0 is one of the key-points for estimating σ stab UD . Wisnom [START_REF] Wisnom | The effect of fibre misalignment on the compressive strength of unidirectional carbon fibre/epoxy[END_REF] estimated a loss of 50% in strength with a φ 0 of 3°with respect to a perfect alignment. Based on this literature, the most accurate value of the initial misalignment is required. The measurement of φ 0 is however difficult to make. Indeed, it is a measurement over a microscopic size width but over a length of several millimetres, so the length-to-width ratio is extremely large.

The distribution of fibre misalignment angle (φ 0 ) was measured using Yugartis' method [START_REF] Yurgartis | Measurement of small angle fiber misalignments in continuous fiber composites[END_REF]. The method consists in cutting samples through the fibres with a specific angle. The cross sections of fibres are therefore elliptic.

By properly polishing the cutting face, the dimensions (height and width) of the ellipses were measured (see Figure 5). The initial misalignment of the fibre was calculated as proposed in [START_REF] Yurgartis | Measurement of small angle fiber misalignments in continuous fiber composites[END_REF]. This method was used in the study of Lee and Soutis [START_REF] Lee | Thickness effect on the compressive strength of T800 / 924C carbon fibre -epoxy laminates[END_REF]. About 1000 or more ellipse lengths were measured for each face (mould or vacuum bagging). The angle distribution was then calculated and sorted into the class interval width at the angle of 0.5°. The distributions were found to be centered at an angle close to the cutting angle.

They were then transformed and centered to zero as shown in Figure 6. The standard deviation for the mould and vacuum bagging faces were found to be 0.61°and 0.83°, respectively, in agreement with values reported in literature [START_REF] Lee | Thickness effect on the compressive strength of T800 / 924C carbon fibre -epoxy laminates[END_REF].

Mechanical tests

The model (Eq.( 3)) requires measurement of the axial elastic modulus of the fibre E f and both the Young's modulus E m and the Poisson's ratio of the matrix, as well as the non-linear shear behaviour of the ply (UD).

Tensile test on matrix

Tensile tests were performed on the sole matrix, following ISO standard [START_REF]Plastics -Determination of tensile properties -Part 2 Test conditions for moulding and extrusion plastics[END_REF], to measure its elastic properties with a universal testing machine (Instron 5567, 10 kN load cell). The longitudinal and transverse strains were measured using a non-contacting 3D optical system (Aramis 5M, GOM, Germany). Two video cameras (8-bit gray levels, 2050x2448 pixels images at 2 Hz) recorded the motion of a random pattern spray painted on the sample surface, by means of a digital image correlation (DIC) technique. Five dog-bone samples (gauge length of 25x5x2 mm 3 ) were tested for reproducibility at 1 mm/min. Load-unload cycles were performed to ensure that the response of the material was elastic.

Tensile test on fibres

Tensile tests were performed on individual fibres following ASTM standard [START_REF] Astm D | Standard Test Method for Tensile Strength anf Young's Modulus for High Modulus Single-Filament Materials[END_REF]. Fibres were extracted from uncooked pre-impregnated carbon plies and isolated from matrix with solvent (acetone). They were then attached on a paper frame to obtain a gauge length of 10 mm. The average apparent diameter was determined as previously detailed using SEM. The frame was clamped on a universal testing machine (MTS, 2 N load cell), and loaded at a constant crosshead-displacement rate of 1 mm/min up to fracture. Tests were performed on fifty fibres for reproducibility. The post-processing of the longitudinal modulus was performed taking into account the compliance of the loading frame.

Experimental results on the axial modulus of fibre (E f = 361 ± 51 GPa) are in good agreement with the datasheet's values (∼ 390 GPa).

Tensile test on a ± 45°laminate: shear behaviour of UD

A central point in the use of the model (Eq.( 3)) is a relevant identification of the non-linear behaviour of the UD under shear loading. Pure shear loading on a UD is not easy to perform as suggested in [START_REF] Astm D ; M | Standard Test Method for Shear Properties of Composite Materials by the V-Notched Beam Method[END_REF]. In such a test, the UD is subjected to local effects that modify the pure shear behaviour. An alternative way is to carry out tensile tests on a +/-45°stacking according to ISO standard [21] (in-plane shear tests). Such tests were performed using a universal testing machine (Instron 5567, 30 kN load cell). In such a laminate, a ply at +45°or at -45°undergoes a quasi-pure shear loading equivalent to what happens at the location of maximum misalignment under compression. 10 mm bi-axial strain gauges (Kyowa) were glued on coupons at 0°and 90°w ith respect to the tensile axis. Ten samples with the dimensions 250x25x2 mm 3 composed of 8 plies were tested for sake of reproducibility.

The shear stress is defined as τ UD = F 2bh with F the applied load, b the width and h the thickness of the coupon. The shear strain is defined as γ UD = ε 0 -ε 90 , with ε 0 the strain given by the 0°strain gauge and ε 90 the strain given by the 90°strain gauge. The resulting τ UD -γ UD graphs are plotted in Figure 7. They are highly reproducible and only a typical curve is represented.

A non-linear regression to extract the RO parameters (Eq. ( 1)) is made for each sample, with a maximum shear strain γ fit = 2% which is applied to all samples. This latter value is chosen since compressive failure occurs at lower values of shear strain. This fit is superimposed on the same Fig. 7.

The resulting RO parameters (γ y UD , G UD , n UD ) are shown in Table 1. Let us highlight that this RO fit is by no means a constitutive equation able to deal with plasticity or possible damage [START_REF] Feld | Micromechanical prediction of UD laminates behavior under combined compression up to failure: Influence of matrix degradation[END_REF] but a convenient and simple way to describe the loading curve for subsequent use in Eq. ( 2).

Of all the information that can be extracted from this tensile test, it is also advantageous to measure the longitudinal modulus (E BB ) of the laminate (± 45°) for the subsequent calculation of the experimental compressive strength of UDs (see § 4.1). This is made by using straightforwardly the axial stress information as well as the 0°strain gauge, in the linear range of the response. E BB was found to be 13.3 ± 0.4 GPa.

All the parameters measured are gathered in Tables 1 and2.

Comparison between the engineering design tool evaluations with experimental results

Experimental determination of the compressive strength by a four-point bending test

Four-point bending tests were carried out according to ASTM standard [START_REF]Standard Test Method for Flexural Properties of Unreinforced and Reinforced Plastics and Electrical Insulating Materials by Four-Point Bending[END_REF] with a universal testing machine (Instron 5567, 30 kN load cell) with a distance of 90 mm between the upper rollers, where pure bending occurs.

Polyethylene plates were put under the rollers to minimise stress concentrations. As detailed in Section 3.1, the following stacking sequence was used:

[+45, 0 9 , -45, 0 11 , -45, 0 9 , +45]. In any case discussions on the choice of the layup including ply orientation, ply thickness are out of the scope of the paper. A specific span (460 mm) was used for these high modulus fibres.

Accordingly, the specimens had the following dimensions: 500x30x10 mm 3 .

Mono-axial strain gauges (10 mm in length -Kyowa) were glued on both the compressive side (strain ε j c ) and the tension side (strain ε j t ). Ten monolithic samples were tested. Five were broken with the mould side in compression, five with the vacuum bag side in compression. Failure occurred between the upper rollers. The force on the assembly and the two strain gauges signals were recorded during loading and synchronised. The force-displacement curves were all linear up to failure.

The Euler-Bernouilli kinematics assumption is considered valid, since the thickness of the specimen is small vis-à-vis its length. Therefore the deformation gradient through the thickness is constant. The position of the neutral plane (where the longitudinal strain is null) vis-à-vis the mid-plane (λ) is

given by the two strain gauges at the top and the bottom of the specimens, and computed as [START_REF] Astm D | Standard Test Method for Tensile Strength anf Young's Modulus for High Modulus Single-Filament Materials[END_REF].

λ = -( h 2 + e g ) ε j c + ε j t ε j c -ε j t (4)
where h and e g are the thickness of the sample and the thickness of the glue between the sample and the strain gauge (taken ~0.08 mm [START_REF] Allix | Modelling and identification of the mechanical behaviour of composite laminates in compression[END_REF]). λ differs slightly from zero. It is always negative and is ∼ 0.1-0.2 mm. As quoted by Allix et al. [START_REF] Allix | Modelling and identification of the mechanical behaviour of composite laminates in compression[END_REF], the origin is to be found in the asymmetry between the longitudinal tensile and compressive elastic behaviours. We follow the identification procedure of the compressive strength proposed by these authors by assuming that the tensile elastic modulus (E T UD ) differs from the compressive one (E C UD ) and that the latter parameter decreases linearly with the compressive strain (initial modulus E C0 UD and slope α). We use their identification procedure of these three parameters by minimising simultaneously (for all measured points during loading) the gap between an estimated normal force and zero and an estimated bending moment and the measured one. We control the quality of the identification by calculating error indicators on the normal and the bending moment and making them very similar (∼ 1‰). In doing this, we use the tensile moduli of the off-axis plies and the position of each block of plies of the same nature vis-à-vis the neutral plane λ. Values of E T UD , E C0 UD and α were found to be 207 ± 3 GPa, 194 ± 3 and 13 ± 2 GPa/%, respectively.

The position of UD where the instability is assumed to take place is taken as the boundary between the topmost UD and the upper off-axis ply. Next, the strain at failure at this position is calculated as the deformation gradient (given by the two strain gauges and the thickness of the specimen) times the position of the UD. Finally, the compressive strength σ fail UD is calculated as the actual compression modulus times the strain at failure in UD. The mechanism responsible for failure was checked to be fibre microbuckling by observing the broken surfaces by OM (fibre kinking) and SEM (tensile and compressive fracture surfaces of broken fibres).

Estimation of compressive strength with the design oriented failure criterion

Using all the previously measured parameters in § 3, the compressive strength σ crit UD was estimated combining the fibre micro-buckling contribution (Table 1) as well as the structural effect (Table 2). For the former, estimations were made assuming an initial misalignment φ 0 to be either one or two standard deviations of the fibre misalignment normal distribution.

Discussion

A comparison between the measured and estimated compressive strength of the UD is made in Table 3 and Figure 8. Let us recall that it is in no case a validation of the failure criterion relevance but rather a comparison on its use in an engineering context. Additional microstructural features such a void content, resin rich regions is beyond the scope of this paper.

The contribution of the micro-buckling effect (see Table 1) is the largest part of the compressive strength prediction. The structural part contributes to ∼ 15 % of the compressive strength (with respect to the experimental values) for the stackings considered. As seen in Eq. ( 3), parameters r g f , E m , ν m , E f and V f play in turn a minor role on the structural part. Therefore, information from suppliers datasheets is sufficient in terms of fibre diameter and volume fraction or materials elastic parameters. For instance the elastic asymetry in compression w.r.t. to tension for the fibres or the matrix is not of high importance. The choice of the characteristic thickness e b , as discussed in Section 2, has, at the contrary, an influence on the structural part.

From the experimental point of view, compressive failure on the mould side provides higher values than on the vacuum bag side (∼ 15 %). According to the model, a different initial misalignment that increases the shear loading of the UD explains this difference. The analytical model looks relevant for quantifying the effect of initial misalignment on compressive strength. According to our measurements, there is only a 0.2°(s.d.) difference between the mould side and the vacuum bag side. This causes an equivalent gap of 20 % on the compressive strength. The model is highly sensitive to initial misalignment, which in turn, is actually highly sensitive to the manufacturing process. These observations are in good agreement with estimations from Ref. [START_REF] Nl Hancox | The compression strength of unidirectional carbon fibre reinforced plastic[END_REF][START_REF] Wisnom | The effect of fibre misalignment on the compressive strength of unidirectional carbon fibre/epoxy[END_REF][START_REF] Soutis | Static Compression Failure of Carbon Fibre T800/924C Composite Plate with a Single Hole[END_REF].

The micro-buckling part of the model is the one that is so dependent on φ 0 . It decreases non linearly with φ 0 . It reduces by ∼ 500 MPa from 1 to 2°( see Fig. 8). With the comparison between experiments and the estimated strength being so dependent on the choice of φ 0 , it is necessary to discuss the choice of a meaningful value or of a range of meaningful values. In most models, the initial misalignment is unfortunately taken as an adjustable or tested parameter to correlate compressive strength estimations with experiments [START_REF] Jelf | The failure of composite tubes due to combined compression and torsion[END_REF][START_REF] Wisnom | The effect of fibre misalignment on the compressive strength of unidirectional carbon fibre/epoxy[END_REF][START_REF] Nl Hancox | The compression strength of unidirectional carbon fibre reinforced plastic[END_REF][START_REF] Grandidier | A fibre direction compressive failure criterion for long fibre laminates at ply scale, including stacking sequence and laminate thickness effects[END_REF]. In several observations, the large numerical sensitivity to the initial misalignment has been experimentally compared on pultruded composites with small initial misalignment angle (misalignment is measured using X-ray scanning tomography) [START_REF] Creighton | The compressive strength of highly-aligned carbonfibre/epoxy composites produced by pultrusion[END_REF]. Creighton et al. [START_REF] Creighton | The compressive strength of highly-aligned carbonfibre/epoxy composites produced by pultrusion[END_REF] observed an increase in 25 % of the compressive strength by reducing the initial misalign-ment from 1.0°to 0.7°. However, samples with a larger φ 0 also contained a larger void content volume fraction that will necessarily also affect the compressive strength of UDs. Therefore, it is somehow difficult to consider this result as a straight validation of the effect of misalignment on the compressive strength prediction. This would be indeed very relevant to mastering a manufacturing process able to control the value of φ 0 without inducing any kinds of other defects to test whether the experimental influence of φ 0 on the compressive strength would follow that of the model. The micro-buckling effect will affect predominantly the fibres that are the most misaligned. The fibre misalignment distribution follows a normal distribution (see Fig. 6). The value usually taken for φ 0 [START_REF] Lee | Thickness effect on the compressive strength of T800 / 924C carbon fibre -epoxy laminates[END_REF][START_REF] Jumahat | Fracture mechanisms and failure analysis of carbon fibre / toughened epoxy composites subjected to compressive loading[END_REF] is one standard deviation of the normal distribution. It represents a coverage probability of 68 %, meaning that 32 % of the fibre population will contribute to the instability. Such a choice of φ 0 is plotted in Fig. 8. For both mould and vacuum bagging faces, the model largely over-predicts the experimental results, as already related in the literature [START_REF] Jumahat | Fracture mechanisms and failure analysis of carbon fibre / toughened epoxy composites subjected to compressive loading[END_REF]. Another choice is to consider that the instability will be triggered by a smaller population of more misaligned fibres. Taking a value of φ 0 corresponding to 2 s.d. of the distribution gives a coverage probability of ∼ 95 %. Such a choice is also plotted in Fig. 8. In that case, the comparison is much better between estimations and measurements. It means that ~5 % of the most misaligned fibres will trigger the micro-buckling instability.

The robustness of the model was tested for some other parameters that were thought of being of some influence on the compressive strength estimation. The compressive strength (micro-buckling of the fibre, σ stab UD ) depends linearly on the shear modulus of the UD (G UD ). A variation of 10 % of this modulus increases or decreases σ stab UD by ∼ 9 %. Its measurement is made without too much uncertainty regarding the behaviour and is reproducible over a large number of samples. On the other hand, an appropriate value for the characteristic thickness (e b ) can be discussed depending on the deformation gradient effect. From an analytical point of view, a variation of 10 % in e b induces an increase or decrease of σ crit UD by ∼ 1 %. Therefore, it is not necessary to pay further attention to the measurement of ply thickness, e ply , at least for monolithic composite laminates [START_REF] Mechin | An Analytical Approach to Estimate the Compressive Strength of Carbon Fibre Reinforced Plastics[END_REF].

Conclusion

The axial compressive strength of unidirectional plies is one of the key design parameters for composites in structural engineering. Several mod-

  els have been developed to describe the specific mechanism of compressive failure. The one used in this study is a failure criterion based on the microbuckling of the fibre, controlled by the shear elasto-plastic behaviour of the matrix. An additional contribution concerns a structural effect linked to the beneficial effect of a strain gradient, for instance induced by a flexural loading. This criterion is design oriented so that it can be used for the fast design loop of composite structures. This paper presented a comprehensive experimental protocol to measure all the necessary input parameters of the criterion, which was lacking. It is by no means a validation of this failure criterion. It allows for the measurement of geometrical parameters such as fibre volume fraction, fibre radius, thickness of the plies and initial fibre misalignment. It also permits experimentally determination of material parameters such as the non linear shear elasto-plastic behaviour of the ply, the elastic properties of the matrix and the elastic longitudinal modulus of the fibre. An epoxy matrix/high modu-lus carbon fibre composite was used as a representative example of this very protocol.The compressive strength was then estimated by means of the measurement of these parameters. A comparison between this estimation and experimental values obtained by four-points bending tests with the same composite material was then made. The estimation of the compressive strength was found to be highly sensitive to the initial misalignment of fibres in the ply, in relative accordance with experimental observations. The choice of a proper value of this initial fibre misalignment required by the model was discussed in view of the fibre misalignment distribution and the experimental values. In no case, this paramater was chosen as a fiting variable. Estimations and experiments are closely related when 2 standard deviations of the fibre misalignment distribution is made. It also suggests that 5 % of the most misaligned fibres trigger the instability.

-Figure 1 .

 1 Figure 1. Schematics of compression failure mechanisms leading to kinking band formation. The blue points correspond to the position of maximum axial stress through fibre's section (compression or tension) while the red surface is the area of the matrix exhibiting plasticity and extending due to the micro-buckling under compression.

Figure 2 .

 2 Figure 2. Schematics for defining the critical thickness (e c ) and the characteristic thickness (e b ) involved in the structural effect (see Eq. (3)) and especially the case for the deformation gradient due to bending. Two examples are shown: (a) [0] 20 laminate, (b) [0 2 ,±45, 0 2 , ±45, 0 4 ] S laminate.

  (3)) and especially the case for the deformation gradient due to bending. Two examples are shown: (a) [0] 20 laminate, (b) [0 2 ,±45, 0 2 , ±45, 0 4 ] S laminate.

Figure 3 .Figure 4 .

 34 Figure 3. Schematics of the experimental protocol proposed to measure all the parameters for the compressive strength model of Eq. (3)

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. OM picture showing the ellipses, formed by the fibres after cutting, used to obtain the fibre misalignement distribution using Yurgartis method

Table 1 .

 1 Experimental parameters for the micro-buckling contribution σ stab UD , see Eq. (2) for details. The non linear shear behaviour of UD is described by a RO fit (see[START_REF] Allix | Modelling and identification of the mechanical behaviour of composite laminates in compression[END_REF]) with parameters G UD , γ UD , obtained from a tensile test on a ±45°laminate. φ 0 is the initial misalignment on both sides of the composite, taken as the standard deviation (s.d.) of the fiber misalignment distribution using Yurgartis method.

	Attention must therefore be paid on convenient and accurate ways to mea-
	sure the initial misalignment and its influence on compressive strength, with-
	out introducing other manufacturing defects such as porosities, that could
	affect the shear behaviour of the ply.

y

UD and n

Table 2 .

 2 Experimental parameters for the structural effect σ stru UD , see Eq. (3) for details. The fibre volume fraction (V f ) and diameter (D f ) are obtained by SEM and digital image analysis. The elastic properties (Young's modulus E m , Poisson's ratio ν m ) of pure matrix, is obtained by tensile testing. The longitudinal elastic modulus of fibre (E f ), is obtained from fibre tensile tests.

		σ fail UD	σ crit UD (1 s.d.) σ crit UD (2 s.d.) Unit
	Mould side	1221 ± 101 1762 ± 71	1264 ± 68 MPa
	Vacuum bag side 987 ± 115	1532 ± 63	1073 ± 54 MPa

Table 3 .

 3 Experimental compressive strength of UD, σ fail UD (obtained by means of a fourpoints bending test) compared to the critical compressive stress estimated σ crit UD for one or two standard deviations (s.d.) of the fibre misalignement distribution.

retrieve/pii/0266353887900169.

Acknowledgements

The authors are grateful to C. Baley (IRDL -UBS) for encouraging discussions, to F. Portanguen, K. Henry, A. Jegat, A. Magueresse, M. Grabow (UBS) for experimental support. We are also grateful to CDK Technologies (Racing yacht builder, France), especially S. Digard who provided all composites coupons and supported this program. Undergraduate students A. Boyer and O. Chabaat (UBS) are also acknowledged for experimental help. Partial fund-

Data availibility

The raw/processed data required to reproduce these findings cannot be shared at this time due to legal or ethical reasons.