N Leboeuf

T Boileau

B Nahid-Mobarakeh

N Takorabet

F Meibody-Tabar

G Clerc

Modeling of PM Synchronous Machines Under Inter-turn Fault

This paper presents a theoretical approach as far as faulty Permanent Magnet Synchronous Machines (PMSM) are concerned. The knowledge of PMSM's behavior under inter-turn fault is very useful in order to check drive performances under fault conditions. Control, stability and fault detection often require faulty models. The aim of this paper is to develop a dynamical model of a PMSM under fault conditions without using circuit models coupled with finite element analysis (FEA). The method exposed in this paper only needs some preliminary FE computations that can be used in order to take into account of different winding distributions.

I. INTRODUCTION

One way to study PMSM under inter-turn short-circuits is to build a improved model. Many authors have developed PMSM models including electrical faults [START_REF] Dehkordi | Real Time Simulation of Internal Faults in Synchronous Machines[END_REF], [START_REF] Bui | Noninvasive Fault Monitoring of Electrical Machines by Solving the Steady State Magnetic Inverse Problem[END_REF] or mechanical faults [START_REF] Ebrahimi | Static Eccentricity Fault Diagnosis in Permanent Magnet Synchronous Motor Using Time Stepping Finite Element Method[END_REF], [START_REF] Tabatabaei | Modeling and Simulation of a Salient-Pole Synchronous Generator With Dynamic Eccentricity Using Modified Winding Function Theory[END_REF]. Such models often require strong hypotheses. Core saturation, leakage inductances are often neglected and winding arrangements are simplified [START_REF] Reichmeider | Internal Fault in Synchronous Machines: Part I: The Machine Model[END_REF]- [START_REF] Shah | Analysis of turn to turn Fault in Surface PM Machines with Multi-Layer Fractional-Slot Concentrated Windings[END_REF]. This leads to a lack of accuracy on results. The results can be strongly improved by using time-stepping finite element method (TSFEM). However, this approach needs to perform a new simulation for every fault configuration and requires too much CPU time. In addition, TSFEM cannot be easily coupled with a faulty controlled model implemented under MATLAB/SIMULINK ® environment. It is also not very appropriate for parametric studies [START_REF] Ebrahimi | Static Eccentricity Fault Diagnosis in Permanent Magnet Synchronous Motor Using Time Stepping Finite Element Method[END_REF], [START_REF] Vaseghi | Inductance Identification and Study of PM Motor with Winding Turn Short Circuit Fault[END_REF], [START_REF] Birolleau | Etude d'actionneur pour la tolérance aux défauts[END_REF].

The presented work concerns non-salient PMSM. The authors present a weak coupled model using FEA to determine elementary inductances of the different coils. An appropriate method gives the right self and mutual inductances of the equivalent external circuit of the faulty machine. Then, thanks to a precise identification of the various parameters (inductances, no-load back-EMF), the study of the behavior of the faulty machine can be easily performed whatever the fault configuration is.

II. MODEL OF PM MOTOR UNDER FAULT CONDITIONS

We assume that saturation of iron materials is negligible. In addition, only the non-salient PM machines are considered. Under these hypotheses, the following model can be adopted.

Consider a 3-phases PM motor with inter-turn short circuit. The main idea consists in considering an additional circuit in the model as shown on Fig. 12 Each phase is composed of series connected coils which need to be properly modeled [START_REF] Vaseghi | Inductance Identification and Study of PM Motor with Winding Turn Short Circuit Fault[END_REF]. The electrical equation of such model can be written of the form: where, 9 4567 and C 4567 can be expressed as follows: 9 and C 6 8 C ' are respectively the phase resistance and the self inductance of the healthy machine. 9 4 andC 4 are the resistance and the self inductance of the faulty winding + and F 7 is its no-load EMF. ' 4(4 ,' +, and ' 46 are respectively the mutual inductances between the winding + and the windings+ 1 , , and -. A 7 and 9 7 are the fault current and the fault resistance (Fig. 1.). . / is the electrical position of the rotor. It has been shown that the parameters of the different parts affected by the shortcircuit are directly linked to the ratio 0 8 1 7 1 2 between the number of shorted turns and the number of total turns of one phase [START_REF] Tallam | Transient Model for Induction Machines With Stator Winding Turn Faults[END_REF], [START_REF] Farooq | Modeling and simulation of stator inter-turns faults in Permanent Magnet Synchronous Motors[END_REF]. This assumption is exact in the case of windings with one slot per pole and per phase and supposing one pole pair machines. In fact, the previous equations %3 * depend on the winding arrangement. In the case of non conventional windings or windings with more than one slot per pole and per phase, a more rigorous approach need to be adopted. The aim of this paper is to develop a general method based on FEA and winding arrangement modeling which leads to determine the parameters of the faulty model for various fault configurations.

III. INDUCTANCE MATRIX CALCULATION

The main idea presented concerns the use of a minimum number of FE simulations in order to obtain the inductance matrix of the equivalent scheme of a PM motor.

A. Healthy conditions

A 4-poles, 1 / -slots PM non-salient machine is considered with 1 6 turns per coil in the case of a concentric winding arrangement. Fig. . shows the example of a 15slots 16poles PMSM. Field calculation is performed with linear FEM. If only one coil 7is supplied by the current A 8 , the fluxes 9 : in the different coils can be determined by:

9 : 8 ; < = > 2 ?@ A2 DB C @ A2 DB D E A 8 13 2 2 2 1 / F
where C G is the axial length of the machine, A is the potential vector, H 6 is the cross section area of the coil. H I and H J denote the go and return areas of the coil. The computation the various inductances are performed under the assumptions presented in section II. The fluxes of the different coils obtained with the supply currents A 8 allow us to determine the inductance C 8 of the coil 7 and the different mutual inductances ' 83: . Self and mutual inductances can be divided in two parts: C 8 3 ' 83: and C K 3 ' K matching respectively with common flux 9 6 and leakage flux 9 L as shown on Fig. 2. Even if 9 L does not cross the air-gap, it contributes to the mutual inductance between the two coils. This is specific to concentric windings. Hence, assuming that the leakage flux does not depend on current level, the energy M in a whole slot with currents A ( and A is given by:

M 8 1 2 2 C K 2 A ( B 1 2 2 C K 2 A B ' K 2 A ( 2 A 6
Consequently, C K and ' K can be easily determined using M3 A ( and A . In the case of concentric windings and considering the symmetry of the machine, the previous calculation is performed once and permutation technique leads to obtain the inductance matrix of the healthy system of 15 independent coils:

B. Windings arrangements

The approach presented above can be performed to obtain the inductance matrix of both single or double layer distributed windings for healthy machines. Here, we present the case of a double layer winding assuming that the same approach can be done for single layer arrangement. We consider a *5-slots 5-poles PMSM with a double layer winding including shortened coils with a coil pitch FW6 as shown on Fig. 6. The whole inductance matrix C is a symmetric and circulant matrix. As each slot contains the wires of two different coils, each coil has two leakage mutual inductances with two other coils sharing the same slot. On Fig. *, we can see that coil 1 has two leakage mutual terms with coils 6 and **. This can be clearly seen on Fig. 5 which shows the representation of the first row of the inductance matrix 52 The location of the mutual leakage terms ' K on a row is depending on the coil pitch. The proposed approach can also be applied for other more complicated winding arrangements such as fractional arrangements. In all cases, the inductance matrix of the 3-phase winding can be obtained according to the winding arrangement:

C 456 8 X Q3Y Z 2 C2 X Q3Y Z [
where X Q3Y Z is the winding matrix.

C. Parameters identification under fault conditions

We consider the case of the 15-slots 16poles PMSM. We assume that 1 7 turns of coil 1 are short-circuited. The ratio 0 becomes 0 \ 8 1 7 1 [START_REF] Vaseghi | Inductance Identification and Study of PM Motor with Winding Turn Short Circuit Fault[END_REF] 2 . Using local considerations, we will modify the expression of the healthy inductance matrix [START_REF] Tallam | Transient Model for Induction Machines With Stator Winding Turn Faults[END_REF] in order to obtain a formulation corresponding to the faulty case. Inside the slot which contains the coil 1 and 2 for instance, the short-circuit leads to consider the configuration shown on Fig. 8 where A ( ]] 8 A ( A 7 and A ( ] 8 A ( according to Fig. 1. We can determine the various leakage inductances by using energy formulation. The whole energy in the slot is given by:

M 8 1 2 2 C K( ] 2 A ( ] B 1 2 2 C K 2 A B 1 2 2 C K( ]] 2 A ( ]] B ' K( ] 3 2 A ( ] 2 A B ' K3( ]] 2 A 2 A ( ]] B ' K( ]] 3( ] 2 A ( ] 2 A ( ]] !1
Using this equation, we can deduce the value of the leakage and mutual inductance terms for a given faulty case. The self and mutual inductances of the additional circuit can be given by: In addition to the self and mutual inductances corresponding to the common flux 9 6 , the self and mutual inductances corresponding to the leakage flux 9 L have to be properly identified using 1.

The leakage inductancesC K( ] and C K( ]] are added in the diagonal terms of the new matrix whereas the mutual leakage inductances ' K(3( ]] ,' ^13 and ' K( ]] 3 are added in the non-diagonal terms. The final matrix (13) for the faulty machine is defined as:

C 8 1 0 \ 2 C ( B C K( ] 1 0 \ 2 ' (3(N B ' K( ] 3 1 0 \ 2 ' (3(O P P 1 0 \ 2 ' (3Q 1 0 \ 2 ' (3 B ' K( ] 3 0 \ 2 1 0 \ 2 C ( B ' K(3( ]] 1 0 \ 2 ' (3 B ' K( ] 3 C ( ' (3(N P P ' (3R ' (3Q 0 \ 2 ' (3 B ' K( ]] 3 T T T T 1 0 \ 2 ' (3(N B ' K( ] 3 ' (3(O ' (3(U P P ' (3 B ' K C ( B C K 0 \ 2 ' (3(N B ' K( ]] 3 0 \ 2 1 0 \ 2 C ( B ' K(3( ]] 0 \ 2 ' (3(N B ' K( ]] 3 0 \ 2 ' (3(O P P 0 \ 2 ' (3Q 0 \ 2 ' (3 B ' K( ]] 3 0 \ 2 C ( B C K( ]] " # # # # # # # $
By this way, an inter-turn fault leads to modify the healthy matrix V of order 1 / into a faulty matrix 1% of order 1 / B 1. The same method can be applied on the double layer presented winding by modifying 5.

1 Finally, using the a modified winding matrix X R3(_ , we can easily identify the various terms of the inductance matrix C 4567 of the faulty equivalent model with the following relation:

C 4567 8 X R3(_ 2 C2 X R3(_ 1*
1 Other parameters such as 9 4567 and F 7 also need to be identified in order to get a complete model. 9 4567 can be determined assuming that resistive effects are directly linked to the number of turns. The no-load back-EMF F 7 is also supposed to be proportional to the number of shorted turns. For the considered 1 7 turns in fault, we have: where a bcd is the magnet flux through each turns of the considered coil, . / is the electrical rotor position with respect to the faulty coil, 4 the number of poles pairs,` the mechanical speed. Spatial harmonics can also be included.

IV. COMPARISONS AND DISCUSSIONS

The parameters identification presented above can be validated with a full finite element computation. It will be also compared to another approach presented in [START_REF] Tallam | Transient Model for Induction Machines With Stator Winding Turn Faults[END_REF], [START_REF] Farooq | Modeling and simulation of stator inter-turns faults in Permanent Magnet Synchronous Motors[END_REF]. Then, the whole developed model has to be tested on the 15-slots 16-poles PMSM.

A. Parameters Comparisons

Fig. 7 shows the comparison between the results obtained by the proposed approach and a full FEA. They are also compared with another approach which does not consider the exact winding topology and using the following parameters [START_REF] Tallam | Transient Model for Induction Machines With Stator Winding Turn Faults[END_REF], [START_REF] Farooq | Modeling and simulation of stator inter-turns faults in Permanent Magnet Synchronous Motors[END_REF]: According to the definition given above, 0 \ and 0 are related by: 0 8 0 \ 2 1 1 6 2 1V

where 1 1 6 2 is the number of coils per phase which is equal to 6 in the considered machine. The various computations are performed for different values of the ratio 0 \ on coil 1. A good agreement is observed between the proposed approach and the FEA method whereas formulas in 16 seem to give approximate results. As mentioned in section II, these results concerns some restricted cases. Nevertheless, in order to be more precise on armature reaction and saturation effects, the real dynamic inductances must be considered.

B. Dynamic simulations

The model presented in section II and the parameters identified in III allow to perform the dynamic simulation of faulty operations. We consider that half of the turns of coil 1 are short-circuited (0 \ 8 2F) by resistance 9 7 g 1`. The PMSM is supplied by a Voltage Source Inverter (VSI) and controlled using a speed scheme at `8 194h. Fig. 6 shows the current waveforms using parameters obtained the proposed method (set 1) and those calculated by 16(set ). The phase currents change with the parameter set. The fault current calculated with set is higher than the other. Due to the lower value of the inductances obtained by 16, fault current also contains non negligible high frequency components in this case. Set lead to higher phase currents in order to maintain the speed reference with the considered fault. In fact, results given by set can only be correct where 9 4 f C 4 2 i. Consequently, we can state that the approximate parameters give wrong estimations of the variables in dynamic conditions. Moreover, conclusions that can be taken with set are biased and can lead to "missed fault" or "false alarms" issues. V. CONCLUSION

The method proposed in this paper provides a new way to simulate PMSM under inter-turn faults. Parameters are properly identified for various fault configurations and the approach can cover any type of windings arrangements. This prevent to consider wrong estimations of the variables regarding dynamic results. Compared to TSFEM, it does not require a large CPU time and the accuracy of the results is not affected. Moreover, thanks to the theoretical expressions, parametric studies can be easily performed.

Fig. 1 .

 1 Fig. 1. The external model 3 4567 8 9 4567 2 A 4567 B C 4567 2 DA 4567 DE B F 4567 1 A 4567 8 A 4 A 5 A 6 A 7 3 4567 8 3 4 3 5 3 6 3 7 F 4567 8 F 4 F 5 F 6 F 7

Fig. 2 .Fig. 3 .

 23 Fig. 2. Considered PM motor and the winding arrangement for a half machine

Fig. 4 .Fig. 5 .

 45 Fig. 4. Double layer winding arrangement for one pole

1

  Additional row and column are inserted in the whole inductance matrix corresponding to the faulty coil. The self and mutual inductances C ( 3 ' (3: A 8 3 2 2 2 1 / are affected by 0 \ :

Fig. 6 .

 6 Fig. 6. Faulty case considered

F

Fig. 7 .Fig. 8 .

 78 Fig. 7. Fault inductances: comparison between FEA, proposed method and other approaches


D Trimarchi

M Vidrascu

D Taunton

S R Turnock

D Chapelle

Wrinkle development analysis in thin sail-like structures using MITC shell finite elements

Keywords: Sail Modeling, MITC Shells, Wrinkling, Wrinkling Index

We propose a method of modelling sail type structures which captures the wrinkling behaviour of such structures. The method is validated through experimental and analytical test cases, particularly in terms of wrinkling prediction. An enhanced wrinkling index is proposed as a valuable measure characterizing the global wrinkling development on the deformed structure. The method is based on a pseudo-dynamic finite element procedure involving non-linear MITC shell elements. The major advantage compared to membrane models generally used for this type of analysis is that no ad hoc wrinkling model is required to control the stability of the structure. We demonstrate our approach to analyse the behaviour of various structures with spherical and cylindrical shapes, characteristic of downwind sails over a rather wide range of shape and constitutive parameters. In all cases convergence is reached and the overall flying shape is most adequately represented, which shows that our approach is a most valuable alternative to standard techniques to provide deeper insight into the physical behaviour. Limitations appear only in some very special instances in which local wrinklingrelated instabilities are extremely high and would require specific additional treatments, out of the scope of the present study.

Introduction

Sail-like thin laminates, and more generally inflatable structures, have been analysed with Finite Element (FE) methods, see e.g. [START_REF] Durand | Unsteady numerical simulations of downwind sails[END_REF][START_REF] Graf | RANSE investigations of downwind sails and integration into sailing yacht design processes[END_REF]. In such studies the fabric has generally been modelled using the membrane model, which assumes that the bending components are negligible.

The constant strain triangle (CST) membrane element has often been adopted for inflatable structures [START_REF] Oden | Finite strains and displacements of elastic membranes by the finite element method[END_REF][START_REF] Tabarrok | Nonlinear analysis of tension structures[END_REF]2] or sail-type analysis [START_REF] Graf | RANSE investigations of downwind sails and integration into sailing yacht design processes[END_REF][START_REF] Roux | Strongly coupled VPP and RANSE code for sailing yacht performance prediction[END_REF][START_REF] Malpede | A fully integrated method for optimizing fiber-membrane sails[END_REF], for which the ratio t/L -where t is the fabric thickness and L a characteristic size -is very small. The success and longevity of this model are partly due to its simplicity; the CST membrane is not without problems, however, and examples can be found where this class of model produces incorrect answers [START_REF] Prukl | The CST -An insidious survivor from the infancy of FEA[END_REF].

More complete iso-parametric models have been adopted for general inflatable structure analysis [START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Bletzinger | Structural optimization and form finding of light weight structures[END_REF]. There are, however, particular combinations of load and boundary conditions for which the membrane problem is ill-posed. When this happens, for example in point-supported structures such as downwind sails (see Figure 1), the membrane model itself is inadequate to represent the structural deformation, as the state of the structure is bending-dominated, or of mixed membrane-bending type. a b In such situations a series of oscillations arise on the deformed fabric surface, which is then said to be "wrinkled", as in Figure 1. The development of such oscillations is buckling related, thus mostly controlled by the bending stiffness. The membrane model being unable to represent such configurations, singularities appear and convergence cannot be reached. Generally, these difficulties are overcome using ad hoc wrinkling models [START_REF] Kang | Finite element analysis of dynamic response of wrinkling membranes[END_REF][START_REF] Lu | Finite element analysis of membrane wrinkling[END_REF], and applications examples for sail deformation analysis can be found in [START_REF] Heppel | Accuracy in sail simulation: Wrinkling and growing fast sails[END_REF][START_REF] Renzsch | Fluid structure interaction simulation of spinnakers -getting closer to reality[END_REF].

Similar behaviour can be identified also along the free edges of the sailwhere the membrane stiffness is ineffective [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF] -and such border singularities are often avoided by introducing cable or beam elements as supporting frames [START_REF] Li | An integrated analysis of membrane structures with flexible supporting frames[END_REF].

An alternative approach is proposed here using a finite element shell model. It is then possible to numerically reproduce the wrinkling development. The suitability of using shell finite elements for planar physical membranes has been shown by [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF][START_REF] Wang | Wrinkling prediction of rectangular shell-membrane under transverse in-plane displacement[END_REF]. In the present paper, however, we analyse a fully curved geometry. The intrinsic ill-posedness of the problem induced by the above-mentioned physical instabilities, deriving from the very small thickness of the fabric, is overcome with the use of Newmark-type dynamic schemes with Rayleigh damping, see e.g. [START_REF] Bathe | Finite Element Procedures[END_REF]. Our work is thus concerned with applying and assessing this type of procedure when used with effective shell elements to model sails undergoing wrinkling. Some specific geometries were selected as representative of downwind sail shapes. These structures were then loaded with constant pressure, the value of which lies in the typical range for downwind yacht sails. Although complex interactions exist between the deformed shape -including the wrinkles -and the distribution of the actual wind loading [START_REF] Trimarchi | The use of shell elements to capture sail wrinkles, and their influence on aerodynamic loads[END_REF], the work presented here focuses on the structural response and provides therefore a basis for the coupled FSI analysis [START_REF] Parolini | Mathematical models and numerical simulations for the America's Cup[END_REF][START_REF] Viola | Sail pressures from full-scale, wind-tunnel and numerical investigations[END_REF][START_REF] Trimarchi | Analysis of downwind sail structures using non-linear shells finite elements -wrinkle development and fluid interaction effects[END_REF].

The aim of this work is then to computationally characterize the nature and the behaviour of the wrinkling development with a particular concern for downwind sail configurations. Section 2 describes the method adopted for the analysis of thin wrinkled shells. In Section 3 this method is assessed based on experimental results and analytical solutions published in the literature [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF], for a flat membrane with prescribed boundary displacements. The nature of wrinkling is then discussed in terms of stability in Section 4, by an eigenvalue analysis of the tangent stiffness matrix of the structure. The next sections are devoted to specific applications of the proposed procedure to sail configurations, with detailed sensitivity studies with respect to various parameters characterizing discretization, geometry and material properties. In Section 5, mesh sensitivity is examined with particular regard to wrinkling generation. Section 6 then presents a parametric analysis for a spherical geometry to investigate the influence of fabric thickness, Young's modulus, and the applied load on the wrinkle wavelength and the number of wrinkles. Then Section 7 examines the effect of shape variations in Gennakerand Spinnaker-type geometries. Finally, we give some concluding remarks in Section 9.

Methods

In this Section we outline the major methodological concepts used in our analysis of wrinkled shells.

Shell kinematical assumptions and finite element approximation

Non-linear iso-parametric four-node shell elements of the MITC family (Mixed Interpolation Tensorial Components) have been used for the finite element analysis. The detailed formulation of MITC shell elements can be found in [START_REF] Bathe | Finite Element Procedures[END_REF] and in [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF], where detailed considerations on the shell physical behaviour are also presented. For completeness, however, we will provide a summary of the kinematical assumptions and of the finite element discretizations considered.

The underlying mathematical shell model is based on the so-called Reissner-Mindlin kinematical assumption, stating that any material line initially orthogonal to the shell midsurface in the undeformed configuration remains straight and unstretched in all deformed states. Note that we do not assume further -as in the Kirchhoff-Love kinematical assumption -that the deformed material line remains orthogonal to the deformed midsurface, see e.g. [START_REF] Batoz | Modélisation des Structures par Eléments Finis[END_REF][START_REF] Nguyen-Thanh | Rotation free isogeometric thin shell analysis using pht-splines[END_REF] and references therein for shell elements based on the Kirchhoff-Love assumption. We thus allow for all types of deformation -including transverse shear -in the material.

The specific shell element considered -namely, the MITC4 element -was selected for its remarkable effectiveness, robustness with respect to various numerical pathologies such as locking affecting shell analysis, and versatility, as demonstrated in [START_REF] Bathe | An evaluation of the MITC shell elements[END_REF][START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF]. In particular, this element has been shown to accurately represent both membrane-and bending-dominated types of behaviour, and to be also effective for extremely thin structures, hence removing the need for switching between so-called "thick and thin" shell elements depending on the application considered.

Denoting the local coordinates by (r, s, z) varying between -1 and +1 inside an element, the Euclidean position vector of a point within this element for a configuration indexed by is given by

X = 4 k=1 λ k (r, s) X (k) + z t (k) 2 V (k) n , (1)
where

X (k) , V (k) n
and t (k) respectively denote the position vector, unit transverse direction vector1 , and thickness parameter at node k in configuration , while λ k is the associated finite element shape function. Here, since we have four-node elements, standard bilinear interpolation functions are used. Then, based on the iso-parametric construction discrete displacements within a total Lagrangian formulation are obtained by subtracting the initial position, viz.

U = 4 k=1 λ k (r, s) U (k) + z t (k) 2 V (k) n -0 V (k) n , (2)
where U (k) then denotes nodal displacements with respect to the initial configuration. Note that, under a small displacement assumption -which we do not use in our modelling -the difference between transverse direction vectors could be parametrized in the standard form

V (k) n -0 V (k) n = -α (k) 0 V (k) 2 + β (k) 0 V (k) 1 , where 0 V (k) 1 , 0 V (k) 2 , 0 V (k) n
defines an orthonormal basis in the reference configuration -hence, 0 V

(k) 1 , 0 V (k) 2
are unit vectors lying in the tangential plane to the midsurface -and (α (k) , β (k) ) denote associated rotation angles. In our case, we consider large displacements and this parametrization by rotation angles is only used in the linearized increments needed in the Newton-type algorithm [START_REF] Bathe | Finite Element Procedures[END_REF].

These kinematical assumptions are then incorporated in a 3D variational formulation with a plane stress assumption to obtain the discrete (finite element) problem. We emphasize that, in MITC elements, special care is exercised to circumvent specific difficulties associated with thin structures [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF]. In particular, various forms of numerical locking must be avoided while preserving the stability in the energy. To that purpose, in MITC elements specific interpolation rules for the various components of the strains are defined based on well-chosen "tying-points". Some other shell finite element procedures instead rely on reduced integration, but specific controls are required in order to avoid the zero energy modes called hourglass modes [START_REF] Zienkiewicz | The Finite Element Method: its Basis and Fundamentals[END_REF]. In our computations we used the MITCNL software developed at Inria, an implementation of the non-linear MITC shell elements.

Pseudo-dynamical scheme

The finite element analysis of thin structures such as sail type laminates is likely to encounter convergence issues. For such materials, typical aspect ratios -namely, thickness over characteristic size -varies between 10 -4 and 10 -3 . Since wrinkling modes are essentially controlled by the bending stiffness -proportional to the cube of the aspect ratio [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF]-this leads to configurations where very little stability prevails, and even some intermediate stages -e.g. in incremental loading -where no stable static solution exists. Various techniques can be adopted in order to overcome such instabilities, and we here choose to employ a pseudo-dynamic procedure with artificial damping. Namely, we compute a static solution by following a dynamical path in which we use the natural modelling of the inertia and stiffness effects, and incorporate some Rayleigh damping defined in an ad hoc manner. Of course, the time variable in this procedure is then a "pseudo-time", which should not be considered as the physical time. The dynamic equation solved is then of the type (see e.g. [START_REF] Bathe | Finite Element Procedures[END_REF])

M Ÿ + C Ẏ + K(Y) = F(t), (3)
where Y denotes the vector of displacement degrees of freedom, F the load vector, M the usual (consistent) mass matrix, K(Y) the non-linear stiffness matrix-operator, and

C = c 1 M + c 2 K 0 , (4)
with c 1 , c 2 constant values and K 0 a reference stiffness matrix to be specified -e.g., the tangent stiffness matrix in the undeformed configuration. For a linear behaviour with K(Y) = K 0 Y we can use modal decomposition based on the eigenmodes of K 0 with respect to M, and the damping ratio for a given natural frequency ω then gives

ξ = c 1 + c 2 ω 2 2ω . (5)
For a linear harmonic oscillator, tuning this coefficient has the effect of modifying the damping characteristics of the system. No damping (ξ = 0) results in infinite oscillations, while 0 < ξ ≤ 1 leads to oscillations damped with a characteristic time constant of 1/(ξω), until no oscillation is encountered for ξ = 1 (so-called "critical damping"). For ξ > 1 the system is overdamped -i.e. takes a longer time to stabilize than for ξ = 1 -and for ξ 1 the corresponding characteristic time is 2ξ/ω. When applying Rayleigh damping, a linear structure would then behave as a set of such harmonic oscillators. Since the damping ratio is a function of the modal frequency -recall (5) -every eigenmode is then damped with a different time constant, where the fastest convergence is obtained for the critical value ξ = 1. This implies that no given eigenmode of eigenfrequency ω i can be damped any faster than with the time constant 1/ω i . Hence, a natural strategy consists in adjusting c 1 and c 2 so that the fundamental eigenmode is critically damped (ξ(ω 1 ) = 1), and in order to control the damping time constants of higher frequency modes based on the characteristic time 2ξ(ω i )/ω i . This leads to

c 1 = ω 1 , c 2 = 1/ω 1 , (6)
corresponding to a global time constant of 1/ω 1 . Of course, this simplified reasoning holds for a linear structure, and in our case it is only approximate, although it can be applied to small oscillations around a given configuration, when using the tangent stiffness operator.

Regarding the time discretization strategy, we use a mid-point rule, namely,

       Ÿn+ 1 2 = Ẏn+1 -Ẏn ∆t Ẏn+ 1 2 = Ẏn+1 + Ẏn 2 = Y n+1 -Y n ∆t = 2(Y n+ 1 2 -Y n ) ∆t (7)
which we then substitute in the dynamical equation (3) written at time "n + 1 2 ". This non-linear implicit equation is solved based on a Newton-Raphson algorithm, in which we also allow for time adaptation when convergence difficulties are encountered in the Newton loop. Typically, at the beginning of the solution process the physical instability and the rapidly-changing stiffness operator tend to produce oscillations, which can only be handled with small time-steps. As the computation proceeds the time step can be progressively increased, until it ideally takes infinite values for a steady equilibrium state. The use of a variable time-step is thus necessary both for handling stability issues and for reducing the computation time.

The resulting global solution algorithm is then summarized in Figure 2.

Measures of wrinkling

Two main measures are generally used in order to characterize surface wrinkling: the number n w -or wavelength L w -of the wrinkles appearing on the deformed surface and their amplitude A w [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF][START_REF] Wang | Wrinkling prediction of rectangular shell-membrane under transverse in-plane displacement[END_REF]. A third measure is here introduced, which describes the global state of a wrinkled section. This is defined as I w = n w •σ, where σ is the standard deviation of the section from a smooth regression curve fitting the wrinkled cross-section. This index thus takes similar values for a large number of small wrinkles or a small number of large wrinkles and of course for a smooth surface the index vanishes. Wrinkling measures were computed with an automatic wrinkling recognition algorithm for planar cross-sections. This algorithm uses a regression of the cross-section with a polynomial of order five, and a re-parametrization with the arc length of the polynomial curve. Consequently, we simply compare the transformed cross-section to the straight arc-length axis as visualized in Figure 3. Calculating the standard deviation σ, the number of wrinkles is determined by following the arc-length and counting a new wrinkle every time the curve intersects the thresholds -δ or +δ, where

δ = max(ασ, σ M IN ),
directly after a zero crossing Z c , for a well-chosen factor α. The minimum height of the threshold σ M IN is calculated as a percentage of the section length, and avoids taking into account excessively small oscillations, which e.g. may be induced by the finite element discretization. In practice, suitable values are α = 0.2, σ M IN = s M AX /1000, where s M AX is the maximum coordinate of the arc-length. Wrinkle evaluations were proposed by [START_REF] Tessler | Nonlinear shell modelling of thin membranes with emphasis on structural wrinkling[END_REF][START_REF] Wang | Wrinkling prediction of rectangular shell-membrane under transverse in-plane displacement[END_REF][START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF]. In this paper we have focused on the work of [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF], which analysed an initially planar rectangular membrane with prescribed shearing displacements on one side, as shown in Figure 4. The membrane was made of a thin film of Kapton R , the physical quantities of which are reported in Table 1. We point out that the aspect In this configuration the membrane is in a shear-dominated mode, hence it is likely to produce several wrinkles oriented at 45 • . It was reported that, although the imposed displacement was increased gradually and monotonically, changes in the number of wrinkles occur suddenly, as the effect of an instability originating near the edge. The number, amplitude and wavelength of such oscillations is dependent upon the material and the magnitude of the prescribed displacement. The measured deformed central cross-section of the membrane is plotted in Figure 5. A discrepancy should however be underlined: the declared length of the membrane is 380 mm, whereas the graph maximum x-value in Figure 5 is only ≈ 360 mm. We should also note that the symmetry properties expected from the loading considered are not entirely obtained in the response. In particular, a rather surprising irregularity in the amplitude of the oscillations is observed near the center for x > 180 mm. We conjecture that this is due to a slip in the displacement prescription, as is clearly noticeable in Figure 4 (right-hand side of top edge).

Validation test case

An initial out-of-plane displacement map -defined as a linear combination of eigenmodes -was applied on the undeformed structure for the numerical results published in [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF]. However, the influence of this map was negligible in the final result, and the need for this initialization procedure is probably due to the use of the static solver implemented in the Abaqus commercial code.

We reproduced this experiment with our proposed method. Lateral sides were left free, whereas all degrees of freedom of the bottom edge were fixed, and displacements were prescribed on the top edge. No initial out-of-plane mapping or pre-stress was needed. Due to the expected 45-degrees shear mode, we directly considered adapted meshes aligned with this privileged direction, see Figure 7a. The pseudo-dynamic algorithm was able to capture the final deformed shape. However, of course since no initial imperfections were applied, out-of-plane oscillations only appear when the top edge displacements reach a given finite value. Mesh sensitivity was explored for eight meshes, the characteristics of which are reported in Table 2. Figure 6a shows the mesh sensitivity in terms of number n w and the proposed wrinkling index I w defined in Section 2.3. Figure 6b reports the mesh sensitivity in terms of the amplitude of wrinkles A w and the standard deviation σ, defined with respect to the undeformed straight configuration of the center cross-section. A significant mesh sensitivity is experienced. This is not surprising since several elements are needed in order to reproduce every wrinkle, the wavelength of which may be very small, especially along the constrained sides. As a general rule of thumb, the mesh should exceed 20 elements per wrinkle in order to capture the correct structural behaviour. When the mesh is fine enough the results converge to the solution; this is confirmed by the good degree of accuracy achieved in the comparison with the experiments (see Section 3.2).

Mesh sensitivity for the validation test case

Comparison with experimental results

The results discussed here are those obtained with a converged mesh (630x150 elements). The general behaviour is well-captured in terms of wrinkle shape, size and orientation. Fine details are also well-represented, such as the smaller oscillations observed in between two large wrinkles near the top and bottom edges, see Figure 4. Compared to the experimental results, however, our numerical solutions feature a more regular response and the predicted number of wrinkles (20) exceeds the experimental value [START_REF] Parolini | Mathematical models and numerical simulations for the America's Cup[END_REF].

Figure 8 compares the experimental and the numerical cross-sections. Since the experimental curve is of non-zero mean, an appropriate shift -of specific value (-7, 0.8) [mm] -was applied to make the curves comparable.

This comparison shows a good accuracy in reproducing the experimental behaviour. The response is particularly accurate both in terms of amplitude and wave-length. Five zones can be identified, as illustrated in Figure 8. In Zone 1 the three "boundary" wrinkles are well-captured. A transition region (Zone 2) is then seen where the amplitude is well-captured, but five wrinkles are found instead of four. In the centre of the membrane (Zone 3), the solution is accurate in terms of wave-length. The amplitude is less accurate, but this could be due to possible imperfections in the experimental setup, as already mentioned, or in the material -constitutive and geometric -properties. Zones 4 and 5 display anti-symmetry -with respect to Zones 2 and 1, respectively -in the numerical solution, compared to which a phase inversion is seen in the experimental curve for x > 250 mm. This further justifies investigating material and boundary conditions effects, as follows.

The effect of the material was then investigated with a new analysis, where the fabric thickness was increased by 10%. The results are displayed in a b . We now have a good representation of the experimental behaviour in the two lateral parts -see new Zones 1 and 3 in this Figure -whereas the numerical wave-length of the wrinkles is still higher -although lower than for the thinner membrane -in the central Zone 2. This is also confirmed by the comparison of the Fourier transforms of the experimental and numerical solutions in Figure 8 and Figure 10, which shows a much narrower frequency content in the numerical solutions, with a peak closer to the center of the experimental frequency range in the thick case. The influence of the boundary condition was then analysed by prescribing a lateral displacement varying linearly from the value δ at the top-left corner to zero at the top-right corner, in order to model a slip in the boundary condition. The resulting wrinkles are shown in Figure 11, where we observe a less regular behaviour in the numerical solution starting around the middle of the cross-section, as for the experimental curve, although the behaviour near the right side is less accurate. Note also in Figure 11 that the Fourier transform of the numerical and experimental solutions are now very similar.

Finally, it appears that our numerical procedure adequately captures the experimental behaviour of [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF], taking into account the various experimental uncertainties.

Comparison with analytical estimates and numerical results

In [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF] some analytical derivations are performed for the case considered, under various simplifying assumptions. In particular, the wrinkles half wavelength and amplitude are respectively estimated as

λ an = πHt 3(1 -ν 2 )γ , A an = 2Ht π • (1 -ν)γ 3(1 + ν) , (8)
with H and L the height and width of the sheet, γ = δ/H the shear deformation. Substituting the above-given numerical values we find λ an = 14.4 mm and A an = 0.36 mm. In our numerical solution we have λ ≈ 16 mm and A ≈ 0.35 mm when measuring near the center, which is in good agreement with the analytical values. Numerical experiments were also performed in [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF] using Abaqus. Three regular meshes were tested, using so-called S4R5 shell four-node elements based on an approximate enforcement of the Kirchhoff-Love kinematical assumption, see Abaqus manual [START_REF]ABAQUS Theory and User's Manual[END_REF] for more detail. The first mesh element size was set slightly smaller than the wrinkle half-wavelength. The finer meshes were obtained by increasing the number of elements by factors two and four, approximately. Figure 12 visualizes the resulting deformed solution, and we also notice in this Figure a peculiarity in the prescribed boundary conditions, namely, that only translations along the x-direction were allowed on the two lateral sides. It is argued in [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF] that these additional constraints may induce two extra oscillations, which are then subtracted from the counted number of wrinkles. In order to recall the existence of these extra oscillations we explicitly list the values of [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF] as "n+2" in Table 3. We see that our converged values -recall Table 2 -are intermediate between the converged values of Table 3 with and without the additional two wrinkles.

We further note that the smaller oscillations observed in the experimental results near the top and bottom edges are much more adequately captured in our numerical solutions than in the Abaqus numerical results given in [START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF].

Wrinkles and equilibrium

As wrinkling is a form of buckling, its development can be analysed in terms of stability. An equilibrium configuration is called stable when after a small disturbance the system returns to its initial configuration. An indifferent equilibrium configuration is encountered when any perturbed configuration is an equilibrium configuration. In the unstable equilibrium configuration, there exist arbitrarily small disturbances under which the system moves to a new (not necessary stable) equilibrium configuration.

The Lagrange-Dirichlet theorem states that a minimum of the potential energy in a conservative mechanical system is a stable equilibrium configuration. As the potential energy for linear analysis is given by W = 1 2 (Y T • K 0 • Y), a correspondence can be found between the positiveness of the tangent stiffness matrix

K T (Y) = ∂K ∂Y ,
and the system stability. This is then represented by the sign of the eigenvalues, as summarized in Table 4.

To illustrate this discussion, the sail-type geometry shown in Figure 13 has been loaded with a typical constant pressure (10 -4 N/mm 2 ) and analysed in order to detect the equilibrium behaviour in relation to the generation of wrinkles. A mesh with about 8000 elements was chosen, with 86 elements per horizontal section and some adequate refinement near the corners. This was judged a sufficient mesh resolution for a limited number of eigenvalues. Of course, after discretization the number of eigenmodes is finite unlike for the continuous case, and is equal to the number of degrees of freedom N dof of the structure. The analysis specifications are summarized in Table 5. The physical quantities selected -namely the thickness t, Young's modulus E, the Poisson ratio ν and the mass per unit volume ρ s -are typical of sail-type materials. Except when otherwise clearly specified, all the numerical values used in the paper refer to those given here. Note that we will always neglect the effect of the sail self-weight, since it is roughly two orders of magnitude below that of the pressure loading. The tangent stiffness matrix was extracted every 20 time steps, and an eigen-analysis was performed for the first 10 eigenmodes of lowest associated frequencies. Figure 14-top plots the horizontal displacement of a sample node located on the structure symmetry line at height z = 750 mm, while Figure 14-bottom shows the evolution of the first eigenvalue. Instability is soon encountered after starting the calculation process, and it only disappears near the end of the time path, when the structure is established in a stable configuration. This instability is the cause of the wrinkling generation, as it appears by inspecting the eigenvectors during the time path. Figure 15 displays the first eigenvector for time steps 120, 236 and 394, all representing unstable modes. In the first case, the unstable mode spreads to the whole sail. As the calculation proceeds, the instability tends to concentrate in the wrinkled zones. At the end of the calculation no unstable mode remains, but the first eigenmodes still represent wrinkling patterns. This indicates The strong link between instability and wrinkling is confirmed when inspecting a number of eigenmodes for a single time step. Taking for instance the 394 th time step, in the middle of the unstable time period, the majority of the first eigenvectors represent the wrinkling deformation components, as shown in Figure 16 where the first six eigenvectors are reported.

Mesh sensitivity in the wrinkling representation

A parametric grid sensitivity study was undertaken in order to assess the adequate grid fineness for accurately capturing the wrinkles. The purpose was also to verify whether the converging behaviour of the finite element solution is affected by the wrinkling development.

Two geometries and 14 meshes per geometry have been adopted for this mesh sensitivity analysis. In the first geometry, a Gennaker-like structure is analysed, and the geometry is described using a spherical sector of radius R = 1000[mm], and of height and width corresponding to angles π/2 and 2π/5, respectively. For symmetry reasons, we will consider in our computations only half of this geometry along the width direction, and we display this reduced geometry in Figure 18. In fact, unless otherwise clearly specified all the geometries visualized in the rest of the paper represent half of the real geometry. Of course, in actual applications, sail shapes will very seldom enjoy such symmetry properties, and at the very least suffer from geometric imperfections, indeed. Nevertheless, we consider here these generic shapes as part of a proof-of-concept study. The second shape analysed has the same spherical part with a cylinder added underneath, in order to produce a Spinnaker-like structure as displayed in Figure 22. The height of this cylinder is equal to the radius of the spherical component. In all the meshes the element size distribution was kept as homogeneous as possible. However, the meshes were designed with finer discretization in regions close to the corners, where wrinkling generation is expected. The characteristics of these meshes are reported in Tables 6 and7. For the Gennaker type geometry, we selected two specific cross-sections for their interest regarding the wrinkling pattern, see Figure 17. The Head section is characterized by the intersection of the deformed shape with a plane, the normal vector of which is tangent to the surface and directed along the bisector line, while the height of the cross-section is defined such that approximately half of the wrinkle length is intersected. The Corner crosssection is defined by intersecting the structure with a plane passing through a point Q located on the bottom edge at the crossing with the bisector line. The plane is normal to a vector, the horizontal projection of which is tangent to the sail at the point Q, but vertically elevated at 45 • upwards.

Figure 19a shows the results in terms of displacements of the selected points shown in Figure 17, and Figure 19b shows the number of wrinkles for the Gennaker geometry. Acceptable results for the prediction of local wrin- kling are obtained starting from the mesh m 45 , which has about 20 elements per wrinkle in the Head section. This result is in line with that obtained in Section 3. Regarding the number of wrinkles captured with different meshes, it is interesting to observe the general behaviour of the curves, where discontinuities are observed when the mesh considered is able to capture new wrinkles of smaller dimensions. This is better seen in Figures 2021 deformed cross-sections are displayed for meshes of increasing accuracy. Note in particular in Figure 21b how two small additional wrinkles are captured in the finest meshes, whereas the overall behaviour is well represented for coarser meshes. points farther away from each other, and developable surfaces -referring to the cylindrical part, here -tend to undergo larger motions. Regarding the wrinkling representation, we note an adequate convergence with the sequence of meshes, except specifically in the Head cross-section. Although the curves in Figure 23 do not look well converged, some limited oscillations of the wrinkling count are due to the definition of the wrinkles recognition algorithm criterion, since in some cases the thresholds are reached or missed by a very small margin, while clear trends can be identified regarding the -hence approximate -number of wrinkles. In fact, this is the case for all sections Region near the half-side of the structure except the Head section, which displays a high variability in the number of wrinkles. This is due to the progressive resolution of small amplitude wrinkles as the mesh becomes finer, see in particular Figure 25. We further discuss this issue in Section 8.

For both the Gennaker and the Spinnaker cases, the mesh sensitivity analysis indicate that the wrinkling is a higher-order phenomenon, hence a relatively coarse mesh is able to represent the general behaviour of the sail, whereas a detailed wrinkling representation requires at least 20 elements per wrinkle.

Influence of material properties variations

A set of tests was carried out using the Gennaker geometry and nine different combinations of material properties, all of which are realistic for sailtype materials. The parametric changes concerned Young's modulus and the thickness. Starting with an initial value, both Young's modulus and thickness parameters were multiplied and divided by a factor two, thus defining the test grid reported in Table 8. Based on the results of the mesh sensitivity analysis, it was then decided to use the mesh referenced as 45 for the following evaluations, since this mesh provides good convergence properties of the solution, recall Figure 19. The same two probe Head and Corner cross-sections were considered as in Section 5. The deformed cross-sections are visualized in Figure 26, and the corresponding wrinkling indicators are plotted in Figure 27. Although no specific relationship can be identified regarding the dependency of these indicators with respect to the parameters considered, this study shows that our procedure can adequately capture the wrinkling behaviour for the whole parametric range of concern. Globally, the amplitude and the number of wrinkles are principally influenced by the thickness parameter, and as expected the thinner the material the larger the number and the amplitude of the wrinkles, for a given pressure load. By contrast, no clear general trend can be observed when varying Young's modulus.

A comparison can then be made when doubling the load. In this case the wrinkling pattern is not strongly affected, and similar features can be identified. This comparison is shown in Figure 28 for cases number 1 and 5, representative of the whole spectrum of analysed cases. We observe that the deformed sections preserve the same essential features, thus with a stronger dependency upon the material properties than on the loading magnitude.

Influence of the geometry

In this Section we investigate the influence of the geometry on the wrinkling development, starting with Gennaker-like shapes of varying widths, and proceeding with Spinnaker-like shapes of various heights. Five different spherical geometries were analysed for this investigation of the effect of the shape on wrinkling development. These geometries were derived from that already analysed in Section 5 by varying the half-width angle φ previously defined as π/5, see Figure 29. The mesh density was kept as constant as possible by preserving the ratio of number of element by width angle, which results in the number of elements per horizontal section EpS listed in Table 9. All physical quantities remained set as in Case 5 of Table 8.

Spherical shape with varying width

The resulting computed deformed Head cross-sections -as identified in Figure 17 -are shown in Figure 30. We first observe that all cross-sections have roughly similar extents, which is informative as regards the location of the wrinkles. More generally, the shapes of the deformed Head and Corner cross-sections -wrinkling included -are quite comparable for the various geometries, except for φ = π/2 in the Head section for which the global curvature is inverted, indeed. This behaviour can be explained by observing that the extension of the wrinkles is in this case very small, hence the crosssection is located very near the attachment and the force transmission area associated with the wrinkles tends to pull the sail downwards.

In the Corner section, the number of wrinkles as well as the average wrinkle amplitude are largely influenced by the base angle. As expected, the wider the angle -i.e. the larger the structure -the higher the number and amplitude of wrinkles, see also Figure 29. Concerning the Head sections, the wrinkles present a more stable behaviour, both in terms of number and of amplitude. The general wrinkling behaviour, which can be observed in the screen shot in Figure 29, is well represented by the wrinkling indicators in Figure 31. In particular the wrinkling index indicates well the large increase in the wrinkling development on the Corner section as the geometry becomes wider.

Spinnaker-like shapes with varying heights

Here we consider Spinnaker-like shapes similar to that defined in Section 5 and vary the height h of the cylindrical part with values h/R = 1; 1.5; 2, in addition to the previous case h/R = 0.5. All other quantities remain unchanged. Results are reported in terms of the deformed sections already shown in Figure 22, where the Head-base section is always defined as the junction of the spherical and cylindrical parts, and the so-called Mid section is taken at half the height of the cylindrical part.

As we increase the height of the cylindrical body, the overall pattern formed by the displacements -namely, along the force application direction in the lower part of the structure, and in the opposite direction in the higher part -is more and more marked, see Figures 32 and33. In addition, in the higher part the spherical shape progressively becomes flatter as the width decreases due to the distance to the lower boundary conditions, while sharp wrinkles appear and contribute to the overall width reduction, indeed. The wrinkling behaviour is analysed in more detail in Figure 34 which provides for each shape and section the number, the average amplitude of wrinkles, recall Figure 3, and the Wrinkling index value. Analysing these indicators in conjunction with the deformed shapes visualized in Figure 32, we see that the mean amplitude of the wrinkles increases with the height of the cylindrical part. Only the Head-Base section constitutes an exception, which occurs because the wrinkling is so developed in this junction region that some wrinkles tend to roll-over, thus inducing very small amplitude oscillations. This reveals an area of very low physical stability, in which the wrinkling is insufficiently captured by the mesh. The response in terms of the number of wrinkles -see Figure 34a -tends to be more stable. Here also, the number of wrinkles for the Head-Base section constitutes an exception, as it rapidly increases then stabilizes around much higher values than in the other sections.

The general behaviour of the curves for the Wrinkling Index -see Figure 34c -follows the number of wrinkles. The response is in this case slightly more regular, and it respects the general wrinkling behaviour observed in Figure 32.

Discussion

First of all, it should be emphasized that the general deformation pattern of a sail is well represented, in particular with the wrinkling developing for Geometrically, sails made of a spherical head and a cylindrical base are those where the wrinkles determine a large change in the global flying shape. In such configurations large wrinkles arise at the junction between the spherical and the cylindrical parts, and this induces a shrinkage which significantly changes the extension of the structure head. Moreover, we have observed a good convergence of the computed wrinkles when adequately refining the meshes, except in the case of the Spinnaker-like shapes with a slender cylindrical body. In particular, for the case considered in Section 5 namely, with h/R = 1 -we have obtained a significant mesh dependence and a lack of resolution of the wrinkles including for the finest meshes considered, recall Table 7 and Figure 25.

In order to further analyse this case we display in Figure 35 the details of the deformed pattern by plotting the distance to the regression curve against the arc-length coordinate, for two of the finest meshes already considered. We see that the numerical solutions consistently identify two main groups of wrinkles, albeit that within each group the individual wrinkles are insufficiently resolved, with just a few nodes per wrinkle even for the finest mesh. This indicates a rather low physical stability -locally -in this configuration due to the shape considered. Moreover, as the wrinkling wavelength is consistently estimated with the three meshes, we can conjecture that increasing the mesh density by a factor 3-5 across the wrinkles would produce an accurate prediction of the phenomenon. Of course, refining the mesh uniformly and in all directions by such a factor would lead to inordinately expensive computations, hence the mesh adaptation should be performed locally and in a non-isotropic manner, which is out of the scope of our analysis.

This phenomenon is even more pronounced for Spinnaker shapes with a higher cylindrical body, as exemplified in Figure 36. In this case, some rather deep and rapidly oscillating wrinkles appear around the middle of the section, and here clearly much more refined meshes would be required to resolve these features, probably also with self-contact analysis -and possibly bifurcation analysis -involved.

Concluding remarks

A dynamic Finite Element method has been used for the analysis of downwind sail-type structures. The main purpose was here to investigate the development of wrinkling, a phenomenon commonly observed in practice 35 The analysis was carried out using non-linear MITC shell finite elements. The main advantage of using a shell model -instead of a membrane model -is that no ad hoc wrinkling model is required, as wrinkles naturally appear on the deformed surface and are controlled by a balance of membrane and bending stresses. The solution was obtained using a dynamic scheme with Rayleigh damping, able to provide sufficient stability in the numerical analysis. Instability can in fact arise due to the effect of the fabric thickness, which is usually very small compared to the other dimensions of the structure.

A validation test case was reproduced and excellent agreement was found between the numerical and the experimental results. Compared to previous numerical solutions the present method has the advantage of not requiring any particular initialization.

Subsequently, realistic sail shapes were considered, with a prescribed pressure loading. In most cases we found that both the overall deformed configuration and the detailed wrinkling pattern could be adequately captured with meshes of reasonable sizes. For just a few geometries, however -corresponding to Spinnaker shapes with a long cylindrical body -we concluded from the detailed convergence analysis that, while the global flying shape was well reproduced, some very fine and highly oscillating wrinkling behaviour was clearly detected but would require a much denser mesh to resolve the phenomenon. For such computations, dynamically adaptive meshes and/or highly intensive computing would be needed, which was out of the scope of this paper.

Finally, we proposed an enhanced Wrinkling Index which allows to detect the general behaviour of a deformed section in terms of the wrinkling development. This index was shown to be a particularly valuable indicator for the mesh sensitivity and the parametric shape variation analysis.

Further developments primarily regard various refinements of the model, in order to represent some even more realistic sail configurations such as with non-isotropic material models, corner reinforced zones, and cable boundary conditions. Of course, another very important perspective concerns the use of such sail models coupled with the wind flow in a fluid-structure interaction framework. Whereas in this work a constant pressure loading has been used to represent the operating condition of the sail, in reality there is a complex interaction between the deformed shape -including wrinkles -and the distribution of the actual wind loading [START_REF] Trimarchi | The use of shell elements to capture sail wrinkles, and their influence on aerodynamic loads[END_REF]. The strategy developed in this work provides a basis for coupled FSI analysis [START_REF] Trimarchi | Analysis of downwind sail structures using non-linear shells finite elements -wrinkle development and fluid interaction effects[END_REF]. Of course, for an actual transient FSI analysis, the structural damping should be modeled and calibrated in a physical manner, and an adequate coupling strategy should be used to ensure the numerical stability of the coupled solution, see e.g. [START_REF] Lombardi | Numerical simulation of sailing boats: Dynamics, FSI, and shape optimization[END_REF][START_REF] Lombardi | A strongly coupled fluid-structure interaction model for wind-sail simulation[END_REF].

Figure 1 :

 1 Figure 1: Examples of wrinkled downwind sails -a: Courtesy of Christian Février and b: Courtesy of Bravura Yachts

2 SolveFigure 2 :

 22 Figure 2: Flow chart summarizing the analysis: eigen-analysis, non-linear finite element solution and time-step adaptation

Figure 3 :

 3 Figure 3: Wrinkling recognition criterion

Figure 4 :

 4 Figure 4: Snapshot of membrane in shear (taken from [29], Fig. 3)

Figure 5 :

 5 Figure 5: Experimental results -a: Measured deformed section with 3 mm imposed displacement (reproduced from [29], Fig. 4) and b: FFT analysis of the deflection

Figure 6 :

 6 Figure 6: Mesh sensitivity for the validation test case -a: Number of wrinkles n w and wrinkling index I w ; b: Amplitude A w and standard deviation σ

Figure 7 :Figure 8 :

 78 Figure 7: a: Adapted mesh with 210x50 elements (symmetry); b: Deformed shape for 3 mm displacement

Figure 9 :Figure 10 :

 910 Figure 9: Comparison of (a) experimental and (b) numerical deformed shapes

Figure 11 :

 11 Figure 11: Effect of the linear constraint on the deformed shape of the membrane in shear -left: Numerical/experimental comparison; b: FFT of the deflections

Figure 12 :

 12 Figure 12: Top: Perspective view of the wrinkle pattern for the Abaqus result (taken from [29], Fig.6); Bottom: Wrinkle pattern for the MITC4 result (real displacements magnified by factor 2.5)

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure13: Gennaker-like geometry considered in the stability analysis, see Table5for detailed specifications

Figure 17 :Figure 18 :

 1718 Figure 17: Definitions used for identifying the Head and Corner sections

Figure 19 :Figure 20 :Figure 21 :

 192021 Figure 19: Gennaker shape -a: Mesh impact on point-wise displacements. b: Mesh sensitivity in the number of wrinkles

Figure 22 :Figure 23 :Figure 24 :

 222324 Figure 22: Spinnaker geometry with the trace of the sections used for the evaluations

Figure 25 :

 25 Figure 25: Spinnaker shape mesh sensitivity -Zoom on the wrinkles of the Head section, as specified by the arrows in Figure 24. Top: Region near the symmetry centerline. Bottom: Region near the half-side of the structure

Figure 26 :

 26 Figure 26: Parametric analysis for Gennaker shape -a: Head and b: Corner deformed sections for test cases referenced in Table 8

Figure 27 : 2 -Figure 28 :

 27228 Figure 27: Parametric analysis for Gennaker shape -a: Number of wrinkles and b: Average amplitude of the wrinkles

Figure 29 :

 29 Figure 29: Five different analysed geometries (symmetry) with the trace of the sections used for the analysis

Figure 30 :

 30 Figure 30: Gennaker shape with varying base angle -a: Head and b: Corner deformed sections with corresponding regressions

Figure 31 :

 31 Figure 31: Gennaker shape with varying base angle -a: Number, b: Average amplitude and c: Wrinkling index of wrinkles for the Head and the Corner sections

Figure 32 :Figure 33 :

 3233 Figure 32: Spinnaker-like deformed shapes with varying height parameter h/R = 0.5; 1; 1.5; 2 (left to right)

Figure 34 :

 34 Figure 34: Spinnaker shape with varying cylinder length -a: Number, b: Average amplitude and c: Wrinkling Index of the identified wrinkles

Figure 35 :

 35 Figure 35: Spinnaker shape with h/R = 1: Head cross-section in arc-length representation for 2 fine meshes

Figure 36 :

 36 Figure36: Spinnaker shape with h/R = 1.5: Head cross-section in arc-length representation for computational mesh

Table 1 :

 1 Kapton R physical quantities reported by[START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF]

Table 2 :

 2 MITC4 results for series of meshes considered

	Subdivisions Elements n w W centre-line [mm]	σ	I w
	65x50	3038	13	0.47	0.250 3.25
	105x50	4998	15	0.47	0.362 5.43
	160x50	7693	19	0.40	0.232 4.41
	210x50	10143	20	0.38	0.178 3.56
	420x75	28638	19	0.40	0.136 2.58
	420x100	42183	21	0.35	0.155 3.23
	630x150	55803	20	0.37	0.140 2.8
	880x150	130673	20	0.37	0.140 2.8

Table 3 :

 3 Abaqus results for membrane in shear (taken from[START_REF] Wong | Wrinkled membranes: Experiments (Part I), analytical models (Part II) and numerical simulations (Part III)[END_REF], Table4)

	N. elements N wrinkles Abaqus
	3960	16+2
	6950	19+2
	13134	19+2

Table 4 :

 4 Equilibrium and eigenvalues of the tangent stiffness matrix

	Matrix property	Eigenvalues	Equilibrium
	Positive definite	λ i > 0 ∀i = 1, ndof	Stable
	Non-positive definite	λ j < 0 for some j	Unstable
	Semi-positive	λ i 0 ∀i and λ j = 0 for some j	Critical

Table 5 :

 5 Parameters used in the eigenvalue analysis

	Radius Width angle N. elts.	t	E	ν	ρ s
	[mm]	[Rad]	[-]	[mm] [N/mm 2 ] [-] [Kg/m 3 ]
	1000	π/5	8000	0.05	376	0.4	1.15

Table 6 :

 6 Meshes specifications for the Gennaker mesh sensitivity analysis. N e is the number of elements, N p the number of nodes and EpS the number of element per horizontal section

	Mesh	m 2	m 3	m 5	m 10	m 15	m 20	m 25	m 30
	N e	49	130	406	1761	4066	7321 11526 16681
	N p	67	158	405	1859	4214	7519 11774 16979
	EpS	4	7	13	28	43	58	73	88
	Mesh m 35	m 40	m 45	m 50	m 55	m 60
	N e 22786 29841 37846 46801 56706 67561
	N p 23134 30239 38294 47299 57254 68159
	EpS	103	118	133	148	163	178

Table 7 :

 7 Mesh specification for the Spinnaker mesh sensitivity analysis. N e is the number of elements, N p the number of nodes and EpS the number of element per horizontal section

	Mesh m 3	m 5	m 10	m 15	m 20	m 25	m 30
	N e	24	104	513	1024 22203 3552 5133
	N p	36	128	565	1283	2330	3687 5250
	EpS	4	8	18	28	38	48	58
	Mesh m 35 m 40	m 45	m 50	m 55	m 60
	N e 7072 9477 12848 16023 20088 23423
	N p 7262 9634 13105 16310 20410 23770
	EpS	68	78	88	98	108	118

Table 8 :

 8 Material test matrix (Analysis code number)

	`````````````È	[N/mm 2 ]	t [mm] 0.05 0.1 0.2
		188	1	2	3
		376	4	5	6
		752	7	8	9

Table 8

 8

	Number of wrinkles				Wrinkle average amplitude
					N w Corner					A w Corner
					N w Head					A w Head
	0 1 2 3 4 5 6 7 8							0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
	0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 t [mm]	800	700	600	300 E [N/mm 400 500 2 ]	200	100	0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 t [mm]	800	700	600	300 E [N/mm 400 500 2 ]	200	100

Table 9 :

 9 Specifications of the analysed geometries

The transverse direction vector is not normal to the midsurface in general deformed configurations


O Osman

C Vézy

J Pivetal

M Frénea-Robin

N Haddour

F Buret

L F Zanini

G Reyne

N M Dempsey

F Dumas-Bouchiat

A NOVEL DEVICE FOR CONTINUOUS FLOW MAGNETIC TRAPPING AND SORTING OF HUMAN CELLS USING FLAT MICRO-PATTERNED NdFeB FILMS

Keywords: Magnetism, Cell separation, Microfluidic

The manipulation of magnetically labeled bio-objects of nano or micrometer sizes is now realizable by combining the magnetic forces with microfluidics. This paper reports a novel method for manipulating, trapping and sorting magnetic beads and magnetically labeled cells under continuous flow. This method uses a microfluidic device with integrated flat micropatterned hard magnetic films which generate high magnetic field gradients.

INTRODUCTION

Magnetic field applications at the microscale are of growing interest in biology [START_REF] Pamme | Continuous sorting of magnetic cells via on-chip free-flow magnetophoresis[END_REF], as magnetic field sources can now be integrated directly into microchips [START_REF] Smistrup | Microelectromagnet for magnetic manipulation in lab-on-a-chip systems[END_REF][START_REF] Choi | An on-chip magnetic bead separator for biocell sorting[END_REF]. Due to the recent progress in the synthesis of nano and micro-superparamagnetic particles which bound to biological objects, a large number of biotechnological devices based on magnetophoresis have been developed. Nevertheless, the generation of magnetic fields and field gradients strong enough to manipulate magnetic micro objects remains a challenge. Several papers report the use of micro-electromagnets or the combination of soft materials and external magnetic field to create magnetic field gradients for manipulating biological objects [START_REF] Derec | Local control of magnetic objects in microfluidic channels[END_REF]. The use of permanentmagnet-based microstructures is very attractive since Joule heating is avoided while no energy source is required.

RESULTS AND DISCUSSION

In our experiments, we used flat chessboard micro-magnet arrays (figure "1") obtained by thermo-magnetic patterning of hard magnetic neodymium iron boron (NdFeB) thick film [START_REF] Dumas-Bouchiat | Thermomagnetically patterned micromagnets[END_REF]. Such patterned films have already been used to trap liposomes containing superparamagnetic nanoparticles [START_REF] Pivetal | Trapping of magnetically-labeled liposomes on flat micro-patterned hard magnetic films[END_REF]. Figure "2" shows HEK 293 cells, which have endocyted 100 nm-magnetic nanoparticles, attracted towards the regions of maximum magnetic field gradient (edges of the squares). This result reveals the applicability of micro-magnets for trapping biological cell and for potential applications in medical diagnostics based on cell endocytic ability. We also demonstrated the efficient trapping of Jurkat cells labeled via the binding of cell antigens with antibody-coated magnetic particles of 1 µm (figure "3"). Due to the flatness of the micro-magnet array, this type of autonomous micro magnetic field source can be easily integrated into a microchannel. The micro-magnet array was first covered with a thin layer of PolyDiMethylSiloxane (PDMS) obtained by spin-coating, the thickness (2 to 12 µm) being controlled by mixing PDMS with heptane to decrease its viscosity and by varying the ratio between both. A PDMS microchannel was sealed onto the PDMS layer using an O 2 plasma (figure "4"). The PDMS layer serves to protect the micro-magnet array from degradation due to contact with liquid and offers a biocompatible interface with the manipulated bio-objects. Moreover, by adjusting the thickness of the PDMS layer, the magnetic force exerted on the particles can be varied. To demonstrate the trapping efficiency of our microfluidic magnetic device, two types of 1 µm-fluorescent beads (magnetic and non-magnetic) were injected. Figure "5" shows magnetic beads trapped at the regions of maximum field gradient while non-magnetic beads moved freely.

CONCLUSION

The lab on chip combination of micro-patterned magnetic films and microfluidic is convenient for trapping and sorting bio-objects. As the microscale patterning of the hard magnetic films generates strong magnetic field gradients (up to 10 6 T/m), cells can be efficiently trapped with reduced magnetic loads. Future experiments will be directed towards comparing our device with commercial kits (eg. Dynabeads ® FlowComp TM ) enabling selective isolation of labeled cells using macro-magnets. As a first step, we succeeded in trapping cells labeled with 30 nm-magnetic nanoparticles using micro-patterned hard magnetic films (figure "6"). Note that due to their very small size, such nanoparticles cannot be attracted using macro-magnets.

*Figure 1 :

 1 Figure 1: (a) Structure of the thermomagnetically patterned micromagnet array -(b) Magneto-optic image (Faraday effect) of a MOIF placed on top of a patterned NdFeB films -(c) Norm of the magnetic flux density in a plane situated at the magnetic film surface

Figure 2 :

 2 Figure 2: Attraction of HEK 293 towards magnetic field gradient maxima. (a) unlabeled cells. (b) Magnetically labeled cells. (c) Magnetically and fluorescently labeled cells. (d) zoom on a group of attracted cells.

Figure 3 :

 3 Figure 3: Attraction of magnetically (a) and fluorescently (b) labeled Jurkat cells towards magnetic field gradient maxima. HEK 293 (unlabeled cells) remained untrapped.

Figure 4 :

 4 Figure 4: schematic view of the microfluidic device.

Figure 5 :

 5 Figure 5: (a) Fluorescence image of magnetic beads trapped at the maximum of the gradient field. (b) Fluorescence image of the non-magnetic particles (bright spots indicated by arrows), moving freely in the microchannel.

Figure 6 :

 6 Figure 6: (a) Overview of the Dynabeads ® FlowComp TM isolation procedure.(b) Attraction of HEK 293 cells labeled with 30 nm-nanoparticles toward magnetic field gradient maxima.