Sensitivity Analysis of the Orthoglide: A

Keywords: parallel kinematic machine, sensitivity analysis, kinematic analysis, kinematic singularity, isotropy

Three-DOF Translational Parallel Kinematic Machine

In this paper, two complementary methods are introduced to analyze the sensitivity of a three-degree-of-freedom (3-DOF) translational parallel kinematic machine (PKM) with orthogonal linear joints: the Orthoglide. Although these methods are applied to a particular PKM, they can be readily applied to 3-DOF Delta-Linear PKM such as ones with their linear joints parallel instead of orthogonal. On the one hand, a linkage kinematic analysis method is proposed to have a rough idea of the influence of the length variations of the manipulator on the location of its end-effector. On the other hand, a differential vector method is used to study the influence of the length and angular variations in the parts of the manipulator on the position and orientation of its end-effector. Besides, this method takes into account the variations in the parallelograms. It turns out that variations in the design parameters of the same type from one leg to another have the same effect on the position of the end-effector. Moreover, the sensitivity of its pose to geometric variations is a minimum in the kinematic isotropic configuration of the manipulator. On the contrary, this sensitivity approaches its maximum close to the kinematic singular configurations of the manipulator.

Introduction

For two decades, parallel manipulators have attracted the attention of more and more researchers who consider them as valuable alternative design for robotic mechanisms. As stated by numerous authors, conventional serial kinematic machines have already reached their dynamic performance limits, which are bounded by high stiffness of the machine components required to support sequential joints, links, and actuators. Thus, while having good operating characteristics ͑large workspace, high flexibility and ma-neuverability͒, serial manipulators have disadvantages of low stiffness and low power. Conversely, parallel kinematic machines ͑PKMs͒ offer essential advantages over their serial counterparts ͑lower moving masses, higher stiffness and payload-to-weight ratio, higher natural frequencies, better accuracy, simpler modular mechanical construction, possibility to locate actuators on the fixed base͒.

However, PKMs are not necessarily more accurate than their serial counterparts. Indeed, even if the dimensional variations can be compensated with PKM, they can also be amplified contrary to with their serial counterparts ͓1͔. Wang and Masory ͓2͔ studied the effect of manufacturing tolerances on the accuracy of a Stewart platform. Kim and Choi ͓3͔ used a forward error bound analysis to find the error bound of the end-effector of a Stewart platform when the error bounds of the joints are given, and an inverse error bound analysis to determine those of the joints for the given error bound of the end-effector. Kim and Tsai ͓4͔ studied the effect of misalignment of linear actuators of a 3-DOF translational parallel manipulator on the motion of its moving platform. Han et al. ͓5͔ used a kinematic sensitivity analysis method to explain the gross motions of a 3-UPU parallel mechanism, and they showed that it is highly sensitive to certain minute clearances. Fan et al. ͓6͔ analyzed the sensitivity of the 3-PRS parallel kinematic spindle platform of a serial-parallel machine tool. Verner et al. ͓7͔ presented a new method for optimal calibration of PKM based on the exploitation of the least error sensitive regions in their workspace and geometric parameters space. As a matter of fact, they used a Monte Carlo simulation to determine and map the sensitivities to geometric parameters. Moreover, Caro et al. ͓8͔ developed a tolerance synthesis method for mechanisms based on a robust design approach.

This paper aims at analyzing the sensitivity of the Orthoglide to its dimensional and angular variations. The Orthoglide is a 3-DOF translational PKM developed by Chablat and Wenger ͓9͔. A smallscale prototype of this manipulator was built at IRCCyN.

Here, the sensitivity of the Orthoglide is studied by means of two complementary methods. First, a linkage kinematic analysis is used to have a rough idea of the influence of the dimensional variations to its end-effector and to show that the variations in design parameters of the same type from one leg to another have the same influence on the location of the end-effector. Although this method is compact, it cannot be used to know the influence of the variations in the parallelograms. Thus, a differential vector method is developed to study the influence of the dimensional and angular variations in the parts of the manipulator, and particularly variations in the parallelograms, on the position and the orientation of its end-effector.

In the isotropic kinematic configuration, the end-effector of the manipulator is located at the intersection between the directions of its three actuated prismatic joints, and the condition number of its kinematic Jacobian matrix is equal to 1 ͓10͔. It is shown that this configuration is the least sensitive one to geometrical variations, contrary to the closest configurations to its kinematic singular configurations, which are the most sensitive to geometrical variations.

Although the two sensitivity analysis methods are applied to a

Stéphane Caro Philippe Wnger

Fouad Bennis Damien Chablat particular PKM, these methods can be readily applied to other 3-DOF Delta-linear PKM such as ones with parallel linear joints instead of orthogonal ones.

Manipulator Geometry

The kinematic architecture of the Orthoglide is shown in Fig. 1. It consists of three identical parallel chains that are formally described as PRP a RR, where P, R, and P a denote the prismatic, revolute, and parallelogram joints, respectively, as shown in Fig. 2.

The mechanism input is made up of three actuated orthogonal prismatic joints. The output body ͑with a tool mounting flange͒ is connected to the prismatic joints through a set of three kinematic chains. Inside each chain, one parallelogram is used and oriented in a manner that the output body is restricted to translational movements only.

The small-scale prototype of the Orthoglide was designed to reach Cartesian velocity of 1.2 m / s and an acceleration of 17 m/s 2 . The desired payload is 4 kg ͑spindle, tool, included͒. The size of its prescribed cubic workspace C u is 200ϫ 200 ϫ 200 mm, where the velocity transmission factors are bounded between 1 / 2 and 2. The three legs are supposed to be identical. According to ͓9͔, the nominal lengths L i and widths d i of the parallelograms, and the nominal distances, r i , between points C i and the end-effector P are identical, i.e.,

L = L 1 = L 2 = L 3 = 310.58 mm, d = d 1 = d 2 = d 3 = 80 mm, r = r 1 = r 2 = r 3 = 31 mm.
As depicted in Fig. 3, Q 1 and Q 2 , vertices of C u , are defined at the intersection between the Cartesian workspace boundary and the axis x = y = z expressed in the reference coordinate frame R b . Q 1 and Q 2 are the closest points to the singularity surfaces. Their Cartesian coordinates, expressed in R b , are equal to ͑-73.21, -73.21, -73.21͒ and ͑126.79, 126.79, 126.79͒, respectively. The parts of the manipulator are supposed to be rigid bodies and there is no joint clearance. The legs of the manipulator, composed of one prismatic joint, one parallelogram, and three revolute joints, generate a five-DOF motion each. Besides, they are iden-tical. Therefore, according to Karouia et al. ͓11͔, the manipulator is isostatic. Thus, the results obtained by the sensitivity analysis methods developed in this paper are meaningful.

Sensitivity Analysis

Two complementary methods are used to study the sensitivity of the Orthoglide. First, a linkage kinematic analysis is used to have a rough idea of the influence of the dimensional variations to its end-effector. Although this method is compact, it cannot be used to know the influence of the variations in the parallelograms. Thus, a differential vector method is used to study the influence of the dimensional and angular variations in the parts of the manipulator, and particularly variations in the parallelograms, on the position and the orientation of its end-effector.

Linkage Kinematic Analysis.

This method aims at computing the sensitivity coefficients of the position of the endeffector, P, to the design parameters of the manipulator. First, three implicit functions depicting the kinematic of the manipulator are obtained. A relation between the variations in the position of P and the variations in the design parameters follows from these functions. Finally a sensitivity matrix, which gathers the sensitivity coefficients of P, follows from the previous relation written in matrix form.

3.1.1 Formulation. Figure 1 depicts the design parameters taken into account. Points A 1 , A 2 , and A 3 are the bases of the prismatic joints. Their Cartesian coordinates, expressed in R b , are a 1 , a 2 , and a 3 , respectively.

a 1 = ͓-a 1 0 0͔ T ͑1a͒ a 2 = ͓0 -a 2 0͔ T ͑1b͒ a 3 = ͓0 0 -a 3 ͔ T ͑1c͒
where a i is the distance between points A i and O, the origin of R b .

Points B 1 , B 2 , and B 3 are the links between the prismatic and parallelogram joints. Their Cartesian coordinates, expressed in R b are

b 1 = ΄ -a 1 + 1 b 1y b 1z ΅ ͑2a͒ b 2 = ΄ b 2x -a 2 + 2 b 2z ΅ ͑2b͒
c 1 = ͓p x -r 1 0 0͔ T ͑3a͒ c 2 = ͓0 p y -r 2 0͔ T ͑3b͒ c 3 = ͓0 0 p z -r 3 ͔ T ͑3c͒
where p = ͓p x p y p z ͔ T is the vector of the Cartesian coordinates of the end-effector P, expressed in R b .

The expressions of the nominal lengths of the parallelograms follow from Eq. ͑4͒,

L i = ʈc i -b i ʈ 2 , i = 1,2,3 ͑4͒
where L i is the nominal length of the ith parallelogram and ʈ • ʈ 2 is the Euclidean norm. Three implicit functions follow from Eq. ͑4͒ and are given by the following equations:

F 1 = ͑-r 1 + p x + a 1 -1 ͒ 2 + ͑p y -b 1y ͒ 2 + ͑p z -b 1z ͒ 2 -L 1 2 = 0 F 2 = ͑p x -b 2x ͒ 2 + ͑-r 2 + p y + a 2 -2 ͒ 2 + ͑p z -b 2x ͒ 2 -L 2 2 = 0 F 3 = ͑p x -b 3x ͒ + ͑p y -b 3y ͒ 2 + ͑-r 3 + p z + a 3 -3 ͒ 2 -L 3 2 = 0
By differentiating functions F 1 , F 2 , and F 3 with respect to the design parameters of the manipulator and the position of the endeffector, we obtain a relation between the positioning error of the end-effector ␦p, and the variations in the design parameters ␦q i .

␦F i = A i ␦p + B i ␦q i = 0, i = 1,2,3 ͑5͒
with

A i = ͓ץF i /ץp x ץF i /ץp y ץF i /ץp z ͔ ͑ 6͒ B i = ͓ץF i /ץa i ץF i /ץb iy ץF i /ץb iz ץF i ץ/ i ץF i /ץL i ץF i /ץr i ͔ ͑7͒ ␦p = ͓␦p x ␦p y ␦p z ͔ T ͑8͒ ␦q i = ͓␦a i ␦h i ␦k i ␦ i ␦L i ␦r i ͔ T ͑9͒
where ␦a i , ␦h i , ␦k i , ␦ i , ␦L i , and ␦r i , depict the variations in a i , h i , k i , i , L i , and r i , respectively, with

h 1 = b 1y , k 1 = b 1z , h 2 = b 2x , k 2 = b 2z , h 3 = b 3x , k 3 = b 3y .
Integrating the three loops of Eq. ͑5͒ together and separating the position parameters and design parameters to different sides yields the following simplified matrix form:

A␦p + B␦q = 0 ͑10͒ with A = ͓A 1 T A 2 T A 3 T ͔ T R 3ϫ3 ͑11͒ B = ΄ B 1 0 0 0 B 2 0 0 0 B 3 ΅ R 3ϫ18 ͑12͒ ␦q = ͓␦q 1 T ␦q 2 T ␦q 3 T ͔ T R 18ϫ1 ͑13͒
Equation ͑10͒ takes into account the coupling effect of the three independent structure loops. According to ͓9͔, A is the parallel Jacobian kinematic matrix of the Orthoglide, which does not meet parallel kinematic singularities when its end-effector covers C u . Therefore, A is not singular and its inverse A -1 exists. Thus, the positioning error of the end-effector can be computed using Eq. ͑14͒.

␦p = C␦q

͑14͒

where

C = -A -1 B = ΄ ץp x /ץa 1 ץp x /ץh 1 ¯ץp x /ץr 3 ץp y /ץa 1 ץp y /ץh 1 ¯ץp y /ץr 3 ץp z /ץa 1 ץp z /ץh 1 ¯ץp z /ץr 3 ΅ R 3ϫ18

͑15͒

represents the sensitivity matrix of the manipulator. The terms of C are the sensitivity coefficients of the Cartesian coordinates of the end-effector to the design parameters and are used to analyze the sensitivity of the Orthoglide.

Results of the Linkage Kinematic Analysis.

The sensitivity matrix C of the manipulator depends on the position of its end-effector.

Figures 4567depict the mean of the sensitivity coefficients of p x , p y , p z , and p, when the end-effector covers C u . It appears that the Besides, it is noteworthy that p x ͑p y , p z , respectively͒ is very sensitive to the design parameters which make up the 1st ͑2nd, 3rd, respectively͒ leg of the manipulator, contrary to the others. That is due to the symmetry of the architecture of the manipulator. Henceforth, only the variations in the design parameters of the first leg of the manipulator will be taken into account. Indeed, the sensitivity of the position of the end-effector to the variations in the design parameters of the second and the third legs of the manipulator can be deduced from the sensitivity of the position of the end-effector to variations in the design parameters of the first leg.

Chablat and Wenger ͓9͔ showed that if the prescribed bounds of the velocity transmission factors ͑the kinematic criteria used to dimension the manipulator͒ are satisfied at Q 1 and Q 2 , then these bounds are satisfied throughout the prescribed cubic Cartesian workspace C u . Q 1 and Q 2 are then the most critical points of C u , whereas O is the most interesting point because it corresponds to the isotropic kinematic configuration of the manipulator. Here, we assume that if the prescribed bounds of the sensitivity coefficients are satisfied at Q 1 and Q 2 , then these bounds are satisfied throughout C u .

Figures 8 and9 depict the sensitivity coefficients of p x and p y to the dimensional variations in the 1st leg, i.e., a 1 , b 1y , b 1z , 1 , L 1 , r 1 , along Q 1 Q 2 . It appears that these coefficients are a minimum in the isotropic configuration, i.e., P ϵ O, and a maximum when P ϵ Q 2 , i.e., in the closest configuration to the singular one. Figure 10 depicts the sensitivity coefficients of p along diagonal Q 1 Q 2 . It is noteworthy that all the sensitivity coefficients are a minimum when P ϵ O and a maximum when P ϵ Q 2 . Finally, Fig. 11 Figures 14 and 15 depict the sensitivity coefficients of p x and p when the end-effector hits Q 2 ͑P ϵ Q 2 ͒. In this case, variations in p x , p y , and p z are coupled. For example, variations in p x are due to both dimensional variations in the 1st leg and variations in the 2nd and the 3rd legs. Besides, the amplification of the dimensional variations is important. Indeed, the sensitivity coefficients of p are close to 2 in this configuration. For example, as the sensitivity However, this linkage kinematic method does not take into account variations in the parallelograms, except the variations in their global length. Thus, a differential vector method is developed below.

Differential Vector Method.

In this section, we perfect a sensitivity analysis method of the Orthoglide, which complements the previous one. This method is used to analyze the sensitivity of the position and the orientation of the end-effector to dimensional and angular variations, and particularly to the variations in the parallelograms. Moreover, it allows us to distinguish the variations which are responsible for the position errors of the endeffector from the ones which are responsible for its orientation errors. To develop this method, we were inspired by a Huang & al. work on a parallel kinematic machine, which is made up of parallelogram joints too ͓12͔.

First, we express the dimensional and angular variations in vectorial form. A relation between the position and the orientation errors of the end-effector is then obtained from the closed-loop kinematic equations. The expressions of the orientation and the position errors of the end-effector, with respect to the variations in the design parameters, are deduced from this relation. Finally, we introduce two sensitivity indices to assess the sensitivity of the position and the orientation of the end-effector to dimensional and angular variations, and particularly to the parallelism errors of the bars of the parallelograms.

Formulation.

The schematic drawing of the ith leg of the Orthoglide depicted in Fig. 2 is split in order to depict the variations in design parameters in a vectorial form. The closedloop kinematic chains

O -A i -B i -B ij -C ij -C i -P, i =1,2,3, j =1,
a i -o = R i ͑a 0 + ␦a i ͒ ͑ 16͒
where a 0 is the nominal position vector of A i with respect to O expressed in R i , ␦a i is the positioning error of A i . R i is the transformation matrix from R i to R b . I 3 is the ͑3 ϫ 3͒ identity matrix and

R 1 = I 3 ͑17͒ R 2 = ΄ 0 0 -1 1 0 0 0 -1 0 ΅ ͑18͒ R 3 = ΄ 0 1 0 0 0 1 1 0 0 ΅ ͑19͒ According to Fig. 17, b i -a i = R i ͑ i + ␦ i ͒e 1 + R i ␦ Ai ϫ ͑ i + ␦ i ͒e 1

͑20͒

where i is the displacement of the ith prismatic joint, ␦ i is its displacement error, ␦ Ai = ͓␦ Aix ␦ Aiy ␦ Aiz ͔ T is the angular variation of its direction, and

e 1 = ΄ 1 0 0 ΅ ͑21͒ e 2 = ΄ 0 0 1 ΅ ͑22͒ ͑j͒ = ͭ 1, if j = 1 -1, if j = 2 ͮ ͑23͒ According to Fig. 18, b ij -b i = R i ͓I 3 + ␦ Ai ϫ ͔͑͑j͒͑d/2 + ␦b i /2͓͒I 3 + ␦ Bi ϫ ͔e 2 ͒ ͑24͒ c ij -b ij = L i w i + ␦L i w i + L i ␦w i ͑25͒
where d is the nominal width of the parallelogram, ␦b i is the variation in the length of link B i1 B i2 and is supposed to be equally shared by each side of B i . ␦ Bi = ͓␦ Bix ␦ Biy ␦ Biz ͔ T is the orientation error of link B i1 B i2 with respect to the direction of the ith prismatic joint, L i is the length of the ith parallelogram, ␦L ij is the variation in the length of link B ij C ij , of which w i is the direction, and ␦w i is the variation in this direction, orthogonal to w i .

According to Fig. 19,

c ij -c i = R i ͓I 3 + ␦ ϫ ͔͑͑j͒͑d/2 + ␦c i /2͓͒I 3 + ␦ Ci ϫ ͔e 2 ͒ ͑26͒ c i -p = ͓I 3 + ␦ ϫ ͔R i ͑c 0 + ␦c i ͒ ͑ 27͒
where ␦c i is the variation in the length of link C i1 C i2 , which is supposed to be equally shared by each side of C i .

␦ Ci = ͓␦ Cix ␦ Ciy ␦ Ciz ͔ T is the orientation error of link C i1 C i2 with respect to link C i P. c 0 is the nominal position vector of C i with respect to end-effector P, expressed in R i , ␦c i is the position error of C i expressed in R i , and ␦ = ͓␦ x ␦ y ␦ z ͔ T is the orientation error of the end-effector, expressed in R b .

Implementing linearization of Eqs. ͑16͒-͑26͒ and removing the components associated with the nominal constrained equation p 0 = R i ͑a 0 + i e 1 -c 0 ͒ + L i w i , yields

␦p = p -p 0 = R i ͑␦e i + i ͑␦ Ai ϫ e 1 ͒ + ͑j͒d/2͑␦ Ai ϫ e 2 ͒ + ͑j͒d/2͑␦␥ i ϫ e 2 ͒ + ͑j͒␦m i /2e 2 ͒ + ␦L ij w i + L i ␦w i -␦ ϫ R i ͑c 0 + d/2͑j͒e 2 ͒ ͑ 28͒
where ␦p is the position error of the end-effector of the manipulator; ␦e i = ␦a i + ␦ i e 1 -␦c i is the sum of the position errors of points A i , B i , and Equation ͑28͒ shows the coupling of the position and orientation errors of the end-effector. Contrary to the orientation error, the position error can be compensated because the manipulator is a translational 3-DOF PKM. Thus, it is more important to minimize the geometrical variations, which are responsible for the orientation errors of the end-effector than the ones, which are responsible for its position errors.

C i expressed in R i ; ␦␥ i = ␦ Bi -␦ Ci is
The following equation is obtained by multiplying both sides of Eq. ͑28͒ by w i T and utilizing the circularity of hybrid product.

w i T ␦p = w i T R i ␦e i + i ͑R i e 1 ϫ w i ͒ T R i ␦ Ai + ͑j͒d/2͑R i e 2 ϫ w i ͒ T R i ͑␦ Ai + ␦␥ i ͒ + ͑j͒␦m i /2w i T R i e 2 + ␦L ij -͕R i ͓c 0 + ͑j͒d/2e 2 ͔ ϫ w i ͖ T ␦ ͑29͒ 3.2.
1.1 Orientation Error Mapping Function. By subtraction of Eqs. ͑29͒ written for j = 1 and j = 2, and for the ith kinematic chain, a relation is obtained between the orientation error of the end-effector and the variations in design parameters, which is independent of the position error of the end-effector.

d͑R i e 2 ϫ w i ͒ T ␦ = ␦l i + d͑R i e 2 ϫ w i ͒ T R i ͑␦ Ai + ␦␥ i ͒ + ␦m i w i T R i e 2

͑30͒

where ␦l i = ␦L i1 -␦L i2 , the relative length error of links B i1 C i1 and B i2 C i2 , depicts the parallelism error of links B i1 B i2 and C i1 C i2 as shown in Fig. 20. Equation ͑30͒ can be written in matrix form:

␦ = J ͑31͒ with J = D -1 E ͑32͒ D = d ΄ ͑R 1 e 2 ϫ w 1 ͒ T ͑R 2 e 2 ϫ w 2 ͒ T ͑R 3 e 2 ϫ w 3 ͒ T ΅ ͑33͒ E = ΄ E 1 ¯Ē 2 Ē3 ΅ ͑34͒ E i = ͓1 w i T R i e 2 d͑R i e 2 ϫ w i ͒ T R i d͑R i e 2 ϫ w i ͒ T R i ͔ ͑35͒
␦ is the orientation error of the end-effector expressed in R b , and

= ͑ 1 T , 2 T , 3 T ͒ T such that i = ͑␦l i , ␦m i , ␦ Ai T , ␦␥ i T ͒ T .
The determinant of D will be null if the normal vectors to the plans, which contain the three parallelograms respectively, are collinear, or if one parallelogram is flat. Here, this determinant is not null when P covers C u because of the geometry of the manipulator. Therefore, D is nonsingular and its inverse D -1 exists.

As ͑R i e 2 ϫ w i ͒ T Ќ R i e 2 , ␦ Aiz , and ␦␥ iz , the third components of ␦ Ai and ␦␥ i expressed in R i , have no effect on the orientation of the end-effector. Thus, matrix J can be simplified by eliminating its columns associated with ␦ Aiz and ␦␥ iz , i =1,2,3. Finally, eighteen variations: ␦l i , ␦m i , ␦ Aix , ␦ Aiy , ␦␥ ix , ␦␥ iy , i =1,2,3, should be responsible for the orientation error of the end-effector.

Position Error Mapping Function. By addition of Eqs.

͑29͒ written for j = 1 and j = 2, and for the ith kinematic chain, a relation is obtained between the position error of the end-effector and the variations in design parameters, which does not depend on ␦␥ i .

w i T ␦p = ␦L i + w i T R i ␦e i + i ͑R i e 1 ϫ w i ͒ T R i ␦ Ai -͑R i c 0 ϫ w i ͒ T ␦ ͑36͒
Equation ͑36͒ can be written in matrix form:

␦p = J pp p + J p = ͓J pp J p ͔ͫ p ͬ ͑37͒ with J pp = F -1 G ͑38͒ J p = F -1 HJ ͑39͒ F = ͓w 1 w 2 w 3 ͔ T ͑40͒ G = diag͑G i ͒ ͑ 41͒ G i = ͓1 w i T R i i ͑R i e 1 ϫ w i ͒ T R i ͔ ͑ 42͒ H = -͓R 1 c 0 ϫ w 1 R 2 c 0 ϫ w 2 R 3 c 0 ϫ w 3 ͔ ͑ 43͒ p = ͑ p1 T , p2 T , p3 T ͒ T ͑44͒ pi = ͑␦L i ,␦e i T ,␦ Ai T ͒ T ͑45͒ ␦L i = ͑␦L i1 + ␦L i2 ͒ / 2

is the mean value of the variations in links

B i1 C i1 and B i2 C i2 , i.e., the variation in the length of the ith parallelogram. p is the set of the variations in design parameters, which should be responsible for the position errors of the endeffector, except the ones which should be responsible for its orientation errors, i.e., . p is made up of three kinds of errors: the variation in the length of the ith parallelogram ͑i.e., ␦L i , i =1,2,3͒, the position errors of points A i , B i , and C i , ͑i.e., ␦e i , i =1,2,3͒ and the orientation errors of the directions of the prismatic joints ͑i.e., ␦ Ai , i =1,2,3͒. Besides, F is nonsingular and its inverse F -1 exists because F corresponds to the Jacobian kinematic matrix of the manipulator, which is not singular when P covers C u , ͓9͔.

According to Eq. ͑36͒ and as ͑R i e 1 ϫ w i ͒ T Ќ R i e 1 , matrix J pp can be simplified by eliminating its columns associated with ␦ Aix , i =1,2,3. Finally, 33 variations: ␦L i , ␦e ix , ␦e iy , ␦e iz , ␦ Aix , ␦ Aiy , ␦ Aiz , ␦l i , ␦m i , ␦␥ ix , ␦␥ iy , i =1,2,3, should be responsible for the position error of the end-effector. Rearranging matrices J pp and J p , the position error of the endeffector can be expressed as:

␦p = J q = ͓J 1 J 2 J 3 ͔͑ q1 q2 q3 ͒ T ͑46͒ with qi = ͑␦L i , ␦e ix , ␦e iy , ␦e iz , ␦ Aix , ␦ Aiy , ␦ Aiz , ␦l i , ␦m i ,
␦␥ ix , ␦␥ iy ͒, and J R 3ϫ33 .

Sensitivity Indices.

In order to investigate the influence of the design parameters errors on the position and the orientation of the end-effector, sensitivity indices are required. According to Sec. 3.1.2, variations in the design parameters of the same type from one leg to the other have the same influence on the location of the end-effector. Thus, assuming that variations in the design parameters are independent, the sensitivity of the position of the end-effector to the variations in the kth design parameter responsible for its position error, i.e., q͑1,2,3͒k , is called k and is defined by Eq. ͑47͒.

k = ͱ ͚ i=1 3 ͚ m=1 3 J imk 2 , k = 1, . . . ,11 ͑47͒
Likewise, r is a sensitivity index of the orientation of the endeffector to the variations in the rth design parameter responsible for its orientation error, i.e., ͑1,2,3͒r . r follows from; Eq. ͑31͒ and is defined by Eq. ͑48͒.

r = arccos tr͑Q r ͒ -1 2 ͑48͒
where Q r is the rotation matrix corresponding to the orientation error of the end-effector, and r is a linear invariant: its global rotation ͓13͔.

Q r = ΄ C zr C
xr = ͱ ͚ j=0 2 J 1͑6j+r͒ 2 , r = 1, . . . ,6 ͑50͒ yr = ͱ ͚ j=0 2 J 2͑6j+r͒ 2 , r = 1, . . . ,6 ͑51͒ zr = ͱ ͚ j=0 2 J 3͑6j+r͒ 2 , r = 1, . . . ,6 ͑52͒
Finally, k can be employed as a sensitivity index of the position of the end-effector to the kth design parameter responsible for the position error. Likewise, r can be employed as a sensitivity index of the orientation of the end-effector to the rth design parameter responsible for the orientation error. It is noteworthy that these sensitivity indices depend on the location of the endeffector.

Results of the Differential Vector

Method. The sensitivity indices defined by Eqs. ͑47͒ and ͑48͒ are used to evaluate the sensitivity of the position and orientation of the end-effector to variations in design parameters, particularly to variations in the parallelograms.

Figures 21͑a͒ and21͑b͒ depict the sensitivity of the position of the end-effector along the diagonal Q 1 Q 2 of C u , to dimensional variations and angular variations, respectively. According to Fig. 21͑a͒, the position of the end-effector is very sensitive to variations in the lengths of the parallelograms, ␦L i , and to the position errors of points A i , B i , and C i along axis x i of R i , i.e., ␦e ix . Conversely, the influence of ␦l i and ␦m i , the parallelism errors of the parallelograms, is low and even negligible in the kinematic isotropic configuration. According to Fig. 21͑b͒, the orientation errors of the prismatic joints depicted by ␦ Aiy and ␦ Aiz are the most influential angular errors on the position of the end-effector. Besides, the position of the end-effector is not sensitive to angular variations in the isotropic configuration.

Figures 22͑a͒ and22͑b͒ depict the sensitivity of the orientation of the end effector, along Q 1 Q 2 , to dimensional and angular variations. According to Fig. 22͑a͒, ␦l i and ␦m i are the only dimensional variations, which are responsible for the orientation error of the end-effector. However, the influence of the parallelism error of the small sides of the parallelograms, depicted by ␦l i , is more important than the one of the parallelism error of their long sides, depicted by ␦m i .

According to Figs. 21 and 22, the sensitivity of the position and the orientation of the end-effector is generally null in the kinematic isotropic configuration ͑p =0͒, and is a maximum when the manipulator is close to a kinematic singular configuration, i.e., P ϵ Q 2 . Indeed, only two kinds of design parameters variations are responsible for the variations in the position of the endeffector in the isotropic configuration: ␦L i and ␦e ix . Likewise, two kinds of variations are responsible for the variations in its orientation in this configuration: ␦l i , the parallelism error of the small sides of the parallelograms, ␦ Aiy and ␦␥ iy . Moreover, the sensitivities of the pose ͑position and orientation͒ of the end-effector to these variations are a minimum in this configuration, except for ␦l i . On the contrary, Q 2 configuration, i.e., P ϵ Q 2 , is the most sensitive configuration of the manipulator to variations in its design parameters. Indeed, as depicted by Figs. 21 and22, variations in the pose of the end-effector depend on all the design parameters variations and are a maximum in this configuration.

Moreover, Figs. 21 and 22 can be used to compute the variations in the position and the orientation of the end-effector with knowledge of the amount of variations in design parameters. For instance, let us assume that the parallelism error of the small sides of the parallelograms ͑␦l i ͒, is equal to 10 m. According to Fig. Let us assume now that the length and angular tolerances are equal to 0.05 mm and 0.03 deg, respectively. Figure 23͑a͒ shows the maximum position error of the end-effector when it follows diagonal Q 1 Q 2 of cube C u . Likewise, Fig. 23͑b͒ shows the maximum orientation error of the end effector along Q 1 Q 2 . On both sides, the error is a minimum when the manipulator is in its kinematic isotropic configuration and is a maximum in Q 2 configuration. Besides, the maximum position and orientation errors of the end-effector are equal to 0.7 mm and 0.4 deg, respectively. These values correspond to the worst case scenario, which is unlikely.

In order to take into account realistic machining tolerances, let us assume that the distribution of length and angular variations is normal.

Figures 24͑a͒ and24͑b͒ illustrate the repartition of the position and the orientation errors of the end-effector in Q 1 configuration, respectively.

Figures 25͑a͒ and25͑b͒ illustrate the repartition of the position and the orientation errors of the end-effector in the kinematic isotropic configuration, respectively.

Figures 26͑a͒ and26͑b͒ illustrate the repartition of the position and the orientation errors of the end-effector in Q 2 configuration, respectively.

In Figures 24͑a͒,25͑a͒, and 26͑a͒ ͑͑b͒, resp.͒, the horizontal axis depicts the Euclidean norm of the position ͑orientation, resp.͒ error of the end-effector and the vertical axis depicts the corresponding probability density function. To plot these figures, we computed the position and orientation errors of the end-effector corresponding to more than 3000 sets of geometric variations following a normal distribution. For example, we can deduce from these calculations the probabilities to get a position error lower than 0.3 mm and an orientation error lower than 0.25 deg in Q 1 , the isotropic, and Q 2 configurations.

According to Table 1, the probability to get a position error lower than 0.3 mm is higher in the kinematic isotropic configuration than in Q 1 and Q 2 configurations. However, the probability to get an orientation error lower than 0.25 deg is the same in Q 1 and the isotropic configurations, but is lower in Q 2 configuration.

Conclusions

In this paper, two complementary methods were introduced to analyze the sensitivity of a three-degree-of-freedom translational parallel kinematic machine with orthogonal linear joints: the Orthoglide. Although these methods were applied to a particular PKM, they can be readily applied to three-DOF Delta-Linear PKM such as ones with their linear joints parallel instead of orthogonal. Indeed, the input-output equations can be set in a very similar way since all Delta-linear PKM have identical leg kinematics, the only difference being in the closure equations ͓14͔.

On the one hand, a linkage kinematic analysis method was proposed to have a rough idea of the influence of the length variations of the manipulator on the location of its end-effector. On the other hand, a differential vector method was developed to study the influence of the length and angular variations in the parts of the manipulator on the position and orientation of its end-effector. This method has the advantage of taking into account the variations in the parallelograms.

According to the first method, variations in design parameters of the same type from one leg to another have the same effect on the end-effector. Besides, the position of the end-effector is very sensitive to variations in the lengths of parallelograms and prismatic joints. The second method shows that the parallelism errors of the bars of parallelograms are little influential on the position of the end-effector. Nevertheless, the orientation of the end-effector of the manipulator is more sensitive to the parallelism errors of the small sides of the parallelograms than to the ones of their long sides. Furthermore, it turns out that the sensitivity of the pose of the end-effector of the manipulator to geometric variations is a minimum in its kinematic isotropic configuration. On the contrary, this sensitivity approaches its maximum close to the kinematic singular configurations of the manipulator.

Therefore, these results should help the designer of the Orthoglide to synthesize its dimensional tolerances. Likewise, these methods can be applied to other Delta-Linear PKM with an aim of tolerance synthesis. Finally, the next steps in our research work are to compare the sensitivity of Delta-Linear PKM to variations in their geometric parameters, and to study the relation between the sensitivity and the stiffness of such manipulators.

Fig. 1 Fig

 1 Fig. 1 Basic kinematic architecture of the Orthoglide

Fig. 4

 4 Fig. 4 Mean of sensitivity of p x throughout C u Fig. 5 Mean of sensitivity of p y throughout C u

  depicts the global sensitivities of p, p x , p y , and p z to the dimensional variations. It appears that they are a minimum when P ϵ O, and a maximum when P ϵ Q 2 . Figures 12 and 13 depict the sensitivity coefficients of p x and p in the isotropic configuration. In this configuration, the position error of the end-effector does not depend on the orientation errors of the directions of the prismatic joints because the sensitivity of the position of P to variations in b 1y , b 1z , b 2x , b 2z , b 3x , b 3y is null in this configuration. Besides, variations in p x , p y , and p z are decoupled in this configuration. Indeed, variations in p x , ͑p y , p z , respectively͒ are only due to dimensional variations in the 1st, ͑2nd, 3rd, respectively͒ leg of the manipulator. The corresponding sensitivity coefficients are equal to 1. It means that the dimensional variations are neither amplified nor compensated in the isotropic configuration.

Fig. 7

 7 Fig. 7 Mean of sensitivity of p throughout C u Fig. 9 Sensitivity of p y to the variations in the 1st leg

  Fig. 12 Sensitivity of p x in the isotropic configuration

  the sum of the orientation errors of the ith parallelogram with respect to the ith prismatic joint and the end-effector; ␦m i = ␦b i -␦c i corresponds to the parallelism error of links B i1 C i1 and B i2 C i2 , which is depicted by Fig. 20.

Fig. 18

 18 Fig. 18 Variations in B i -B ij -C ij chain

Fig. 21

 21 Fig. 21 Sensitivity of the position of the end-effector along Q 1 Q 2 : "a… to dimensional variations and "b… to angular variations

Fig. 23

 23 Fig. 23 Maximum position "a… and orientation "b… errors of the end-effector along Q 1 Q 2

NomenclatureR

  b ͑O , x , y , z͒ ϭ reference coordinate frame centered at O, the intersection between the directions of the three actuated prismatic joints R P ͑P , X , Y , Z͒ ϭ coordinate frame attached to the end-effector R i ͑A i , x i , y i , z i ͒ ϭ coordinate frame attached to the ith prismatic joint, i =1,2,3 p = ͓p x p y p z ͔ T ϭ vector of the Cartesian coordinates of the end-effector, ex-pressed in R b ␦p = ͓␦p x ␦p y ␦p z ͔ T ϭ position error of the end-effector, expressed in R b ␦ = ͓␦ x ␦ y ␦ z ͔ T ϭ orientation error of the end- effector, expressed in R b i ϭ displacement of the ith prismatic joint ␦ i ϭ displacement error of the ith pris- matic joint L i ϭ theoretical length of the ith parallelogram A i , B i , C i ϭ depicted inFig. 1 a i ϭ distance between points O and A i r i ϭ distance between points P and C i b 1y , b 1z ϭ position errors of point B 1 along y and z axes, respectively b 2x , b 2z ϭ position errors of point B 2 along x and z axes, respectively b 3x , b 3y ϭ position errors of point B 3 along x and y axes, respectively

Fig. 26

 26 Fig. 26 Repartition of the position error "a… and the orientation error "b… of the end-effector in Q 2 configuration

h 1 ϭ b 1y k 1

 1 ϭ b 1z h 2 ϭ b 2x k 2 ϭ b 2x h 3 ϭ b 3x k 3 ϭ b 3y d i ϭ nominal width of the ith parallelogram ␦L i ϭ variation in the length of the ith parallelogram ␦L ij ϭ variation in the length of link B ij C ij , j =1,2 ͑see Fig. 2͒␦b i ϭ variation in the length of linkB i1 B i2␦c i ϭ variation in the length of linkC i1 C i2 ␦l i ϭ parallelism error of links B i1 B i2andC i1 C i2 ␦m i ϭ parallelism error of links B i1 C i1 and B i2 C i2 w i ϭ direction of links B i1 C i1 and B i2 C i2␦w i ϭ variation in the direction of linksB i1 C i1 and B i2 C i2 ␦e i ϭ sum of the position errors of points A i , B i , C i ␦ Ai = ͓␦ Aix ␦ Aiy ␦ Aiz ͔ T ϭ angular variation in the direction of the ith prismatic joint ␦ Bi = ͓␦ Bix ␦ Biy ␦ Biz ͔ T ϭ angular variation between B i1 B i2and the direction of the ith prismatic joint␦ Ci = ͓␦ Cix ␦ Ciy ␦ Ciz ͔ T ϭangular variation between the end-effector and C i1 C i2 ␦␥ i = ͓␦␥ ix ␦␥ iy ␦␥ iz ͔ T ϭ sum of the orientation errors of the ith parallelogram with respect to the ith prismatic joint and the end-effector DOF ϭ degree of freedom PKM ϭ parallel kinematic machine

  yr ͑C zr S yr S xr -S zr C xr ͒ ͑C zr S yr C xr + S zr S xr ͒ S zr C yr ͑S zr S yr S xr + C zr C xr ͒ ͑S zr S yr C xr -C zr S xr ͒

	-S yr	C yr S xr	C yr C xr

΅ ͑49͒ such that C xr = cos xr , S xr = sin xr , C yr = cos yr , S yr = sin yr , C zr = cos zr , S zr = sin zr , and

Table 1 Probabilities to get a position error lower than 0.3 mm and an orientation error lower than 0.25 deg in Q 1 , the isotro- pic, and Q 2 configurations

 1

			Configuration
		Q 1	Isotropic	Q 2
	Prob͑ʈ␦ p ʈ 2 ഛ 0.3 mm͒ Prob͑␦ ഛ 0.25 deg͒	0.8468 0.9691	0.9683 0.9690	0.7276 0.9453

Acknowledgment

This research, which is part of ROBEA-MP2 project, was supported by the CNRS ͑National Center of Scientific Research in France͒.


Matthieu Sacher
email: matthieu.sacher@ecole-navale.fr

Mathieu Durand

Elisa Berrini

Frédéric Hauville

Régis Duvigneau

Olivier Le Maitre

Jacques-André Astolfi

Élisa Berrini

Olivier Le Maître

Flexible hydrofoil optimization for the 35th America's Cup with constrained EGO method

Keywords: Constrained optimization, Gaussian process model, hydrofoil, cavitation, fluid-structure interaction

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

Introduction

Created in 1851, the America's Cup trophy is one of the oldest international sports competition. The race opposes two sailing boats, representing a defender and a challenger yacht clubs, that race one against the other to be the winner of the America's Cup match races. In each edition, specific rules are provided by the defender in accordance with the Deed of Gift, which define the format of the regatta, the location of the races and restrictions on the type of boats allowed to compete. When a challenging yacht club wins the match races, it becomes the defender and gains the stewardship for the next America's Cup edition. The 35 th edition, which took place in June 2017, has been raced on so-called AC50s catamarans with size fixed to 50 feet, equipped with hydrofoils and wing-sails. The hydrofoils are expected to be a key component of the AC50s performance, and the competing teams are dedicating a growing fraction of their design effort to their optimization in comparison to older yacht design approach [START_REF] Richards | America's Cup downwind sails-vertical wings or horizontal parachutes?[END_REF]]. The optimization aims at designing hydrofoil with a minimal drag force over a whole range of boat speeds (BS), while ensuring a sufficient lift force for the boat to fly, especially at the lowest BS. The design of the hydrofoil is restricted by the "AC Class Rule" documents. In particular, the rule number 11.11 states: "Daggerboard components shall be rigidly fixed to each other and the daggerboard shape shall not be adjusted while racing." One possible way to comply with this rule, while allowing for some adaptivity of the shape with the BS, is to design a hydrofoil composed of rigidly mounted components, undergoing large elastic deformations under the hydrodynamical loads and without crew adjustments [START_REF] Hueber | A Passively Morphing Trailing Edge Concept for Sailing Hydrofoil[END_REF]. For instance, a reduction of the hydrofoil camber, with increasing BS and hydrodynamic loads, can be obtained using a flexible trailing edge.

Drag reduction by shape optimization is an important research topic, with many applications in domains ranging from aeronautic to wind-energy recovery. Multiple studies have been proposed, based on physical models with different levels of complexity. For instance, the geometrical optimization of an airfoil with a two-dimensional steady flow model is considered in [START_REF] Ribeiro | An airfoil optimization technique for wind turbines[END_REF] when the unsteady flow situation is treated in [START_REF] Srinath | Optimal aerodynamic design of airfoils in unsteady viscous flows[END_REF]. Complex optimizations for a range of Mach numbers of a two-dimensional airfoil have been performed by [START_REF] Li | Robust airfoil optimization to achieve drag reduction over a range of Mach numbers[END_REF], and the robust optimization for uncertain flows and geometries is considered in Papadimitriou and [START_REF] Papadimitriou | Aerodynamic shape optimization for minimum robust drag and lift reliability constraint[END_REF]. The numerical optimization of hydrofoils is to our knowledge much rarer, and the hydrofoil literature generally focuses on the stability [START_REF] Inukai | Development of a single-handed hydrofoil sailing catamaran[END_REF], wake energy reduction [START_REF] Kandasamy | Multi-fidelity optimization of a high-speed foil-assisted semi-planing catamaran for low wake[END_REF], and free surface elevation [START_REF] Duvigneau | Hydrodynamic design using a derivative-free method[END_REF] aspects.

This paper reports the development and application, in partnership with the Groupama Team France (GTF) design team, of an Efficient Global Optimization (EGO) procedure tailored to the flexible hydrofoil optimization problem. We restrict ourselves to the optimization of a two-dimensional hydrofoil consisting of rigid forward and trailing parts, connected by a deformable element. The objective of the optimization is the reduction of the hydrodynamic drag at several BS and lift conditions, with cavitation constraints. Cavitation is indeed an important concern for the hydrofoil performance [START_REF] Sedlar | Numerical and experimental investigation of three-dimensional cavitating flow around the straight {NACA2412} hydrofoil[END_REF], but its numerical prediction remains a difficult problem, as shown for instance in [START_REF] Leroux | A joint experimental and numerical study of mechanisms associated to instability of partial cavitation on two-dimensional hydrofoil[END_REF], Coutier-Delgosha et al. [2007], [START_REF] Ducoin | Numerical and experimental investigation of hydrodynamic characteristics of deformable hydrofoils[END_REF], [START_REF] Akcabay | Cavity induced vibration of flexible hydrofoils[END_REF]. These difficulties explain that cavitation aspects are usually not considered in hydrofoil shape optimization, unless when the objective is precisely to delay the cavitation, such as in [START_REF] Wei | Numerical research on the performances of slot hydrofoil[END_REF].

The optimization of flexible hydrofoils raises many challenges. First it requires multi-design variables to prescribe the foil geometry at rest (unloaded shape), on the one hand, and the elastic properties and the internal structure of the deformable parts, on the other hand. Second, the evaluation of the performance of a given design requires the resolution of a fully nonlinear fluid-structure interaction (FSI) problem, possibly combining different boat speeds, with a numerical cost that inherently limits the number of designs that can be evaluated. Finally, for a given design, the consideration of non-cavitating flow requires an appropriate treatment of constraints that are costly to evaluate. As a matter of fact, the numerical optimization requires dedicated techniques to circumvent the computational complexity of the problem. In the present work, the numerical modeling involves the coupling of the two-dimensional vortex lattice method solver XFOIL [START_REF] Drela | XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils[END_REF]], to compute the flow, with the structural ARA software developed by K-Epsilon [START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF], to compute the hydrofoil deformations. Even for this two-dimensional modeling, the computational times prevent the direct optimization based on the FSI solution with standard methods.

We then rely on a surrogate model to reduce the computational burden of optimizing complex systems with costly objective function estimation [START_REF] Simpson | Metamodels for Computer-based Engineering Design: Survey and recommendations[END_REF]. Specifically, we use Gaussian processes [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF]] surrogate models with Efficient Global Optimization (EGO) strategies [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]], that have been previously applied to aerodynamic drag reduction [START_REF] Jeong | Efficient Optimization Design Method Using Kriging Model[END_REF], vibration reduction for rotating aircrafts [START_REF] Glaz | Helicopter Vibration Reduction throughout the Entire Flight Envelope Using Surrogate-Based Optimization[END_REF], optimization of FSI problems [Aghajari and Schäfer 2015] and sail trimming optimization [START_REF] Sacher | Efficient optimization procedure in non-linear fluid-structure interaction problem: Application to mainsail trimming in upwind conditions[END_REF]]. The classification approach for discontinuous or binary constraints proposed in [START_REF] Basudhar | Constrained efficient global optimization with support vector machines[END_REF] is extended to the treatment of the cavitation constraints in the optimization procedure. In our approach, The admissible domain is approximated using Least-Squares Support Vector Machine (LS-SVM) [START_REF] Suykens | Least Squares Support Vector Machine Classifiers[END_REF] regression from the previous observations. An original contribution of the work is a new probabilistic treatment of the constraints which is designed to mitigate the risk of exploring unfeasible solutions (which yield a cavitating flow).

This optimization procedure used in the present work is briefly described in Section 2, while Section 3 discusses the formulation of the optimization problem for the flexible hydrofoil, including the definition of the optimization variables and a brief discussion of the numerical models. Section 4 presents the optimization results, emphasizing on impact of the flexibility of the performance of the optimal hydrofoil. In particular, these performances are contrasted with the case of a rigid hydrofoil optimized for the same conditions. We have tried to provide as much as possible quantitative results, but for confidentiality reasons related to the partnership with the GTF design team, some values can not be provided and only trends are reported. Finally, conclusions and recommendations for future development are given in Section 5.

GP-based constrained optimization

The objective is to estimate the solution x opt of the following abstract optimization problem

min x∈Ω P(x), s.t. Q(x) ≥ 0,
where x ∈ R d is the vector of control variables, Ω the optimization domain, P : Ω → R the objective function and Q : x → R m is the vector of constraints. Even in the unconstrained case, m = 0, finding the global optimum of P can be very costly, in particular when its evaluation is numerically expensive. The use of surrogate models for P is then a classical solution to reduce the computational burden of complex systems optimization [START_REF] Simpson | Metamodels for Computer-based Engineering Design: Survey and recommendations[END_REF]]. We use Gaussian processes (GP) [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF]] which, owing to their statistical nature, provide both a prediction of the objective function and a measure of the uncertainty (variance) in the prediction. These features are appealing in optimization, as they can be used to derive rigorous optimization strategies based on the maximization of the Expected Improvement (EI) criterion, leading to methods referred globally as the Efficient Global Optimization (EGO) [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]] methods. In Section 2.1 we review the construction of the GP model for P and the resulting optimization strategy in the unconstrained case. The constrained case is considered in Section 2.2.

Unconstrained EGO using GP

Consider a set of n training inputs points X n = {x 1 , . . . ,x n }, each associated to a noisy observation y i of the objective function. Specifically, it is assumed that y i = P(x i ) + ε i , where the ε i are Gaussian measurement noises, assumed for simplicity independent and identically distributed with variance σ 2 . The GP construction considers that P is a realization of a zero-mean multivariate Gaussian process with covariance function C P . In this work, we consider the multidimensional squared exponential covariance functions defined by

C P (x,x ; Θ) . = θ 1 d i=1 exp       -(x i -x i ) 2 2l 2 i       + θ 2 , (1)
where Θ = {θ 1 ,θ 2 ,l 1 ,l 2 , . . . ,l d } is a vector of hyper-parameters to be inferred from the observations. From the conditional rules of joint Gaussian distributions [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], the best prediction P n (x) of P(x), i.e. the mean of y, and the prediction variance σ2 P (x) are given by

P n (x) = k T (x) C(Θ) + σ 2 I -1 Y n , (2)
σ2 P (x) = κ(x) + σ 2 -k T (x) C(Θ) + σ 2 I -1 k(x). ( 3
)
where C(Θ) ∈ R n×n is the covariance matrix with entries C i, j (Θ) . = C P (x i ,x j ; Θ), 1 ≤ i, j ≤ n, κ(x) . = C P (x,x; Θ), k(x) is the covariance vector between the observations in X n and x, and I is the identity matrix of R n . The hyper-parameters Θ and noise variance σ 2 can be determined by maximizing the logmarginal likelihood [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]] using an evolution strategy algorithm [START_REF] Hansen | The CMA Evolution Strategy: A Comparing Review[END_REF]]. More details on GP meta-models can be found in [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

Let xn be the optimum of P n (x). It is expected that xn ≈ x opt if the approximation error P n -P is small enough. The advantage of minimizing P n instead of P is that GP models are usually inexpensive to evaluate compared to the original objective function. To control the error in the approximation, one proceeds sequentially by adding progressively new observation points in the area of interest. A deterministic optimization procedure would choose the next point x n+1 as the optimal point of P n . However, the GP model provides probabilistic information that can be exploited to propose more robust strategies based on merit functions, which combine the prediction and its variance. In this work, we use the Augmented Expected Improvement (AEI) merit function [START_REF] Huang | Global Optimization of Stochastic Black-Box Systems via Sequential Kriging Meta-Models[END_REF]], which estimates the expected progress in the objective, taking into account the noise in the observed values and the prediction variance:

AEI(x) = EI(x)             1 - σ σ2 P (x) + σ 2             , (4)
where the Expected Improvement EI(x) is defined by

EI(x) = σP (x) u(x)Φ (u(x)) + φ (u(x)) , (5)
u(x) = P n (x * ,n ) -P n (x) σP (x) , (6)
with Φ and φ the cumulative and density functions of the standard Gaussian distribution, and x * ,n ∈ X n is the current effective best solution (see [START_REF] Huang | Global Optimization of Stochastic Black-Box Systems via Sequential Kriging Meta-Models[END_REF]). The optimum x n+1 of the AEI is added to X n , and P is evaluated at the new point providing y n+1 . Setting n ← n+1, a new iteration can start updating the GP model with the new observation. Overall, each iteration requires one computation of the objective and the resolution of two optimization problems: a first one for the hyper-parameters of the GP model, and a second one to find the AEI optimum.

The iterations of the GP-based optimization problem are continued until a stopping criterion is satisfied or the resources allocated to the optimization have been exhausted.

2.2. EGO under constraints, a classification based approach EGO methods with inequality constraints was considered in [START_REF] Schonlau | Computer Experiments and Global Optimization[END_REF]. The key idea was to rely on m additional surrogates to estimate the constraints Q i (x). For Gaussian Process models, one can easily determine the probability P i (Q i ≥ 0|x) that a constraint Q i is satisfied at x and, assuming the independence of the constraints probability, to come up with the consolidated probability

P(Q ≥ 0|x) = Π m i=1 P i (Q i ≥ 0|x
). This probability is used to modify the unconstrained AEI criterion and favor feasible regions:

AEI Q (x) = AEI(x)P(Q ≥ 0|x).
Although effective in many problems, this GP modeling of the constraints faces several limitations. First, its computational cost increases with the number m of constraints and can be an issue for problems with large m. Second, the approximation by GP models assumes a sufficient smoothness of the Q i , which must be evaluable for almost every x ∈ Ω. This rules out the case of binary constraints (feasible / infeasible) and models having no solution for Q(x) 0. Finally, the approach is also limited to situations where the feasible domain has an explicit representation in terms of the constraints. The optimization of hydrofoils reaches some of these limitations as it is difficult to express feasible geometries in terms of constraints on design parameters, while some values may lead to uncomputable solutions.

Classification methods recently proposed in [START_REF] Basudhar | Constrained efficient global optimization with support vector machines[END_REF] are better suited to deal with discontinuous or binary constraints in a GP-based optimization procedure. The binary classifier considers two classes C + and C -over Ω, corresponding to the feasible and unfeasible domains respectively. Each x i of X n is equipped with a value z i = ±1 depending on its membership C ± . To predict the class of a new point x we introduce a classification function h : x ∈ Ω → R, such that z(x) = signh(x). A Least-Squares Support Vector Machine (LS-SVM) [START_REF] Suykens | Least Squares Support Vector Machine Classifiers[END_REF] is used to construct h. The LS-SVM method extends the original Support Vector Machine (SVM) classifier [START_REF] Vapnik | The nature of statistical learning theory[END_REF] to quadratic penalization, resulting in a linear (but non-sparse) system to be solved (see below). Further details and discussion on SVM and LS-SVM methods can be found in [START_REF] Cawley | Leave-one-out cross-validation based model selection criteria for weighted LS-SVMs[END_REF].

The LS-SVM method [START_REF] Suykens | Least Squares Support Vector Machine Classifiers[END_REF]] is a linear classifier, for C + and C -, in a feature space induced by the transformation φ : Ω → F :

h(x) = w T φ(x) + b. (7)
Here, w and φ(x) are the weights and features vectors, while b ∈ R is a constant. The feature space F is generated by a reproducing kernel r :

Ω × Ω → R representing the inner product between images in F of vectors: r(x,x ) = φ(x) • φ(x ). We use the classical Gaussian kernel, r(x,x ) = exp - x -x 2 2λ 2 , (8)
with scale factor λ ∈ R to be adjusted. The LS-SVM parameters (w,b) satisfy the primal constrained optimization problem min w,b,e

1 2 w 2 + γ 1 2 n i=1 e 2 i , s.t. z i = w T φ(x i ) + b + e i , i = 1, . . . ,n, (9)
with trade-off parameter γ ∈ R + (to be fixed) and relaxation variables e i allowing for miss-classification. By writing the Lagrangian of the optimal problem, we introduce the α i ∈ R, that are the Lagrange multipliers of the constraints. Denoting

z = (z 1 • • • z n ) T
, the optimality conditions of the Lagrangian are used to derive a linear system for the dual model parameters b

and α = (α 1 • • • α n ) T ,         R + γ -1 I 1 1 T 0                 α b         =         z 0         , (10)
where R ∈ R n×n is the kernel matrix with R i, j = r(x i ,x j ). This system is solved via a Cholesky factorization [START_REF] Cawley | Leave-one-out cross-validation based model selection criteria for weighted LS-SVMs[END_REF]], and h in ( 7) is expressed in terms of dual model parameters to obtain

h (x) = n i=1 α i r(x i ,x) + b. (11)
The LS-SVM classifier depends on two parameters, γ and λ, to be fixed when solving (10). We rely on a Leave-One-Out procedure to determine γ and λ that minimize the Predicted Residual Sum-of-Squares (PRESS) criterion [START_REF] Allen | The relationship between variable selection and data augmentation and a method for prediction[END_REF]]. A simple expression for the predicted residuals is proposed in [START_REF] Cawley | Efficient leave-one-out cross-validation of kernel fisher discriminant classifiers[END_REF], which is presently minimized for λ and γ with the CMA-ES algorithm [START_REF] Hansen | The CMA Evolution Strategy: A Comparing Review[END_REF]].

Finally, the LS-SVM binary classifier is extended to a probabilistic classification, relating h in (11) to the probability of the class C + , denoted P(C + |x). A comparison of several probability models for the LS-SVM classification is provided in [START_REF] Van Calster | Proceedings, Part II, chap. Comparing Methods for Multi-class Probabilities in Medical Decision Making Using LS-SVMs and Kernel Logistic Regression[END_REF]. We use the sigmoid function [START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF]] and expressing the probability of C + as

P C + |x = 1 1 + exp (Ah (x) + B) . ( 12
)
The parameters A and B of the sigmoid are determined by minimizing the probability of misclassification, see [START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF][START_REF] Lin | A note on Platt's probabilistic outputs for support vector machines[END_REF]]. In practice, the probability P goes to 1 (resp. 0) as the classifier is certain that x is feasible and belongs to C + (resp. unfeasible and belongs to C -), while a value of P = 1/2 denotes a complete uncertainty in the classification. This can occur because x is far from any observations in X n or close to the interface between the two classes.

Following [START_REF] Schonlau | Computer Experiments and Global Optimization[END_REF]] and [START_REF] Basudhar | Constrained efficient global optimization with support vector machines[END_REF], the probability of feasibility / evaluability is used to derive from the AEI a merit function and select a new point x n+1 that yields the highest expected improvement while having a high probability of feasibility / evaluability. To this end, [START_REF] Basudhar | Constrained efficient global optimization with support vector machines[END_REF] uses a modified version of the sigmoid function in (12) that accounts for the distances to the closest classified observations (in X n ). While improving the misclassification rate, this modification results in a discontinuous merit function which prevents the use of gradient-based optimization tools to determine x n+1 . Further, their strategy requires the addition of d + 1 new observation points per iteration (to maintain isotropy in the observation points distribution) with an associated cost deemed too important in our application. Instead, we use the original definition (12) to extend the unconstrained AEI and we consider a sequential infilling strategy that follows the EGO iterations, by selecting alternatively one of the present definitions for x n+1 :

x n+1 = arg max x∈Ω AEI(x)P(C + |x), (13)
x n+1 = arg max x∈Ω AEI(x) s.t. P(C + |x) ≥ ρ, (14)
x n+1 = arg min x∈Ω P n (x) s.t. P(C + |x) ≥ ρ, (15)
x n+1 = arg max x∈Ω AEI(x)P(C + |x) 1 -P(C + |x) . ( 16
)
The first definition in (13) corresponds to the extension of the AEI favoring points with high chance of feasibility. The second expression ( 14) maximizes the original AEI definition by enforcing a minimal probability ρ of feasibility; we use ρ = 0.5. In ( 15), the predicted cost is directly minimized and enforced by ρ. Finally (16) combines the feasibility and infeasibility probabilities to favor areas where the classification is the most uncertain (P ∼ 0.5), with the objective to improve the exploration of the feasible domain boundaries. In practice, the new points x n+1 are determined using CMA-ES algorithms without [START_REF] Hansen | The CMA Evolution Strategy: A Comparing Review[END_REF] or with constraints Arnold and Hansen [2012], depending on the considered definition. As a final note, we remark that the optimization procedure can generate points that are found unfeasible in the sense that P(x n+1 ) cannot be computed. This is typically the case when x n+1 corresponds to an impossible geometrical configuration, or to an extreme situation for which the numerical code is not able to converge. It that case, the point is deemed infeasible, setting z n+1 = -1, but no value of the objective function is provided and the construction of the GP model for P simply disregard the missing data. As a result, the GP model of P and the LS-SVM classification can involve different numbers of observations.

Hydrofoil optimization setup

This Section introduces the formulation of the constrained optimization problem for the flexible hydrofoil, as defined in collaboration with GTF. Note that some details of the optimization problem are omitted for confidentiality reasons. We also briefly discuss the fluid and structural solvers used in this work.

Problem formulation

The optimization of the hydrofoil uses control variables x which parameterize the geometry and of the hydrofoil at rest, i.e. in the absence of fluid loads, and the elastic properties of the element that gives flexibility to the trailing edge (see below). The optimization concerns the minimization of the hydrofoil drag force over a set of 4 sailing conditions, with an additional set of 4 constraints designed to prevent the cavitation of the flow. Each condition corresponds to a prescribed boat speed and an hydrodynamic lift force developed by the hydrofoil. Note that the 4 drag points and the 4 constraints correspond to different conditions so the evaluation of the objective function P and the vector of constraints Q at given x requires solving 8 fluid-structure interaction problems. The optimization problem is written as

min x∈Ω P(x) = 4 i=1 w i C d i (x), s.t.                  Q 1 (x) = -Cp min 1 (x) ≤ λ 1 Q 2 (x) = -Cp min 2 (x) ≤ λ 2 Q 3 (x) = -Cp min 3 (x) ≤ λ 3 Q 4 (x) = -Cp min 4 (x) ≤ λ 4 . (17)
Here, the C d i are the drag force coefficients, the w i > 0 are prescribed weights, the Cp min i are the minimum pressure coefficients, and the λ i are the cavitation numbers defined by

λ = p 0 -p v 1 2 ρU 2 0 . (18)
We have denoted U 0 and p 0 the reference velocity (that is, the boat-speed) and pressure, ρ the fluid density, and p v the saturated vapor pressure. The reference pressure is taken as p 0 = p atm + ρgh, where p atm is the atmospheric pressure, g is the gravity acceleration, and h is the immersion depth of the hydrofoil. The minimum of the pressure coefficient is defined by Cp min = p min -p 0 0.5ρU 2 0

, where p min is the minimal pressure over the hydrofoil surface. The simplest criterion to prevent cavitation is to ensure that the minimal pressure remains higher than the vapor pressure: p min ≥ p v . This condition can be expressed in terms of cavitation number λ to obtain the conditions on the minimum pressure coefficient Cp min to formulate the constraints in (17).

Hydrofoil parametrization

The geometric and elastic characteristics of the flexible hydrofoil are defined starting from a baseline geometry provided by GTF. The baseline geometry was produced by an optimization procedure using an evolutionary method [START_REF] Bäck | An Overview of Evolutionary Algorithms for Parameter Optimization[END_REF] and assuming a complete rigidity of the hydrofoil (no elastic effects). The geometry at rest of the flexible hydrofoil is defined by rotating the sections of the baseline geometry (rigid hydrofoil), without changing the sections thickness. We denote p(s) the law of the sections rotation angle, where 0 ≤ s ≤ 1 is the reduced coordinate along the chord of the baseline geometry. In practice we restrict the rotation to the last 30% of the hydrofoil chord (p(s) = 0, s ≤ 0.7) and use a smooth B-Spline interpolation [START_REF] Piegl | The NURBS Book (2Nd Ed[END_REF]Tiller 1997, Wang et al. 2006] of p(s) between 4 control points 0.7 < s 1 < • • • < s 4 = 1, as illustrated in Figure 1. We shall refer to the rotation angles p 1,...,4 at the control points as the geometrical parameters. Regarding the flexibility of the hydrofoil, the results presented in the following consider a perfectly rigid hydrofoil except for an elastic element located in the backward part of the hydrofoil. As illustrated in Figure 2, the location of the elastic element is prescribed by the reduced coordinates of its starting and ending points 0 < s 1 < s 2 < 1 (resp. 0 < ds 1 < ds 2 < 1) on the lower (resp. upper) surface of the baseline hydrofoil. The deformable element is made of a homogeneous material having a Young modulus E LST ; it is further reinforced by a beam (shown in black in Figure 2) on the pressure surface (lower side). The elastic beam has a progressive Young modulus varying linearly between s 1 and s 2 , with mean value E Beam and slope C Beam . The use of a single, homogeneous, flexible element is supported by preliminary optimization studies briefly presented in the Appendix at the end of the paper. Overall, the parametrization of the deformable part involves 7 parameters (4 geometric coordinates and 3 elastic properties) globally referred to as the structural parameters.

Fluid-structure interaction solver

Given a value of the geometric and structural parameters, the evaluation of the objective function and the constraints of the optimization problem calls for the resolution of 4 FSI problems (one for each condition of the objective function). We assume that these problems have steady solutions for all values of the parameters and we rely on a quasi-steady approach to treat the coupling between the structural and fluid solvers.

For the flow solver, we rely on XFOIL [START_REF] Drela | XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils[END_REF]] which has been developed for many years and is widely used for many aerodynamic applications. While based on a two-dimensional airfoil panel code method, XFOIL can provide viscous and inviscid analyses. It incorporates a two-equation integral formulation of the viscous boundary layer and the approximate e N envelope method [START_REF] Van Ingen | The eN method for transition prediction: historical review of work at TU Delft[END_REF]; it allows the prediction of transition points and separation bubble, through the resolution of the boundary layer and the transition equations which is based on a Newton method. Further details on XFOIL are given in [START_REF] Drela | XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils[END_REF] and it is compared with other CFD methods and solvers in [START_REF] Morgado | {XFOIL} vs {CFD} performance predictions for high lift low Reynolds number airfoils[END_REF]. We only mention that the setting of our optimization problem imposes to enforce a prescribed lift force in the flow computation. This is achieved in XFOIL by determining the Angle of Attack (AoA) that yields the requested lift force. Overall, XFOIL is simple to use and has a reasonable computation cost on a classical workstation, in most of the cases, making it a natural candidate to be coupled with a structural solver (see for instance [MacPhee and Beyene 2013] for a recent FSI study of wind turbine blades using XFOIL).

The ARA software developed by K-Epsilon was used to solve the structural model of the flexible hydrofoil. The solver was originally developed for FSI computations of sails [START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF]] and proposes various structural elements [START_REF] Durand | Interaction fluide-structure souple et légère, application aux voiliers[END_REF]]. In this work, the elastic part of the hydrofoil is modeled by 2D Linear Strain Triangles (LST) [START_REF] Pedersen | Some properties of linear strain triangles and optimal finite element models[END_REF]] and Timoshenko beam elements (see Figure 2). The nonlinear equilibrium solution is computed by a Newton method with Aitken relaxation.

A quasi-monolithic algorithm [START_REF] Durand | Interaction fluide-structure souple et légère, application aux voiliers[END_REF]] is finally applied to couple the structural and flow solvers. Briefly, in this algorithm, the resolution of the structural problem is nested inside the iterations of the nonlinear steady flow solver. This ap-proach preserves the convergence rate and stability properties of the monolithic approach. The resulting coupled solver is finally driven by a utility that computes the set of FSI solutions required to estimate the objective function value and its constraints (17). This utility is itself nested into the optimization driver that decides of the sequence of new optimal parameters x n+1 to be evaluated.

Results and discussions

We now present the results of the optimization problem. In particular, we contrast the optimal flexible hydrofoil with its rigid counterpart optimized for the same set of conditions and constraints. This comparison concerns the characteristics and performances of the two optimal hydrofoils. We do not provide a comparison in terms of computational times because the two optimizations use different strategies and have significantly different complexities. Indeed, the optimization of the rigid hydrofoil although involving a larger number of optimization variables has an objective function and constraints much faster to compute than in the flexible case (which requires the resolution of FSI problems), enabling the application of evolutionarybased methods [START_REF] Bäck | An Overview of Evolutionary Algorithms for Parameter Optimization[END_REF].

In the following, the rigid hydrofoil corresponds to the optimized baseline geometry provided by GTF, assuming a perfectly rigid hydrofoil, while the unloaded geometry refers to the rest geometry resulting from our optimization of the flexible hydrofoil, that is, the rigid geometry with the rotations of the trailing edge sections (but without elastic deformations). We stress that the structural and geometrical parameters characterizing the flexible hydrofoil and are not changing with the boat speed. Thus, the rigid, unloaded, and in fact, all flexible geometries are coinciding if the elastic element is infinitely stiff, because the rigid geometry is optimal in the absence of elastic deformations. For a finite elastic stiffness, one needs, in general, to jointly optimize the geometrical and structural parameters of the flexible hydrofoil to improve the performances compared to the rigid case. Optimizing only for the structural (resp. geometrical) parameters usually leads to a degradation of the performance, compared to the rigid hydrofoil case, unless the geometrical (resp. structural) parameters are fixed close to their joint optimization solution values.

Selection of the design variables

The computational time for the optimization of the flexible hydrofoil is an important aspect, and several tests were made to determine an appropriate trade-off between the computational and parametrization complexities, before arriving at the set-up described in Section 3.2 with the up to 11 parameters to be optimized (4 geometrical and 7 structural). The number of optimization variables can be further reduced, fixing, for instance, the material properties, the geometry of the deformable element, or changing the number of control points in the parametrization of the rotation rule. These reductions yield different complexities of the optimization problem with different computational times as a result. However, reducing the number of optimization variables also impacts the performance of the optimal design if the fixed parameters are not set close to their optimal values. Figure 3 shows the evolutions with the computational time of the estimated minimum of the objective function for different numbers of optimization variables considered in the optimization problem. The plot shows that when nine variables are considered the optimization needs ≈ 16 times more resources to converge than for only five optimization variables, with close optimal objective values in this case. The case of two optimization variables is even less computationally demanding, but leads to a noticeable drop in the optimal performance because of inappropriately set values for the parameters disregarded in the optimization problem. These effects can be better appreciated from on the case of 9 optimization variables, corresponding to the 4 geometrical and 5 of the 7 structural parameters. E Beam and C Beam have been fixed close to their optimization solution values.

Optimal shapes and pressure distributions

Figures 4 and5 show for 20 and 40 kts of boat-speed respectively, the shapes of the rigid and flexible hydrofoils. These shapes are depicted at their effective angle of attack giving the prescribed lift force. The pressure coefficient distributions are also reported in the top part of the Figures. At low boat-speed (20 kts), in Figure 4, the rigid and deformed hydrofoils have noticeable differences in term of angle of attack: the rigid one is higher. Regarding the pressure coefficient distributions, the rigid case presents a minimum at the leading edge of ≈ -2.25, whereas in the flexible case the minimum peaks to Cp ≈ -1.1 only, at the same location. This can be explained as a less cambered geometry generally requires a higher angle of attack to achieve the same lift force, while the minimum pressure coefficient usually increases with the angle of attack. Anyway, the results indicate a poorer tolerance to cavitation for the rigid hydrofoil in these conditions. For 40 kts of boat-speed, shown in Figure 5, the two angle of attack are negative and the rigid one presents again the largest deviation from zero angle of attack. The deformation of the trailing edge is now pronounced, in the flexible case, as it can be appreciated using dashed black contour which corresponds to the unloaded flexible geometry (which comes from the GTF optimizations with flap laws). The trailing edge deflection is ≈ 9.5 deg. The minimum pressure coefficients are now located approximately at mid-chord with Cp min ≈ -0.5 in the two cases. The distortions in the pressure coefficient distribution observed at ≈ 0.8 percent of the chord in the case of the flexible hydrofoil are due to the transitions between the rigid parts and the deformable element constituting the flexible hydrofoil.

The flow around the flexible hydrofoil at 40 kts of boatspeed has also been computed with a URANS transition model in order to validate the pressure coefficient distribution predicted by XFOIL. For that purpose, a structured mesh of the fluid domain around the hydrofoil has been generated using ≈ 110 000 elements; the flow was subsequently solved using FLUENT with the Transition SST turbulence model [START_REF] Menter | A correlationbased transition model using local variables-Part I: model formulation[END_REF]]. The computed pressure coefficient distribution is reported in Figure 6, together with the XFOIL prediction, highlighting the excellent agreement between XFOIL and the Transition SST flow model. Small differences in the pressure coefficients are visible only on the pressure side at ≈ 90% of the chord length, and the magnitude of the differences in the global hydrodynamic loads are 2.6% and 4.2% for the lift and drag coefficients respectively. Thus, the pressure coefficients computed by XFOIL can be considered reliable, even in the area of the deformable element. Figure 7 shows the evolutions of the camber along the chord for the rigid, unloaded and the loaded flexible hydrofoils at several boat-speeds. These geometrical quantities are computed based on the distances in the direction normal to the deformed chord, following the classical conventions. Note that because of the changes of the displacement of the trailing edge and the change in chord length with the boat-speed, for the flexible hydrofoil, the camber law in the forward part before the flexible element is also affected. The flexibility is seen to have an important effect for a boat-speed larger than 20 kts. In addition, the shapes of the camber distribution of the flexible hydrofoil differ significantly from the rigid one for all the boat-speeds shown. Specifically, the maximum camber of the flexible foil at 20 kts of boat-speed is more than twice that of the rigid one, while at 40 kts of boat-speed it is less by roughly 20%. These evolutions of the apparent camber distribution with the boatspeed (and load in general) is one of the hydrodynamical advantages brought by the flexibility.

Angles of attack and cavitation criterion

The effective angle of attack of the optimal hydrofoils is reported in Figure 8 as a function of the boat-speed. The previous observations regarding the angle of attack are confirmed. The angle of attack of the rigid hydrofoil almost linearly decreases from ≈ 3.9 deg at 20 kts to ≈ -0.9 deg at 30 kts of boat-speed, and subsequently continues to decrease, but at a lower rate, when the boat-speed increases further. For the flexible case, the maximum of angle of attack is only ≈ 1.8 deg at 20 kts and decreases to ≈ -0.1 deg at 25 kts of boat-speed. Then, unlike in the rigid case, the flexible hydrofoil keeps an almost constant effective angle of attack between 0 and -1 deg for the rest of boat-speed range, instead of decreasing continuously in the rigid case. Again, the flexibility can explain this behavior, as the deformation of the trailing edge (and reduction of the camber) discharges the foil when the boat-speed increases. This discharge calls for fewer changes in the angle of attack compared to the rigid case. Note that this sort of automatic trimming, through flexibility, may also be beneficial by saving some trimming efforts to the crew members. Another positive impact of the flexible hydrofoil can be appreciated in Figure 9 which reports the minimum pressure coefficients for the two optimal hydrofoils, as a function of the boat-speed. The dashed black curve corresponds to the critical value based on the cavitation number (18). We thus expect the optimal hydrofoils to have a minimum pressure coefficient greater than this critical value over the whole range of boatspeeds considered. However, we recall that the cavitation criterion is enforced for a finite set of only 4 distinct values of the boat-speed. For the maximal boat-speed shown, 40 kts, the two optimal foils, in fact, violate by a small margin the criteria based on the cavitation number, which indeed is not explicitly enforced for this boat-speed. As lower boat-speeds are considered, till 25 kts, the two foils satisfy the non-cavitating criterion by a large margin, suggesting that any constraints in this range of boat-speeds would not be active. Finally, when the boatspeed approaches 20 kts the behavior of the minimum pressure coefficient differs between the two hydrofoils: while for the flexible hydrofoil Cp min remains away from its critical value, it decreases sharply and even becomes significantly lower than its critical value in the rigid case. Comparing the trends of the two hydrofoils, it can be concluded that the flexibility yields a design with minimum pressure coefficient much less dependent on the boat-speed, and therefore having lower chances of violating the non-cavitating conditions for a boat-speed other than in the constraints.

BS [kts]

AoA [deg]

Drag performances

To complete the comparison between the flexible and rigid hydrofoil, we report in Figure 10 the drag coefficients as a function of the boat-speed. The actual values are not shown, for confidentiality reasons, but the plot allows to appreciate the improvement brought by the flexible hydrofoil. Except for the lowest boat-speed (20 kts), the flexible hydrofoil has a drag coefficient C d consistently lower than for the rigid one. The improvement is particularly significant at 25 and 40 kts of boatspeed, with roughly 20% drag reduction at the highest boatspeed. Again, these improvements come from the flexibility which allows the hydrofoil to naturally adapt its camber with the boat-speed, with a thinner geometry at rest and lower variations of the angle of attack in the boat-speed range.

For a boat-speed between 20 and 25 kts, the hydrofoil drag is mainly depending on the location of the transition point on the suction side, denoted x tr . As the boat-speed is lowered, the angle of attack increases to satisfy the required lift force and the transition point x tr moves toward the leading edge with a sharp increase in the friction drag (see Figure 10). Note that in XFOIL the location of the transition point is computed whenever it appears to be in the first 40% of the chord. Otherwise, it is set at this maximal location. To highlight the importance of the transition point on the drag of the hydrofoil, Figure 11 presents the evolutions of the flexible hydrofoil drag and the location of the transition point as functions of the boat-speed. The sharp drop in the hydrofoil drag, around 22.5 kts of boat-speed, is clearly related to the displacement of the transition point away from the leading edge. Note that by using different weights and adding more conditions in (17), one could further improve the hydrofoil drag at a low boat-speed, by delaying the displacement of the transition point, but to the detriment of the drag at higher boat-speeds.

Optimal hydrofoil at 30 kts

To illustrate the importance of the selected conditions used in the definition of the objective function, we optimize another rigid hydrofoil based on the minimization of the drag and noncavitating constraint in (17) based on the unique condition corresponding to 30 kts of boat-speed. The design variables are reduced to the 4 geometrical angles p 1,...,4 and we refer to this optimized hydrofoil as the Rigid 30 kts one. 30 kts can be further appreciated from Figure 12, where the two distributions of pressure coefficients Cp and the shapes plotted at the same angle of attack for comparison purposes are shown. The main difference in the pressure coefficient is the smoother character of the distribution in the Rigid 30 kts case, which can be explained by the absence of the localized elastic deformations present in the flexible case. However, the differences between the two distributions are small as one could have expected from the similarity of the shapes and effective angle of attack. However, as soon as the boat-speed departs from 30 kts, the flexible hydrofoil exhibits a lower drag compared to the Rigid 30 kts one. This can be appreciated in Figure 13 which reports the differences in the drag coefficients, with positive values in favor of the flexible case. Note that these consistent improvements were not necessarily expected as the Rigid 30 kts case does not enforce any non-cavitating conditions at boat-speeds other than 30 kts. Therefore, the Rigid 30 kts optimization may have produced a cavitating design with lower drag for other boat-speeds. Note also the higher drag for the flexible hydrofoil, at 20 kts of boat-speed, where the cavitation criterion is actually not satisfied by the rigid solution.

Conclusions

An advanced multi-design constrained optimization method has been proposed in this paper for the design of a two-dimensional hydrofoil with a flexible trailing edge. Surrogate-model based optimization, with inequality constraints, have been used to enable the optimization of the hydrofoil at a reasonable computational cost. The surrogates are constructed and updated using a limited number of carefully selected resolutions of the fully nonlinear fluid-structure interaction problem. The nonlinear FSI solver used in this work couples a two-dimensional vortex lattice method code for the flow around the hydrofoil with a static structural solver, based on linear stress triangle and Timoshenko beam elements, for the hydrofoil deformations. The optimization involves 4 parameters describing the rotation of a base geometry, and seven elastic parameters describing the properties of the elastic bound between the rigid leading part and trailing edge of the hydrofoil. These parameters are sought to globally minimize the hydrofoil drag forces in selected conditions (boat speeds and lift forces), while constraints were introduced to ensure noncavitating conditions.

The performances of the resulting optimal flexible hydrofoil have been compared to the rigid case optimized for the same set of conditions. It was found that the flexible hydrofoil performs globally better over the whole set of conditions, compared to the rigid one. In fact, the flexible hydrofoil was found to have drag forces lower than that of the rigid one the whole range of boat speeds, except around 20 knots where the two cases have essentially the same drag force. The improvement in the drag force brought by the flexibility is particularly noticeable for the highest boat speeds (≥ 35 knots) and around 25 knots. A closer inspection of the optimal solutions revealed that not only the flexible hydrofoil performs better for most conditions, but its angle of attack varies with the boat speed in a tighter range to achieve the prescribed lift conditions, compared to the rigid case. In particular, the flexible hydrofoil has a much smaller angle of attack at low boat speed, which immediately translates into a higher minimum of pressure and less susceptibility to cavitation.

Overall, the optimization clearly evidences the various advantages of considering a hydrofoil with a flexible trailing edge. In particular, it demonstrates that, if correctly designed, the flexible trailing edge allows for an auto-adaption of the hydrofoil camber line and angle of attack with the boat speed, with a global reduction of the drag and lower cavitation risk.

Future developments will concern the optimization of a threedimensional flexible hydrofoil with the use of more advanced flow models to account for the three-dimensional effects, turbulence, and more complex cavitation criteria. We are currently exploring multi-fidelity strategies [START_REF] Park | Remarks on multi-fidelity surrogates[END_REF] to tempered the increase in the computational costs in the three-dimensional case. Another aspect requiring further investigation is the dependence of the transition point displacement as a function of the hydrofoil angle of attack at lower boat speed range. Indeed, we believe that the flexible trailing edge could be effective in producing enough lift for the platform to take off at a lower boat speed.

Appendix: Preliminary study

We report in this appendix an optimization result to support the definition of the elastic hydrofoil used in the paper. Specifically, we considered a more general definition of the elastic hydrofoil using six consecutive homogeneous elastic elements in the elastic part ranging in the last 30 % of the hydrofoil chord (see Figure 14). The elastic elements have a fixed size and their respective Young moduli E i=1,...,6 are optimized for the same objective function and constraints in (17), together with 3 geometrical parameters (rotations p 1,...,3 defined in Section 3.2), leading to an optimization problem with 9 design variables. The Young moduli are also constrained to remain in a finite (positive) range of available materials. The optimal Young moduli E 1,...,6 are reported (as a function of the chord length coordinate), in Figure 15. The values have been divided by a reference value for confidentiality reasons. It is seen that the optimal solution corresponds a localized flexibility, located at ≈ 85 % of the chord length, and with the first and last elastic elements of the flexible hydrofoil having Young moduli equal to the maximum value allowed. In other words, the optimal design presents very stiff upstream part and trailing edge, separated by a significantly softer region. An idealized version of this solution is the case of a single homogeneous elastic element described in Section 3.2 and depicted in Figure 2.

Figure 1 :

 1 Figure 1: Illustration of the unloaded geometry definition using a rotation of the baseline geometry. The top plot shows the B-Spline approximation of the rotation rule p(s) based on the rotation angles p 1 to p 4 at the control points. The bottom plot compares the baseline (red) and unloaded flexible geometries, the green vertical line indicating the chord position s = 0.7.

Figure 2 :

 2 Figure 2: Illustration of the structural parameters.

Figure 3 :

 3 Figure 3: Convergence of the objective function minimum with the computational time, for different numbers of optimization variables considered.

Figure 4 :

 4 Figure 4: Shapes and distributions of Cp at 20 kts of boat-speed for the rigid and flexible hydrofoils.

Figure 5 :

 5 Figure 5: Shapes and distributions of the pressure coefficient Cp at 40 kts of boat-speed for the rigid and flexible hydrofoils.

Figure 6 :

 6 Figure 6: Distributions of the pressure coefficients Cp predicted by XFOIL and with the URANS transition models.

Figure 7 :

 7 Figure 7: Flexible hydrofoil camber distribution at several boat-speeds. The rigid hydrofoil case is also shown for comparison.

Figure 8 :

 8 Figure 8: Effective angle of attack as a function of the boat-speed (BS).

Figure 9 :

 9 Figure9: Minimum pressure coefficients Cp min as a function of the boat-speed (BS). Also shown is the cavitation number λ given by (18).

Figure 10 :Figure 11 :

 1011 Figure 10: Drag C d of the optimal hydrofoils as a function of the boat-speed (BS).

Figure 12 :

 12 Figure 12: Shapes and pressure coefficient Cp distributions for the flexible and Rigid 30 kts solutions at 30 kts of boat-speed.

Figure 13 :

 13 Figure 13: Difference ∆C d in the drags of the Rigid 30 kts and flexible optima, as a function of the boat-speed (BS).

Figure 14 :

 14 Figure 14: Flexible hydrofoil composed of 6 elastic elements.

Figure 15 :

 15 Figure 15: Optimal distribution of the Young modulus (arbitrary units) with the chord length.

Table 1

 1

			which reports the
	objective function values and associated computational times
	for the different numbers of design variables to optimize. This
	underlines the importance of identifying the most influent vari-
	ables of the optimization problem. In the following, we focus
	Dimension	min P	Time [h]
	2	0.006527	3
	5	0.006520	36
	9	0.006518	844

Table 1 :

 1 Best objective function values, at convergence, and corresponding computational times for different numbers of optimization variables.

Table 2 :

 2 Table2compares the angle of attack (AoA), minimum pressure coefficient Cp min and drag C d of the flexible and Rigid 30 kts hydrofoils at 30 kts of boat-speed. The Rigid 30 kts case has a lower incidence, with a drag reduction of roughly 1.2% compared to the flexible case. The two Cp min are relatively close and far from the criterion based on the cavitation number λ. The closeness of the two hydrofoils at this boat-speed of Comparison of the hydrodynamical characteristics of the Rigid 30 kts and flexible solutions at 30 kts of boat-speed.

		Rigid 30 kts	Flexible	Difference
	AoA	-0.36 •	-0.81 •	0.55
	Cp min	-0.605	-0.627	0.035
	C d	-	-	-0.012