Keywords: Reconfigurable Cable Parallel robot, Modulization, Optimization, Derivative-free method

Four Cables-Suspended Parallel Robot with

Reconfigurable platform H. EL GOLLI a , A.-A. KUMAR c , J-F ANTOINE b , G. ABBA a a. Ecole Nationale d'Ingénieurs de Metz, Université de Lorraine, LCFC, F-57000 Metz, France ; e-mails : hanene.el-golli@univ-lorraine.fr , gabriel.abba@univ-lorraine.fr b. IUT Nancy-Brabois, Dépt. GMP, Université de Lorraine, LCFC, F-54000 Nancy, France ; e-mail : jean-francois.antoine@univ-lorraine.fr, c. University of Luxembourg, 6, rue Richard Coudenhove-Kalergi, L-1359, Luxembourg; e-mail : atal.kumar@uni.lu

Résumé

Avec la progression scientifique des systèmes automatisés, les robots parallèles à câbles présentent des excellentes performances grâce à leur faible coût et de

Abstract:

The design of manipulators has progressed within few last decades. Recently, Cable-Suspended/Driven Parallel Robots are progressively used in industrial applications due to their large payload capacity, low cost, and good ability to perform several tasks requiring a large workspace. For these applications, the main objective is to pursue imposed trajectories. However, despite the advantage of the considerable workspace, the design of the suspended four-cables parallel robots provides an uncontrolled or limited rotational motion. In order to ensure control over the end effector's position and orientation, a self-reconfigurable four-Cable parallel robot with a compact structure is developed. For this purpose, the typical fixed anchor points will be translated on an inclined guide rail. Accordingly, this proposed work aims to select the appropriate anchor points which will be associated with each position of the desired trajectory. This suitable solution can be performed by an optimization approach by considering several criteria such as the robot equilibrium, cables tension, and the desired orientation.

Eventually, the numerical simulation of a reconfigurable robot is performed and analyzed. The results verify the feasibility of a 6-DoF for the parallel reconfigurable robot with four suspended cables.

Introduction

Cable-Driven Parallel Robots (CDPR) are parallel robots where the mobile platform is suspended by flexible cables instead of rigid elements [START_REF] Newman | Static analysis and dimensional optimization of a cable-driven parallel robot[END_REF]. This feature gives CDPR excellent performance in terms of large payload capacity, expanded workspace, accelerations, and speeds compared to conventional parallel robots [START_REF] Gosselin | Cable-driven parallel mechanisms: state of the art and perspectives[END_REF].

Recently, several CDPRs have been proposed and developed by the robotics community. These robots have been used in many engineering fields such as for urban construction tasks [START_REF] Pinto | A cable-driven robot for architectural constructions : a visual-guided approach for motion control and path-planning[END_REF], video recording applications [START_REF] Cone | Skycam: An Aerial Robotic Camera System[END_REF], and for radio telescopes [START_REF] Nan | The five-hundred-meter Aperture Spherical radio Telescope (FAST) project[END_REF]. Similarly, cable robots have demonstrated their adaptation for medical tasks such as haptic applications, mini-invasive surgery [START_REF] Boschetti | A haptic system for robotic assisted spine surgery[END_REF], and physical rehabilitation exercises [START_REF] Zanotto | Sophia-3: A semi adaptive cable-driven rehabilitation device with a tilting working plane[END_REF], and for industrial tasks such as motion simulators [START_REF] Miermeister | The CableRobot simulator large scale motion platform based on Cable Robot technology[END_REF], and manipulating parts in the manufacturing process [START_REF] Kumar | Design and Control of an underactuated cable-driven robot for agile handling of parts in a manufacturing line[END_REF]. However, most of these applications need single operation tasks which lead to the use of several CPRs with different configurations. Therefore, in the last few years, some research dealt with CDPRs which are modular and reconfigurable [START_REF] Gagliardini | A Reconfiguration Strategy for Reconfigurable Cable-Driven Parallel Robots A Reconfiguration Strategy for Reconfigurable Cable-Driven Parallel Robots A Reconfiguration Strategy for Reconfigurable Cable-Driven Parallel Robots[END_REF], [START_REF] Mishra | Optimizing Cable-Routing for Reconfigurable Cable-Driven Parallel Robots Optimizing Cable-Routing for Reconfigurable Cable-Driven Parallel Robots[END_REF] to be used in various situations. In this context, several works have dealt with manual reconfigurable robots such as reconfigurable cabledriven parallel robots (RCDPRs) ReelAx [START_REF] Mishra | Optimizing Cable-Routing for Reconfigurable Cable-Driven Parallel Robots Optimizing Cable-Routing for Reconfigurable Cable-Driven Parallel Robots[END_REF] and the NIST Robocrane [START_REF] Zhao | Typical configuration analysis of a modular reconfigurable cable-driven parallel robot[END_REF]. Other research deals with self-RCDPRs like SEGESTA prototype [START_REF] Bruckmann | A new force calculation algorithm for tendon-based parallel manipulators[END_REF], and the CoGiRo prototype [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable Suspended Parallel Robots[END_REF]. Nevertheless, all previous works, investigate RCDPRs by moving the attachment points associated to the fixed base's pulley. Therefore, to improve the mobility of this kind of CDPR without the need for more cables or multiple actuators, the self-reconfiguring task will be dependent on the anchor points attached to the moving platform.

In this work, we are interested in modeling the static behavior of the RCDPRs, then searching for a suitable position for the anchor point to have a secure motion during the process. The rest of the paper is structured as follows: Following this introduction, the description of the RCDPR is presented. The kinematic and static models of this type of parallel robot are formulated in section 3. Section 4 is dedicated to formulating and discussing the optimization results of the robot. Finally, the conclusion is presented in section five.

Structure Design

The proposed suspended CDPR is composed of a moving platform connected to a fixed base by four cables of variable lengths, used to control the pose of the end-effector. The cables are attached to the fixed structure of the robot at points 𝑨 𝒊 which represents the center of the attachment points, and connected to the moving platform at the anchor points 𝑩 𝒊 .

Referring to

Inverse Kinematic Model

The Inverse Kinematics model of the RCDPR determines the lengths of the cables, for a given pose 𝝌 and a given configuration of the moving platform.

The expression for calculating the cable length (li) which are considered as a straight-line segment is given by:

𝒍 𝒊 = ‖𝑨 𝒊 𝑩 𝒊 ‖ = 𝜌 𝑖 𝒅 𝒊 = ‖𝒑 -𝒓 𝒊 -𝒂 𝒊 ‖ (1)
Where 𝐫 𝐢 = 𝑟 𝑖 𝐁 𝐢 𝐂 𝐢 is the displacement vector of the anchor points.

Static Model

During the motion, external forces 𝑾 𝒑 = [ 0 0 -𝑚 𝑝 𝑔 0 0 0] 𝑇 are exerted on the platform due to its weight with payload and forces produced on task execution. Moreover, each cable produces forces on the platform proportional to its tension 𝑻 = [𝑻 𝟏 𝑻 𝟐 𝑻 𝟑 𝑻 𝟒 ] 𝑇 (cable weight is neglected). According to Newton-Euler formulation, the static equilibrium is denoted by:

𝕎 𝑻 + 𝑾 𝒑 = 𝟎 (2) Where 𝕎 = [𝑾 𝟏 𝑾 𝟐 𝑾 𝟑 𝑾 𝟒 ] 𝑾 𝒊 = [ 𝒅 𝒊 (𝑩 𝒊/𝑹 𝟎 -𝑷) × 𝒅 𝒊 ]
4

Reconfiguration Method 4.1 Planning Method

To reach the desired position, a new configuration is realized by translating two or more cables from their first anchor points to new locations. Accordingly, a scheme (Fig. 2) of a RCDPR control is given, to determine the optimal reconfiguration.

Optimization

Fig. 3: Flow chart of the proposed optimization method

An optimization approach for mechanism reconfiguration is proposed in Fig. 3. It aims to expand the Derivative Free Method allowing the search for the suitable configuration in each position. Starting from searching initial population 𝑰 𝟎 , the individuals are evaluated and then operators are applied to interpolate and evaluate the new population.

Design variables I

The design vector of the proposed approach is composed of three parts, cables tensions 𝑻, to avoid slacking during the displacement anchor points' displacement 𝑟 𝑖 , to select appropriate configuration, and orientation of the moving platform 𝜽 to improve pose control. Therefore, the design vector is written as: (𝑟 𝑚𝑖𝑛 = 0) ≤ 𝑟 𝑖 ≤ (𝑟 𝑚𝑎𝑥 = 𝑑 𝑚𝑥 ) And, the constraint imposed to guarantee better performance of the robot is to have always positive cables tension to avoid sagging cables. 0 < 𝑇 𝑚𝑖𝑛 ≤ 𝐓(𝑖) ≤ 𝑇 𝑚𝑎𝑥 Then, to examine if the desired orientation has occurred, the function to minimize error is used as follows (3)

𝑓 = min

Objective function

The criterion used in this approach aims to enhance the static equilibrium of the robot during its motion over the entire workspace. This criterion is mathematically formulated by the expression presented in section 3.2:

𝐹 1 (𝐼) = 𝕎 𝑻 + 𝑾 𝒑 (4)

Results

The feasible robot reconfigurations have been identified by mapping all the workspace. For each position, a configuration is selected, aiming to minimize the error between the calculated orientation of the platform and the required one. These optimal solutions are determined using MATLAB software. The algorithm was implemented with a function having a tolerance of 1e-16, and a maximum number of iterations equal to 10000. Fig. 4 shows the suitable position of the anchor points 𝑩 𝒊 . This figure presents a mapping value going from 𝑟 𝑚𝑖𝑛 to 𝑟 𝑚𝑎𝑥 along the plane (𝑿, 𝒀) for z = 0.5m, and (𝑿, 𝒁) for y=0.5m.

It should be noted that the anchor point position is near the center of mass of the platform, if the cable is the shortest length, inversely, the anchor point position is close to the border if the cable is of the maximum length. The motion of the robot is demonstrated by positive cable tension which is represented in Fig. 6 where the maximum tension can achieve 50N near the attachment point for that cable.

Conclusion

The first part of this work focuses on the fundamental characteristics related to the kinematic and static modeling of the new RCDPR. Then, an optimization approach specific to the position of the anchor points of the RCDPR is proposed in order to ensure stability and the maximum accuracy of the endeffector. This has allowed identifying the desired optimal solution. Results demonstrate that the new reconfigurable robot is most effective in terms of minimizing the inclusion of the end-effector during the displacement process with positive cable tensions and user safety.

Fig. 1 .Fig. 1 :

 11 Fig. 1:a. Geometric Parameters of the Ludic design b. Geometric Parameters of the platform

Fig. 2 :

 2 Fig. 2: Scheme of the optimal reconfiguration planning method

  𝑇 1 , 𝑇 2 , 𝑇 3 , 𝑇 4 , 𝑟 1 , 𝑟 2 , 𝛼, 𝛽, 𝛾] Where 𝑟 3 = 𝑑 𝑚𝑥 -𝑟 1 and 𝑟 4 = 𝑑 𝑚𝑥 -𝑟 2 and

Fig. 4 :

 4 Fig. 4 : Anchor points positions a. on the plan (𝑿, 𝒀) b. on the plan (𝑿, 𝒁)

Fig. 5 :Fig. 5 Fig. 6 :

 556 Fig. 5 : Platform Orientation on the plan (𝑿, 𝒀) Fig. 5 illustrates that the orientation of the platform is almost equal to the required orientation in the center of the workspace and the error of 0.2° on the boundary near the base's pulley.

Fig. 7 :Fig. 8 :Fig. 8

 788 Fig. 7: Desired trajectory


Keywords: Cycloidal propulsion, Experimental optimization, Multi-objective surrogate-based optimization

Abstract :

In recent years, innovative naval propulsion systems have been investigated thanks to the growing development of unmanned underwater vehicles. Cycloidal propellers are promising alternative concepts to the usual screw propellers. As bio-inspired technology, these systems use mechanical energy from unsteady hydrodynamic forces generated by blades oscillation like natural marine animal swimmers.

As an academic platform, the French Naval Academy Research Institute developed a large-scale experimental cycloidal propeller with the aim of running diverse pitch motions to evaluate performances of cross-flow propellers. Blades' pitching is here performed by servo-motors in order to control each blade independently. While common cycloidal propellers use mechanical blade actuators which restrain the blade motion possibilities, this blade-controlled platform allows new investigations of interesting research area in marine propulsion.

The platform is widely instrumented with load and torque sensors to measure instantaneous hydrodynamic forces during the rotation of the blades. Experiments, performed at the Ifremer current flume tank, first reveal that for classical sinusoidal pitch laws, performances are depending on the operating point (advance parameter λ) : the higher the advance parameter, the lower the sinusoidal amplitude must be for a better efficiency. These results are detailed in a previous paper written by Fasse et al. These results confirm the requirement of an adaptable pitch control for cycloidal propeller to boost their performances regarding the operation mode. Thus this paper deals with an experimental optimization, based on surrogates models (Efficient Global Optimization). This optimization is undertaken to surpass the performance of the propeller with parameterized pitch laws. This parameterization consists on modeling the pitch law by a spline defined through relevant control points corresponding to the optimization parameters. This method authorizes a wide range of possible motion taking account of the platform speed limits. Multi-objective optimization is performed for total thrust and efficiency maximizing. Results show that a trade-off is necessary between thrust and efficiency. Finally, Pareto Front is determined to get the best compromise of thrust and efficiency of the propeller.

Introduction

Cycloidal propellers are characterized by the rotation of several blades around an axis perpendicular to the advance direction (opposed to conventional propellers for which the rotation axis is generally aligned with the advance direction). This first rotation is called main rotation, given by the rotational speed Ω and related to the azimuth angle θ (Ω = dθ/dt). Each blade is also rotating around their own axis during the main rotation of the whole system (like helicopter's blades). This secondary rotation is mainly called pitch rotation, related to the blade pitch angle β or ϕ (depending of the definition of the angle). This particular motion creates strong unsteady hydrodynamic forces which produces lift and drag during the main rotation. The combination of the two rotations reproduces the motion of natural marine swimmers, especially studied by Triantafyllou et al [START_REF] Triantafyllou | Review of Experimental Work in Biomimetic Foils[END_REF]. For propulsive purpose, the part of hydrodynamic force in the direction of the ship advance produces the thrust force F X whereas the part in the perpendicular direction is called side force F Y . These propellers can easily orientate the thrust over 360 • by shifting the blade pitch law over the main rotation angle. For this type of propulsion, two kinematic modes are commonly defined according to the advance parameter λ which is the ratio of the ship advance speed V a and the blade peripheral speed V r = ΩR, R is the propeller radius:

λ = V a V r (1)
Epicycloidal mode is defined for λ<1 (Figure 1). For this mode, rotational speed is higher than advance speed. This mode produces a high thrust and strong manoeuvrability for low advance speed (limited by the maximal main rotation speed). On the other hand, trochoidal mode defined for λ>1 (Figure 2), is used to reach higher advance speed but with a lack of efficiency during starting phases (not adapted when V a is near zero). The blade motion is completely different between these two modes. Indeed for epicyloidal mode, the blade chord is roughly following the tangent of the main rotation (because ΩR is higher than V a ), whereas for trochoidal mode the blade oscillates around the advance direction. As an interesting marine propulsion system, cycloidal propeller have been studied by the IRENav with the developpment of an blade-controlled plateform [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF]. Thanks to an electrical blade-command design, this platform allows to reproduce blades' kinematic from all conceivable movements (trochoidal or epicycloidal kinematics). This electrical blade-command also allows authors to perform pitch law optimization to improve performances of current cycloidal propellers. The aim of this article is to present an experimental optimization using a Gaussian process based method, coded by Sacher [START_REF] Sacher | Efficient optimization procedure in non-linear fluid-structure interaction problem: Application to mainsail trimming in upwind conditions[END_REF], to maximize both thrust and efficiency of cycloidal propellers. The optimization uses experimental measurements collected from the instrumented blade-controlled patform. Experiments are performed at the Ifremer current tank. These experimental facilities are described in the first section as well as the parameterization of the pitch laws. The second section deals with the optimization method overview. Finally results of sinusoidal laws are summarized and compared to the results of the optimization for λ=1.2 simultaneously from thrust and efficiency aspects.

2 Optimization problem

Context

Cross-flow propeller, also called cycloidal propeller, generates a vectorial thrust which can be adapted in norm and direction. The main control parameter of this type of propulsion is the way the blades are moving through the incoming flow which is called pitch law. Indeed, cycloidal propeller uses lift forces generated by the blades through the flow during the whole rotation. So the pitch law is essential to the design and the performance of a cycloidal propeller, therefore it is the purpose of the optimization presented in this paper. Most of existing cross-flow technologies use a mechanical system to operate the blade pitch motion as is the case for the renowned Voith-Schneider Propeller [START_REF] Schneider | Blade Wheel[END_REF] for which blades are actuated by a rod/crank mechanism. More recently Roesler [START_REF] Roesler | Experimental performance of a novel trochoidal propeller[END_REF] designed an experimental platform with a system of crankshaft and concentric collars to generate trochoidal motions and they performed measurements at the UNH testing facility in USA. Mechanical system implies a symmetrical blade pitch law to ensure the stability and the load balancing of the propeller. This is also the simplest way to design a mechanical cycloidal propeller. But because of the complexity of the flow and the high perturbed wake, the blade behaviour has no reason to be symmetrical between up and downstream flow. Indeed during the upstream halfrotation, the flow viewed by the blades is calm and undisturbed, whereas in the downstream half-rotation the blades encounter the wake of upstream other blades regularly agitated by vortex or free surface suction phenomenon.

To deal with this challenge, authors choose to go through an experimental optimization of the cycloidal propeller pitch law using a non-symmetrical parameterization and an experimental platform able to run this kind of law.

Blade-controlled platform

To perform the optimization, the experimental blade controlled cross-flow propeller, called SHIVA developed by the authors, is operating at the Ifremer wave and current flume tank. Figure 3 gives an overview of the platform and its positioning at the flume tank.

SHIVA is composed of a triangular frame [F] in which three blades [I] are mounted. These blades are rotating thanks to three servo-motors [E] that gives them independent motions. The outer frame [D] allows the platform to be elevated above the tank in order to dive just the blades Each blade is independently actuated by a 0.8kW auxiliary motor [E]. The blades rotate around the quarter of their chord, close to the aerodynamic center for the NACA 0018. This rotation is called blade pitch and is referred to as the oriented angle φ between the tank flow direction (X-axis) and the blade chord. The pitch angle can also be defined as the oriented angle β between the ortho-radial line of the main rotation and the blade chord (see Figure 4). These two pitch angles are linked by the relation θ = φ -β. The azimuth position of a blade θ is the position of its chord quarter on the main rotation disk, and for an arbitrary time ∆t, θ=Ω∆t where Ω is the main rotational speed. Expression of the pitch angle in function of the azimuth position is the pitch law. These laws are tabulated as discrete functions on an embedded micro-controller unit which send regulation orders directly to the speed controllers for each blades in order to follow the reference pitch law. These laws are loaded on the platform via Ethernet liaison from a computer located alongside the tank. Pitch laws can be loaded even during the rotation which helps to gain precious time during optimization procedure.

Because the pitch motion is electric and no more mechanical, a delay of the real blade position can occur when the motors can't follow the reference pitch according to their limits in speed and torque. Figures 5 gives the pitch tracking for two laws: a sinusoidal law at the left and an epicycloidal law at the right. These pitch tracking measurements have been performed during in water experiments with a tank speed V =0.8 m/s and a rotational speed Ω set for the desired λ value (Ω=8.68 RPM corresponding to λ=2.2 and Ω=21.2 RPM corresponding to λ=0.9). The pitch law for trochoidal mode is given by the φ angle (pitch around the inlet flow direction) whereas the pitch is expressed by β angle in epicycloidal mode (pitch around the tangent of the main rotation). The left axis (in red) of each sub-figure gives the difference between the reference and the real pitch. These results show that for trochoidal mode a maximal error of 1 • occurs when the gradient (dφ/dθ) is highest. For epicycloidal mode, there is a blade reversal at θ=180 • , which characterizes epicycloid motions. This quick motion is not well following by the motors because of the speed limits (the error reaches almost 100 • ). After this huge error, the blade keeps a delay which is never caught. A detailed study of the pitch tracking error is given by Fasse in [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF]. These limitations in speed involve the use of a law's classification during the optimization procedure (See section 2.3). Indeed, if the pitch tracking error is too high, the pitch law is classified "non-doable". Moreover, with current motorization only trochoidal laws are performed on the platform SHIVA during experiments.

In addition, SHIVA is widely instrumented to measure time-dependant hydrodynamic loads and blade angular positions. Rotary incremental encoders provide angular position of each blade with a precision of 0.088 • . They are used simultaneously for the control loop as feedback and for loads measurements as projection angle. The azimuth position is also recorded by a rotary incremental encoders with a precision of 0.014 • . One blade is instrumented with an embedded 5-components load-cell [G] which measure local hydro-dynamic loads in the direction of the blade's chord T b (measuring range of 200 N) and normal to the blade N b (measuring range of 900 N). This load-cell also measures the torque C b applied on the blade by hydrodynamic loads (measuring range of 60 Nm). Four fixed 3-components load-cell [H] located between the outer frame I-beams and the tank I-beams provide global solicitations of the propeller in the X,Y and Z directions (measuring range of 5000 N in the X and Y direction and 7000 N for the Z direction). At least, a torque sensor [B] is installed between the main motor and the speed reducer and measures the torque applied on the main driveshaft C tot (measuring range of 20 Nm). Due to the wider measuring range of these fixed load-cell, in comparison with the embedded load-cell, post-processing of the loads is done with local hydrodynamic loads (for more precision, refer to [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF]).

The embedded load-cell measures temporal voltages which are highly disturbed by high frequencies of electromagnetic environment mainly due to the Pulse Width Modulation of speed controllers. Low-pass filter (50Hz cut-frequency) is thus applied on temporal signals. Then voltages are converted into loads thanks to the transfer matrix determined by calibration. Finally, thanks to the acquisition of angular positions, forces in blade reference frame (F N b , F T b ) are projected into rotor reference frame (F N , F T ) and absolute reference frame (F X , F Y ) by equations ( 2) and ( 3).

F T F N = cos(β) -sin(β) -sin(β) -cos(β) F T b F N b (2)
F X F Y = cos(θ) sin(θ) -sin(θ) cos(θ) F T F N (3)
Then instantaneous forces are phase-averaged over at least 30 revolutions to one period, shifted of π 3 and 2π 3 , and summed to give the total force acting on the three blades. Lastly, total thrust coefficient, torque coefficient and hydrodynamic efficiency are calculated as the non-dimensional averaged values:

C F X = F Xtot 0.5ρSV 2 C T = C tot 0.5ρSDV 2 η = F Xtot V C tot Ω (4)
Where F Xtot and C tot are respectively the mean of total thrust and total torque (measured on the main torque sensor) over 360°, S is the swept frontal area (S = Dl) and V the flume tank flow velocity. These experimental performances are measured and spent through the optimization code to evaluate the objective function. To perform the optimization, the choice of the parameters have to be suited for maximizing the performance of cycloidal propeller. The blade motion is directly linked to the hydrodynamic force generation. And because SHIVA is designed to run all type of blade pitch motion, the parameterization consists in configure the pitch laws with several parameters.

Parametric model

To proceed the optimization procedure, the first step is to develop an accurate parameterization. Indeed, Efficient Global Optimization methods, used here, allows to decrease the number of optimization step in comparison with genetic methods. But more the number of parameters is high and more the number of optimization step is required. The method presented in this paper consists in defining the pitch law by a spline passing through four control points which are given by 3 parameters (inspired from Abbaszadeh [START_REF] Abbaszadeh | Experimental investigation of an optimised pitch control for a vertical-axis turbine[END_REF]). Figure 6 illustrates the parameterization for two random laws. The four control points are the first black triangle, the red square, the blue circle and the last black triangle. These points are defined with three parameters:

1. The first parameter x 1 gives the abscissa of the two central points (red square and blue circle in the Figure 6) by setting the distance of this abscissa and the central position θ mean :

x red square = θ mean -x 1

x blue circle = θ mean + x 1 2. The second parameter x 2 give the ordinate of these two central points:

y red square = x 2 ; y blue circle = -x 2
3. The third parameter x 3 allows an offset of the law for the zero crossing at the beginning and the end of the spline (the last black triangle is the 2π-offset of the first one).

x black triangle = x 3 ; y black triangle = 0

For the presented results, the central position θ mean is set at 170 degrees. This position involves a non-symmetrical pitch law which is the intent of this optimization as explained in the context. The nonsymmetrical construction of the law also implies a different gradient between the start (θ=0 • ) and the end of the spline (θ=360 • ). To fix this problem, the spline is defined on 10 periods (5 before and 5 after) to ensure the continuity at θ=0 • and θ=360 • (with an error of 10 -7 ).

Each parameter can move through the following ranges:

• ∆x 1 = [0.5 : 120]

• ∆x 2 = [0.5 : 80]

• ∆x 3 = [-5 : 40]
This method of parameterization implies a weakness which occurs when two points of the spline are to close, the amplitude of the spline can be very high. To avoid this, a classification is used and taking into account during the optimization procedure (detailed in Section 3.2). This classification consists on defining two limits (high and low) to the spline and classify the parameterized spline "doable" if all the points of the spline are between these two limits. These limits are represented in green on the left of the Figure 6 and fixed at φ = 50 • and φ = -80 • for the present experiments. These limits have been chosen after experiments considerations: high amplitude laws imply an acceleration for auxiliary motor which is to high (refers to Section 2.2) and in the half downstream it has been sensed that the pitch must be higher than the upstream to give the ability of the blade to "cut" the wake.

For hundred arbitrary laws have been generated over the range of the 3 parameters to estimate the classification model. The right sub-figure of Figure 6 represents this classification model with the iso-surface of 50%-probability of the law to be "doable". This figure shows in the front the area of laws which are considered "doable" and in the back those which have high probability to be "non-doable". During the optimization process, the classification is taking into account and so the area of interest for the new optimized point is always located inside the front area of 50%-probability to be "doable" (see Section 3.2). In addition, if a new optimized law presents a pitch error too high, it is classified "non-doable" and the classification model is updated.

Multi-objective Surrogate-based Optimisation

We consider abstract optimization problems, where several objectives have to be minimized simultaneously over a design variable Ω ⊂ R d :

min x∈Ω f 1 (x), . . . ,f m (x) s.t. g(x) ≤ 0 (5)
where x = (x i ) 1≤i≤d is a vector of design variables, f = (f j ) 1≤j≤m is a vector of objective functions to be minimized (f i : Ω → R), and g = (g k ) 1≤k≤p is a vector of inequality constraints (g k : Ω → R).

The existence of an optimal solution, minimizing all objectives at once is usually not granted. This leads to the search for an optimal set of solutions, called the Pareto front.

In the present work, individual evaluations of the objective functions are assumed to be very expensive (experiments results). A common approach is then to use surrogate models in place of f to reduce the computational burden related to the evaluations of f [START_REF] Simpson | Metamodels for Computer-based Engineering Design: Survey and recommendations[END_REF]. Gaussian processes (GP) [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF] are presently considered, which, owing to their statistical nature, provide for each objective function, both a prediction value and a measure of the uncertainty (variance) in this prediction. These features are appealing in the optimization context, as they can be exploited to derive rigorous optimization strategies, by evaluating sequentially the computer models at design vectors that maximize a so-called merit function [START_REF] Picheny | A benchmark of kriging-based infill criteria for noisy optimization[END_REF].

In mono-objective problems, GP-based approaches are globally referred to as Efficient Global Optimization (EGO) [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF], where the merit function is based on an Expected Improvement (EI) criterion, expressing a trade-off between sampling in promising regions and exploring in unsampled regions. EGO has been successfully applied to complex optimization problems, such as non-linear fluid-structure interaction problems [START_REF] Sacher | Efficient optimization procedure in non-linear fluid-structure interaction problem: Application to mainsail trimming in upwind conditions[END_REF][START_REF] Sacher | Flexible hydrofoil optimization for the 35th America's Cup with constrained EGO method[END_REF] or RANS computations [START_REF] Meliani | Multi-fidelity efficient global optimization : Methodology and application to airfoil shape design[END_REF].

Over the last few years, surrogate-based approaches have also been proposed to address the multiobjective problem. Indeed, several GP-based multi-objective optimization strategies propose to extend the EI infilling criterion of Jones [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]. The definition of the improvement can then be considered over the Pareto front, using for instance, Maxi-min distance [START_REF] Svenson | Multiobjective optimization of expensive black-box functions via expected maximin improvement[END_REF], Euclidean distance [START_REF] Keane | Statistical improvement criteria for use in multiobjective design optimization[END_REF], or the well-known Hypervolume infilling criterion [START_REF] Luo | Kriging model based many-objective optimization with efficient calculation of expected hypervolume improvement[END_REF]. In the present work, the latter criterion is used to solve the blade pitch law optimization problem.

In Section 3.1, we summarize the construction of a GP model for an objective function f j . Definitions of multi-objective surrogate-based merit functions, that are used in the present optimization under nondoable observations, are then provided in Section 3.2.

GP surrogate construction

In the present Section 3.1, the GP surrogate construction of an objective function f j is provided. The method is exactly the same for the GP surrogate construction of an inequality constraint function g k .

We consider a set of n training points X = (x 1 , . . . ,x n ), each in Ω. The training points are associated to the vector Y (j) = (y

(j) 1 , . . . ,y (j)
n ) of noisy observations of the objective function f j . It is assumed that y

(j) i = f j (x i ) + ϵ (j)
i , where the ϵ (j) i are independent and identically distributed Gaussian random variables having zero-mean and variance σ 2 ϵ j . The GP construction considers that f j (x) is a realization of a zero-mean multivariate Gaussian process with covariance function C j . We consider here the multidimensional squared exponential covariance functions defined by:

C j x,x ′ ; Θ j . = θ 1 d i=1 exp -(x i -x ′ i ) 2 2l 2 i + θ 2 (6)
where

Θ j = θ (j) 1 , θ (j) 2 , l (j)
1 , . . . ,l (j) d

is the vector of covariance hyper-parameters to be inferred from the Y (j) observations. From the conditional rules of joint Gaussian distributions [START_REF] Williams | Gaussian processes for machine learning[END_REF], the best prediction f j (x) of f j (x), i.e. the mean of y (j) , and the prediction variance σ j 2 (x) are given by:

f j (x) = k T j (x) C j (Θ j ) + σ 2 ϵ j I -1 Y (j) (7)
σ j 2 (x) = κ j (x) + σ 2 ϵ j -k T j (x) C j (Θ j ) + σ 2 ϵ j I -1 k j (x) (8)
In ( 7) and ( 8) we have denoted C j ∈ R n×n the symmetric covariance matrix of the training points,

κ j (x) . = C j (x,x; Θ j ), k j (x) . = (C j (x,x 1 ; Θ j ) • • • C j (x,
x n ; Θ j )) T the covariance vector between the observations in X and x, and I the identity matrix of R n . The hyper-parameters Θ j and noise variance σ 2 ϵ j can be determined by maximizing the log-marginal likelihood (see [START_REF] Williams | Gaussian processes for machine learning[END_REF] for more details).

Merit functions under non-doable design vectors

Surrogate-based optimization methods rely on the sequential construction of statistical surrogate models, using training sets of computed objective and constraint function values, that are refined according to a prescribed infilling strategy (i.e. merit functions [START_REF] Picheny | A benchmark of kriging-based infill criteria for noisy optimization[END_REF]). At each iteration of the iterative surrogatebased optimization, a new design vector x n+1 is thus added to X , and finally f and g are computed. A new iteration can then start by updating surrogate models, and the iterative process is repeated until a stopping criterion is satisfied or the resources allocated to the optimization have been exhausted.

In real world applications, however, this sequential optimization procedure may stop prematurely if an objective function cannot be computed at a proposed x n+1 . Indeed, such a situation may occur when the search space encompasses design points corresponding to a non-physical configuration or an ill-posed problem, for example. Sacher [START_REF] Sacher | A classification approach to efficient global optimization in presence of non-computable domains[END_REF] have proposed to use a classification model to learn non-doable design vector areas alongside the surrogate-based optimization and adapt the infilling strategy accordingly. In short, a probabilistic classification model is built using the union of doable and non-doable training sets. The classifier is then incorporated in the surrogate-based optimization procedure to avoid proposing new design vectors in the non-doable domain while improving the classification uncertainty if needed.

In the present work, two new design vectors are determined under the constraint to be in the doable domain and added to X at each new iteration of the iterative surrogate-based optimization. Specifically, the EIMH and MCPF criteria are considered (see next Sections 3.2.1 and 3.2.2) to allow respectively, the improvement of the computed Pareto front and the search of the best compromise in this Pareto Front.

Expected Improvement Matrix Hypervolume criterion (EIMH)

Multi-objective EI criteria have usually similar processes, where in first, an improvement function over the Pareto front is defined. In second, the expected value of the improvement function is computed by integrating it over the non-dominated region. Zhan [START_REF] Zhan | Expected improvement matrix-based infill criteria for expensive multiobjective optimization[END_REF] have proposed the concept of Expected Improvement Matrix (EIM) to address the usage of multi-objective EI criteria in real world applications, especially when dealing with many objective problems. This EIM concept uses a combination of t × m simple 1-D integrations to compute a cheap-to-evaluate and still efficient multi-objective EI (t is the number of non-dominated design vectors and m is the number of objectives). In [START_REF] Zhan | Expected improvement matrix-based infill criteria for expensive multiobjective optimization[END_REF], the scalar Expected Improvement function EI(x) is expanded into a 2-D matrix EIM, in which,

EI q j (x) . = f q j -f j (x) Φ f q j -f j (x) σ j 2 (x) (9)
+ σ j 2 (x)ϕ f q j -f j (x) σ j 2 (x)
where j = 1, . . . ,m and q = 1, . . . ,t. The element EI q j (x), in EIM, corresponds to the EI of the studying point x beyond the qth non-dominated front point in the jth objective. From the EIM matrix, Zhan [START_REF] Zhan | Expected improvement matrix-based infill criteria for expensive multiobjective optimization[END_REF] have then proposed three infilling criteria (Euclidean distance, Maximin distance and Hypervolume improvement) to aggregate the EIM into a scalar value to measure the overall improvement of the studying point compared against the Pareto front. In the present work, the Hypervolume-based EIM criterion is applied, leading to the selection of a design vector variable x EIM H . For more details on the EIMH criterion, interested readers are referred to [START_REF] Zhan | Expected improvement matrix-based infill criteria for expensive multiobjective optimization[END_REF].

Median Compromise of the Pareto Front (MCPF)

Solving multi-objective optimization problems leads to the search of the optimal set of solutions (i.e. Pareto front). In most real-word applications, a design vector from the Pareto set has to be selected to provide an optimal design choice. This selection is thus a compromise choice, which is usually based on the knowledge of the designer. With the aim of avoiding this designer compromise selection, the MCPF criterion is used in the present work, and is given by:

x M CP F . = arg min x∈X P F   m j=1 arg sort f P F j , f j (x) 2   (10)
where, arg sort(f ,f ) . = {q ∈ {1, . . . ,t} | f < f q ∈ f } and f P F j is the vector of objective function j predictions on the Pareto set X P F , which is composed of t non-dominated design vectors. As can be seen in Eq. ( 10), the median trade-off is selected by considering the rank, from the Pareto set, of each objective function prediction. This definition allows a design vector compromise to be determined without the need for weights or scaling of objective function.

Experimental measurements and optimization results

First experiments are undertaken at the Ifremer flume tank for sinusoidal pitching laws with the SHIVA platform. The measurements have been done for various advance parameter values to evaluate the performances of the cycloidal propeller over a range of operating conditions. Then the optimization procedure is carried out for λ=1.2.

Sinusoidal pitch results

All the results with sinusoidal laws are measured with the following experimental conditions: configuration of SHIVA with a diameter D=0.8 m and blade chord of c=0.35 m leading to a solidity of σ=2.625. The tank speed V is set at 0.8 m/s corresponding to a global Reynolds number Re = V D ν =640 000. The rotational speed Ω is then adjust to the desired λ value (see the Table 1). Measurements are registered on 30 rotations of the platform after waiting that the flow is well established. Blade-control parameters are adjusted in function of the rotational speed to ensure a pitch track with a maximal error under 1 degree. Signals from the instrumentation are converted in force and projected to calculate the thrust force and main torque. These signals are then phase-averaged over the 30 periods and non-dimensionalized according equations 4. More details of the signals treatment for these experiments is providing by Fasse [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF].

λ 0.8 1 1.2 1.4 1.6 1.8 2
Global results of sinusoidal laws are given in the Figure 7. The thrust coefficient is represented by colored dotted lines with square markers whereas torque coefficient is represented by colored dashed lines with star markers. Both variables are referred with the left axis. The efficiency is represented by colored lines with circle markers and referred with the right axis. Results reveal that depending on the value of the advance parameter, each pitch law is interesting with a maximal efficiency reached for different values of λ. This first conclusion corroborates the benefit of an adaptable pitch system to adjust the pitch to the desired operating mode. It also means that during the acceleration phase, it is preferable to start with high amplitude pitch law then decrease the pitch amplitude. However, the results for the Sinus40 law show that for λ<1, a trochoidal law is not well-suited for the epicycloidal mode. The transition between epicycloidal and trochoidal mode is complex in term of pitch law shape and requires a specific investigation that is not detailed in this paper. Furthemore, because motors of the platform have encountered some limitations in epicycloidal (explained in Section 2.2), the study presented here is focused on the trochoidal mode. Concerning the thrust, the Sinus10 law has the higher values for all the advance parameters. It shows that for manoeuvre which requires high thrust, the pitch amplitude must be decreased, balancing a reduction Sinusoidal pitch law are suited for mechanic pitch systems as periodic and symmetrical blade motion (as studied by Bose [START_REF] Bose | Experimental performance of a trochoidal propeller with high-aspect-ratio blades[END_REF] and Roesler [START_REF] Roesler | Experimental performance of a novel trochoidal propeller[END_REF]). But there is no reason why the pitch law definition should be symmetrical. Indeed, the complexity of the flow through cross-flow propellers implies a strong discrepancy between upstream and downstream halves. It is why an experimental optimization is proceeded with a relevant parameterization of the pitch law (detailed previously).

Optimized pitch results

Experimental conditions for the optimization are as follows: λ=1.2, tank speed V =0.8 m/s, SHIVA rotational speed Ω=15.92 RPM. Before tackle the optimization, a sampling of 35 pitch laws is realized through the parametric space using Latin Hypercube Sampling method [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. Then approximately 50 optimization steps have been performed during 8 hours of experiments (approximately 10 minutes by step). One step of optimization consists on:

1. Generate the pitch law with the three parameters (x 1 , x 2 , x 3 ).

2. Upload the pitch law on the SHIVA platform and start the measurement of force signals after waiting few minutes that the flow is well established.

3. After a duration corresponding to 30 rotations of SHIVA, measurements are post-treated to calculate the value of objective function (here is the mean thrust coefficient C F X and the efficiency η)

4. These new data are computed to update the two meta-models of thrust and efficiency and the Pareto Front.

5. The optimizer gives back two new triplet of parameters after computing the EIMH and the MCPF criterion to respectively search for the best improvement for both objectives and search for the best compromise on the Pareto Front.

6. Return to the stage one. Figure 8 gives the two meta-models in the 3d parametric space at the end of optimization. Octahedrons are the experimental points located in function of the parameters which were used to generated the pitch law. As mentioned earlier, all experimental points are inside the front area delimited by the black 50%probability iso-surface. Octahedrons are colored with their objective value (total thrust coefficient at the left and efficiency at the right). Slices show the meta-model inferred from experimental points and so colored in function of the objective value. Lastly, for respectively the thrust and the efficiency, the red iso-surface represents the pockets of the meta-model for which the thrust coefficient is higher than 0.75 (absolute value) and respectively higher than 0.63 for the efficiency. Because the optimizer uses a minimization procedure during the optimization, objective function are given with negative sign. These results show that the pocket of maximal thrust is not located in the same area than this for efficiency. It proves that for λ=1.2, the blade motion which produces the maximal thrust is not this which is the best efficient. The law which maximize the thrust is generated by the triplet (48.88, 7.37, 39. [START_REF] Zhan | Expected improvement matrix-based infill criteria for expensive multiobjective optimization[END_REF]) whereas this which maximize the efficiency is given by the triplet (37.07, 17.92, 8.93). According to the value of parameter x 2 , which deals with the amplitude of the law, optimization shows that the law with the best efficiency has an amplitude around 20 • whereas for thrust maximizing the amplitude of the law is lower. This result confirms observations about sinusoidal laws for which Sinus20 have a better efficiency at λ=1.2 than Sinus10 but a lower thrust coefficient (see Figure 7).

These two meta-models reveal that there is a compromise between thrust maximization and hydrodynamic efficiency. The Median Compromise of the Pareto Front criterion, which is based on ranking trade-off between multiple objectives, permits to search for the best law satisfying both objectives without regarding objective values but only the rank of each objective through all measured laws. The best compromise is the law for which the sum of the rank is minimum. Thus no weighting of objectives is necessary.

Figure 9 presents the Pareto Front obtained from the experimental optimization. The little colored circles correspond to each generated laws and are plotted in function of the experimental results in thrust coefficient and efficiency pairs. Black squares correspond to a virtual front estimated by the optimizer based on the experimental results. Finally, sinusoidal results for λ=1.2 are displayed to compare performances.

As expected, there are significant performance benefits to be gained with optimized law for both objectives. The efficiency gain is around 8% between Sinus20 and the Pareto Front (η Sinus20 = 0.57 and η Optim = 0.64) while the advantage in thrust is over 11% between Sinus10 and the Pareto Front (CF X Sinus10 = 0.67 and CF X Optim = 0.78).

Figure 10 gives general results (pitch law φ and blade forces F X and F Y ) for the 20-sinusoidal law (in dark blue) and three optimized laws: with the best thrust (in orange), the best efficiency (in red) and the best compromise (in light blue). The right sub-figure shows the experimental measurement of hydrodynamic forces along θ positions (each 30 • ). These three laws are also represented in the Figure 9 as well as the Sinus20 law. Concerning the pitch law, optimization procedures reveal two main differences with sinusoidal laws:

• Non-zero pitch at θ=0 • (dephasing effect)

• Higher pitch for downstream positions than for upstream positions (non-symmetrical effect)

The first effect is due to the x 3 parameter which forces the law to cross the φ=0-line at θ=x 3 (θ=39.19 • for the thrust maximization law and θ=8.93 • for the efficiency maximization law). While sinusoidal law passes through the point (θ=0,φ=0), optimized laws have slight incidence at θ=0 • . This leads to a reduction of the hydrodynamic forces near θ=0 • especially side force component F Y . For these positions, thrust force is close to zero, so the reduction of side force have a favorable effect on the efficiency. Indeed, by decreasing the hydrodynamic loads at these positions the total torque required to turn the whole platform is also reducing. The non-symmetrical effect has two different impacts regarding the optimized laws. For the F Xmax law, the thrust peak in the upstream is shifted towards θ=90 • (whereas the peak is reached near θ=120 • for the two other laws). In the downstream half, the thrust peak is even higher than the upstream half, contrary to the η max and the Sinus20 laws for which the downstream thrust is lower than in upstream. But the thrust increasing comes with a side-force inflation which induces a diminution of efficiency. For the η max law, the non-symmetrical effect leads to significant decrease the side-force F Y in the downstream half in comparison with the Sinus20 and other optimized law, while thrust force is quite similar to the Sinus20 law. And so the required torque to the main rotation is lower than the classic sinusoidal law which implies a better efficiency. Finally the law which have the best compromise between thrust and efficiency has a pitch shape located between F Xmax law and η max law. Concerning the thrust, this law is similar to the F Xmax law in the upstream and almost for the downstream. But the significant gain is on the side-force which is clearly lower than for the F Xmax law, resulting to a better efficiency. This trade-off between efficiency and thrust can be modified into the search of the law giving the best efficiency for a given value of thrust. Indeed, for propulsion purpose, the idea is to reach a desired advance speed while overcoming the hull resistance. So knowing this hull resistance, which depends on the hull shape, dimensions and the advance speed, this procedure is very interesting for the design of vertical axis propeller. To go further, a generalization of this optimization for different λ values ensures an optimal law for a wide range of velocities and operation modes.

Conclusion and Perspectives

This paper presented an experimental optimization procedure of a cycloidal propeller blade pitch laws. Indeed, for this specific type of propulsion, characterized by the rotation of blades around an axis perpendicular to the ship advance speed, the blade motion has a huge impact on the propeller performances. Traditional cycloidal propellers have a mechanical system which deals with the blade pitch rotation. These systems, represented by sinusoidal pitch law in trochoidal mode (λ >1 systems), show limitations in term of efficiency. A fully electric system, as the SHIVA platform described in this paper, can surpass the performances by managing non-symmetrical blade pitch laws. By using a suitable parameterization of the pitch laws allowing a wide diversity of pitch motion, the multi-objective optimization is performed on experimental measurements of the SHIVA platform at the Ifremer current tank. Both thrust and efficiency maximization are chased to satisfy two different objectives. Results shows that a trade-off is necessary concerning the pitch law shape to maximize the thrust or the efficiency. A Pareto Front is thus determined between these two objectives to search for the best compromise law. The results and the comparison with sinusoidal results show a significant improvement of both thrust and efficiency for the best compromise. Finally, because experimental optimization takes time, the procedure have been carried out for only one λ value. As perspectives, it will be interesting to repeat this optimization for a range of λ values to compare the performances of the cycloidal propeller for many operating points. Computational Fluid Dynamic calculations using FINE/Marine commercial code (2D Unsteady Reynolds Averaged Navier-Stokes) are currently performed with the experimental facilities dimensions. The same optimisation procedure could be numerically achieved to compare the optimization results. An other perspective is to upgrade the SHIVA platform with the installation of a Artificial Intelligence to deal with the automatizing of the optimization. Indeed, here the optimization is performed step by step during long measurements. But because the platform is instrumented with embedded load sensors, it could process an automatic optimization for various experimental conditions (such as diverse tank or rotational velocities).

Figure 1 :Figure 2 :

 12 Figure 1: Epicycloidal path of a blade

Figure 3 :

 3 Figure 3: Blade-controlled platform and its positioning at the current flume tank

Figure 4 :Figure 5 :

 45 Figure 4: Representation of angles for an arbitrary position θ and an arbitrary relative velocity W

Figure 6 :

 6 Figure 6: At the left, pitch law parameterization method. At the right, classification iso-surface of "doable" pitch law in the parametric space.

Figure 7 :

 7 Figure 7: Global performances (Thrust and Torque coefficients at the left scale and Efficiency at the right scale) for four amplitudes of sinusoidal laws.

Figure 8 :

 8 Figure 8: Meta-models for Thrust (left) and Efficiency (right) at the end of multi-objective optimization.

Figure 9 :

 9 Figure 9: Pareto front obtained from experimental optimization and compared to sinusoidal pitch results.

Figure 10 :

 10 Figure 10: Detailed results of hydrodynamic forces generated by three different laws (Sinus20, Optimized law maximizing the trust F Xmax and Optimized law maximizing the efficiency η max ) over azimuth positions.

Optimization of a marine cycloidal propeller using experimental blade-controlled platform

	G. FASSE a,1 , M. SACHER b , F. HAUVILLE a , G. GERMAIN c
	a. Naval Academy Research Institute (IRENav), BCRM Brest, CC600, 29240 Brest Cedex 9, France
	b. ENSTA Bretagne, CNRS UMR 6027, IRDL, 2 rue Francois Verny, 29806 Brest Cedex9, France
	c. Ifremer, Hydrodynamic & Metocean Service F-62200 Boulogne-sur-Mer, France
	1. Mail: guillaume.fasse@ecole-navale.fr

Table 1 :

 1 Ω (RPM) 23.87 19.09 15.92 13.64 11.94 10.61 9.55 8.68 7.96 Corresponding table between the advance parameter λ and the rotational speed Ω.

	.0	2.2	2.4


Rajasekar Gopalsamy
email: rajasekar.gopalsamy@univ-eiffel.fr

Nicolas Chevaugeon
email: nicolas.chevaugeon@ec-nantes.fr

Olivier Chupin
email: olivier.chupin@univ-eiffel.fr

Jean-Michel Michel Piau
email: jean-michel.piau@univ-eiffel.fr

Ferhat Hammoum
email: ferhat.hammoum@univ-eiffel.fr

Hammoum Ferhat

Lip-field approach for modeling strain softening in viscoelastic materials

Keywords: fracture, localization, damage, viscoelasticity, Lip-field

come   L'archive ouverte pluridisciplinaire

Introduction

Viscoelastic materials like asphalt, biological tissues,wood and polymers have numerous applications in engineering. Almost all materials exhibit viscoelastic behavior when the considered time scale is large enough. The properties of these materials are highly rate and temperature dependent that hence largely affect their failure behavior [START_REF] Lakes | On Poisson's Ratio in Linearly Viscoelastic Solids[END_REF][START_REF] D'amico | Moving cracks in viscoelastic materials : Temperature and energy-release-rate measurements[END_REF]. Design for mitigation of fracture in viscoleastic materials is an important problem for e.g. in case of asphalt and pavement construction [START_REF] Pirmohammad | Fracture Behavior of Asphalt Materials, Structural Integrity[END_REF].

A detailed review of viscoelastic fracture was presented by Knauss [START_REF] Knauss | A review of fracture in viscoelastic materials[END_REF]. Several works have been done in the past for the modeling of fracture in viscoelastic solids. Schapery [START_REF] Schapery | Correspondence principles and a generalized J integral for large deformation and fracture analysis of viscoelastic media[END_REF][START_REF] Schapery | On some path independent integrals and their use in fracture of nonlinear viscoelastic media[END_REF] used the correspondence principle [START_REF] Brinson | Polymer Engineering Science and Viscoelasticity[END_REF] to define a generalized time dependent J-integral for fracture in viscoelastic material analogous to the elastic case. Based on this theory, a new model was proposed in [START_REF] Dubois | Modelling of the crack growth initiation in viscoelastic media by the Gθv-integral[END_REF] to study the creep fracture in viscoelastic materials.

Several numerical methods have been developed to study the failure behavior of different materials caused by crack propagation. They can be broadly classified into discrete fracture mechanics based approach and damage mechanics based approach. In discrete approach, X-FEM [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] has gained a lot of attention in the past. X-FEM approach was used in [START_REF] Zhang | Modeling inclusion problems in viscoelastic materials with the extended finite element method[END_REF][START_REF] Zhang | Numerical study on deformations in a cracked viscoelastic body with the extended finite element method[END_REF] to address linear viscoelasticity problems with inclusions and cracks. Damage mechanics based approach can further be broadly categorised as Cohesive Zone Model (CZM) [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF] and Continuum Damage Mechanics (CDM) [START_REF]Introduction to Continuum Damage Mechanics[END_REF][START_REF] Lemaitre | How to use damage mechanics[END_REF]. Both approaches account for the micro-damage before fracture that can be associated to the nucleation, coalescence and growth of voids to form macro-cracks. Some works on use of CZM to model viscoelastic fracture can be found in [START_REF] Yoon | Damage dependent constitutive behavior and energy release rate for a cohesive zone in a thermoviscoelastic solid[END_REF][START_REF] Rahulkumar | Cohesive element modeling of viscoelastic fracture : application to peel testing of polymers[END_REF]. CDM is another popular approach for modeling damage in viscoelastic materials. Some of the commonly used CDM approaches are non-local integral damage models [START_REF] Pijaudier-Cabot | Non-local damage theory[END_REF][START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF], phase field damage models [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture : Variational principles and multi-field FE implementations[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF] and Thick Level Set(TLS) approach [START_REF] Moes | A level set based model for damage growth : The thick level set approach[END_REF]. For e.g. Shiferaw et al [START_REF] Shiferaw | Development of a damage viscoelastic model using the thick level set approach to fracture : 1D modeling and comparison to uniaxial tension stress tests on bituminous specimens[END_REF] used the TLS approach for modeling damage in viscoelastic materials for a wide range of temperature. Phase field models have also been used to model cracks in viscoelasticity [START_REF] Dammaß | A unified phase-field model of fracture in viscoelastic materials[END_REF][START_REF] Shen | Fracture of viscoelastic solids modeled with a modified phase field method[END_REF][START_REF] Yin | Viscoelastic phase-field fracture using the framework of representative crack elements[END_REF]. Another class of method that was recently developed under the formalism of CDM was the Lip-field approach [START_REF] Moes | Lipschitz regularization for softening material models : the Lipfield approach[END_REF][START_REF] Chevaugeon | Lipschitz regularization for fracture : the Lip-field approach[END_REF] to fracture. The main difference of this approach from the phase field approach lies in the fact that damage energy functional is not a function of the damage gradient.

Fracture behavior of viscoelastic materials (e.g. bituminous materials) is a complex phenomenon due to their highly rate-sensitive and temperature-dependent behavior. Modeling of such fracture behavior in viscoelasticity for a wide range of loading and temperature conditions can still be considered to be a relatively open area. This lays the objective of the present paper. Initially a viscoelastic model based on the formalism of generalized standard materials [START_REF] Germain | Continuum thermodynamics[END_REF][START_REF] Halphen | Sur les materiaux standards généralisés[END_REF] is proposed by defining the respective free energy and viscous dissipation potentials. In particular, the Generalized Kelvin-Voigt (GKV) model [START_REF] Brinson | Polymer Engineering Science and Viscoelasticity[END_REF] is used to represent the viscoelastic behaviour. The temperature effects are taken care of by the Time Temperature Superposition Principle (TTSP) [START_REF] Williams | The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass-forming Liquids[END_REF][START_REF] Brinson | Polymer Engineering Science and Viscoelasticity[END_REF]. Then an isotropic damage variable is being introduced into the model based on the formalism of CDM. Further, an incremental potential is proposed by gathering the respective free energy and dissipative potentials. This leads to the variational formulation of the problem under consideration, with the minimization of the incremental potential providing the solution of the mechanical problem. Finally, the mechanical problem is regularized by the introduction of a length scale using the Lip-field approach. We keep the discussion of the numerical aspects of the Lip-field approach to minimal and the details can be found in [START_REF] Moes | Lipschitz regularization for softening material models : the Lipfield approach[END_REF][START_REF] Chevaugeon | Lipschitz regularization for fracture : the Lip-field approach[END_REF].

The paper is organized as follows : The Lip-field regularized viscoelastic damage model is described in the next section. In Section 3, the numerical aspects are briefed. Numerical results are then presented in Section 4 for different loading rates and temperature. Finally, the paper concludes with summary and possible future works.

The mechanical model

Consider the deformation of a body initially occupying a domain Ω through a displacement field u. We assume small, quasi-static deformations under isothermal conditions. The Cauchy stress and small strain tensor is denoted as σ and ε.

ε(u) = 1 2 (∇u + ∇ T u) (1)
Regarding the boundary conditions, the displacement is applied on a part of the boundary Γ u . On the rest of boundary Γ N , zero traction forces are assumed. We also assume that there is no body force. To be kinematically admissible at any given time t, the displacement field must belong to U (t) :

U (t) = {u ∈ H 1 (Ω) : u = u d (t) on Γ u } (2)
The equilibrium condition then reads

Ω σ : ε(u * ) dΩ = 0, ∀ u * ∈ U * U * = {u ∈ H 1 (Ω) : u = 0 on Γ u } (3)
The equilibrium and the kinematic equations Eq. ( 2) and ( 3) must be complemented with the constitutive equations. For the constitutive behaviour, we consider the Generalized Kelvin Voigt (GKV) model [START_REF] Brinson | Polymer Engineering Science and Viscoelasticity[END_REF].

The schematic of GKV model in 1D case is shown in Fig. 1, where ε 1 , ε 2 , ..., ε n are the internal variables associated with the viscous strain in each Kelvin Voigt (KV) unit. E 0 , E 1 , .., E n and τ 1 , ..., τ n represent the spring constants and retardation times of the dashpots respectively. We assume that the Poisson's ratio ν is constant for all Kelvin-Voigt units (and time independent).

E 0 E 1 E 2 E n τ 1 τ 2 τ n ε 2 ε n ε 1 ε Figure 1 -Generalized Kelvin Voigt (GKV) model
To account for the effects of temperature, time-temperature superposition principle [START_REF] Brinson | Polymer Engineering Science and Viscoelasticity[END_REF] is used. This allows to use the Williams-Landel-Ferry (WLF) [START_REF] Williams | The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass-forming Liquids[END_REF] equation to find the shift factor α T as follows :

log(α T ) = -C1(T -T ref ) C2 + T -T ref (4)
where C1 and C2 are the material constants. The shift factor proportionally alters the retardation time to account for the changes in temperature as follows :

τ i (T ) = α T (T )τ i (T ref ).
We will use the formalism of Generalized Standard Materials [START_REF] Germain | Continuum thermodynamics[END_REF][START_REF] Halphen | Sur les materiaux standards généralisés[END_REF] to describe the constitutive behaviour. The model is characterised by the free energy potential ψ(ε(u), ε 1 , ..., ε n ) and the (viscous) dissipation potential ϕ( ε1 , .., εn ). We use the implicit time discretisation. Consider the displacement and internal variables (viscous strains) known at some time t m . Then finding the solution pair (u m+1 , ε 1,m+1 , ..., ε n,m+1 ) at the next time instant t m+1 = t m + ∆t amounts to the following minimization problem [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites : I. Incremental variational principles[END_REF].

(u, ε 1 , ..., ε n ) = arg min u ′ ∈Um ε 1 ′ ,..,εn ′ F (u ′ , ε 1 ′ , .., ε n ′ ; u m , ε 1,m , .., ε n,m , ∆t) (5)
where U m is the short hand notation for U (t m+1 ) and we dropped the m + 1 indices(for the rest of the paper) to simplify the notations. We will also drop the constants and known variables from previous time step in F whenever possible for clarity. The expression for F is follows [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites : I. Incremental variational principles[END_REF] :

F (u, ε 1 , .., ε n ) = Ω ψ(ε(u), ε 1 , ..., ε n ) + ∆t ϕ ε 1 -ε 1,m ∆t , .., ε n -ε n,m ∆t dΩ (6)
where the time derivative of the viscous strains

ε i (i = 1, 2.., n) is replaced by the differential quotient (ε i -ε i,m )/∆t.
The free energy and dissipation potential are given by

ψ = 1 2 ε - n i=1 ε i : E 0 : ε - n i=1 ε i + n i=1 1 2 ε i : E i : ε i (7) ϕ = n i=1 1 2 εi : τ i E i : εi (8)
where the linear elasticity tensor E i = 2µ i I + λ i 1 ⊗ 1 is defined in terms of Lame's constants λ i and µ i associated to i th Kelvin-Voigt unit , with I and 1 being the fourth and second order identity tensors respectively. The convexity of the free energy and viscous dissipation potential ensures that the viscous dissipation is positive, thus satisfying the second law of thermodynamics.

F (u, ε 1 , .., ε n , d) = Ω g(d)ψ(ε(u), ε 1 , ..., ε n ) + ∆t g(d)ϕ ε 1 -ε 1,m ∆t , .., ε n -ε n,m ∆t dΩ (9)
In our case, we consider the free energy of the springs as the driving force for the damage. So the evolution of the damage is given by the minimization of a different incremental potential F d [START_REF] Moes | Lipschitz regularization for softening material models : the Lipfield approach[END_REF][START_REF] Chevaugeon | Lipschitz regularization for fracture : the Lip-field approach[END_REF]]

F d (u, ε 1 , .., ε n , d) = Ω g(d)ψ(ε(u), ε 1 , ..., ε n ) + Y c h(d) dΩ (10)
where h(d) is a softening function that controls the damage evolution behaviour. As will be seen, both g(d) and h(d) will take user defined convex function so that F d preserves convexity w.r.t damage variable d. Therefore, the solution of the mechanical problem (u, ε 1 , .., ε n , d) at time t m+1 can be formulated as the minimization of the two different incremental potentials Eq.( 9), [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] as follows :

(u, ε 1 , ..., ε n ) = arg min u ′ ∈Um ε 1 ′ ,..,εn ′ F (u ′ , ε 1 ′ , .., ε n ′ , d) (11)
d = arg min d ′ ∈An∩L Ω,lc F d (u, ε 1 , .., ε n , d) (12)
The damage irreversibility constraint is given by the space A m .

A m = {d ∈ L ∞ (Ω) : d m ≤ d ≤ 1} (13)
The damage irreversibility constraint ensures that the dissipation due to damage is positive. The damage softening models exhibit spurious mesh dependency problems [START_REF] Jirásek | Objective modeling of strain localization[END_REF]. Hence in our case a length scale l c is introduced into the model through the Lipschitz regularized space L Ω,lc as explained below : The Lipschitz constant associated to any given damage field d over the domain Ω is the minimum M such that the following holds

|d(x) -d(y)| ≤ M dist(x, y) ∀x, y ∈ Ω (14)
where dist(x, y) is the minimal length of the path inside Ω joining x and y (the distance is considered infinite if the two points cannot be connected inside Ω). The value M defined above is denoted lip(d).

This leads to the following definition of the Lipschitz regularized space [START_REF] Moes | Lipschitz regularization for softening material models : the Lipfield approach[END_REF] L

Ω,lc = d ∈ L ∞ (Ω) : lip(d) ≤ 1 l c (15)
We consider the following energy degradation function and softening function for the rest of the study.

g(d) = (1 -d) 2 + α(1 -d)d 3 (16) h(d) = 2d + 3d 2 (17)

Numerical aspects

The mechanical problem under considered is discretized by means of finite elements in space and finite difference in time. The incremental potential F (Eq. ( 9)) is convex w.r.t u, ε 1 , .., ε n for a fixed damage field d and the incremental potential F d (Eq. ( 10)) is convex w.r.t to damage field d for fixed u, ε 1 , .., ε n . Hence, a staggered scheme is used in which the minimization (given by Eq. ( 11) and ( 12)) is performed alternatively until convergence. It is to be noted that while performing the minimization of F (Eq. ( 11)), the damage field is kept fixed (and equal to its latest available value) and while performing the minimization of F d (Eq. ( 12)), the deformation u and viscous strains ε 1 , ..., ε n are fixed (and equal to the latest available value).To accelerate the minimization process for finding the damage field, we employ a technique detailed in [START_REF] Moes | Lipschitz regularization for softening material models : the Lipfield approach[END_REF][START_REF] Chevaugeon | Lipschitz regularization for fracture : the Lip-field approach[END_REF]. Briefly, it contains three steps :

1. First a purely local minimization is performed without the Lipschitz constraint to find the trial damage d as follows :

d = arg min d ′ ∈Am F d (u, ε 1 , .., ε n , d) (18)
2. Then the set of active nodes is found. Active nodes are defined as the nodes where the d / ∈ L Ω,lc . d is assigned d in all the nodes except for the active nodes.

3. Finally, a global minimization for the damage variable given by Eq. ( 12) is performed only in the subdomain occupied by active nodes.

Step 2 is done efficiently using a fast marching technique described in [START_REF] Chevaugeon | Lipschitz regularization for fracture : the Lip-field approach[END_REF]. For the Step 3, cvxopt [START_REF] Andersen | The cvxopt linear and quadratic cone program solvers[END_REF] library of the Python has been used to impose the damage irreversibility (A m ) and Lipschitz constraint (L Ω,lc ). This library allows to impose the damage constraints in discrete form through first and second order cone constraints.

Simulation results

A TDCB (Tapered Double Cantilever Beam) geometry shown in Fig. 2 is used for the study. This geometry has already been used in [START_REF] Chevaugeon | Lipschitz regularization for fracture : the Lip-field approach[END_REF] for studying fracture in elasticity using Lip-field approach. This geometry ensures that the crack growth is relatively stable as the imposed displacement increases.

i,ref (s) N/A 1E-5 1E-4 1E-3 5E-3 1E-2 1E-1 1 10 1E2 5E2 1E3 Table 2 -Parameters for GKV model Y c (N/m 2 ) ν α C 1 C 2 ( • C) 30 .2 .1 29.095 211.60 Table 3 -Material parameters
The dimension for the TDCB geometry is listed in Table 1. The loading in y-direction is applied through linearized rigid body motion on the boundary two holes. We make the assumption of plane strain and isothermal conditions. For the GKV model (viscoelastic model) , the parameters are listed in Table 2 and these values were extracted from [START_REF] Shiferaw | Development of a damage viscoelastic model using the thick level set approach to fracture : 1D modeling and comparison to uniaxial tension stress tests on bituminous specimens[END_REF].The retardation times of the GKV model are listed for the reference temperature T = 20 • C. The other material parameters of the model are listed in Table 3, where ν is the Poisson's ratio. It is to be noted that for the present numerical study Y c is considered constant irrespective of the loading rate and temperature. However this is not the case as observed in experiments. For the assumed plane strain condition, the Lame's constants are found as follows :

λ i = E i ν/(1 + ν)/(1 -2ν) and µ i = E i /2/(1 + ν).
The TDCB geometry is spatially discretized through linear triangular elements (number of nodes = 8112 and number of triangular elements = 15944). We use l c = 5 mm to ensure that there are enough elements (10 elements) in the crack band. The time step dt is set to 1e -4 s for all the simulations. The crack length is found approximately as Ω d/l c dΩ.

Influence of loading rate

Viscoelastic models are rate dependent. We study the influence of the imposed displacement rate on the TDCB specimen to study the rate sensitivity of the considered model. The simulation was performed at T = T ref = 20 • C (shift factor α T = 1). We consider the following three imposed loading rates for the study : 1 mm/s, 0.1 mm/s, 0.01 mm/s. Figure 3 shows the force-displacement curves and crack length vs displacement curve for the considered loading rates. It can be seen from the Figure 3a as the imposed displacement rate increases, the load carrying capacity of the specimen increases and the displacement at rupture decreases, which is a typical response of the viscoelastic materials undergoing fracture. Moreover, Figure 3b shows that the crack grows faster as the imposed displacement rate increases. This can be explained as follows : a) at extreme loading rates (very low and very high loading rates) the viscous dissipation is negligible and most of the dissipation happens through fracture. b) Higher the loading rate higher the stiffness of the model (glassy region) compared to low loading rates (rubbery region). Both in glassy and rubbery regions the model behavior is similar to behavior in elasticity with negligible viscous dissipation, except for the transition region where the viscous dissipation is significant. Hence at extreme loading rates, the model behavior is similiar to elasticity but with different stiffness and hence differenct crack speeds. Snapshots for the stress in y-direction and damage field (crack band) for the imposed displacement rate of 1 mm/s at u imposed = .02 mm is shown in Figure 4. It can be seen that the fracture results in mode-I crack pattern with stress being maximum close to the crack tip. It can be seen from Figure 4a that the crack initially grows faster and then it slows down after approximately three quarters of the maximum size of crack.

Pa

Influence of temperature

We now study the influence of temperature on the fracture behaviour. We consider the following three different temperatures for the study : Figure 5 shows the force displacement curves and crack-length vs imposed displacement for different temperatures for fixed imposed displacement rate of .1 mm/s. It can be seen that as the temperature increases the load carrying capacity of the model decreases and the rate of increase of crack length decreases. Moreover, it can be observed that increase in temperature has similar (qualitative) effect as decrease in loading rate, which is typical to viscoelastic materials.

Conclusion

The paper described a viscoelastic damage model being implemented within the variational framework and regularized using the Lip-field approach. The numerical results produced for various loading rates and temperature conditions exhibit the similar qualitative behavior observed in viscoelastic fracture.The novelty of the present work lies in the use of Lip-field approach for the first time to model damage in viscoelastic materials. This can be considered as a first step to model damage softening in viscoelastic materials.

The future work will involve implementation of a viscoelastic damage model with unilateral effects within the Lip-field approach and to calibrate the model with experimental results. This will help to use the model for some real world applications.

  An isotropic (internal) damage variable d is introduced into the GKV model within the framework of Continuum Damage Mechanics (CDM) by use of the effective stress principle. d = 0 indicates virgin material and d = 1 indicates crack. In particular, an energetic degradation function g(d) diminishes the materials capability to store the free energy and dissipate viscous energy to yield the modified incremental potential as follows (with g(d) : [0, 1] → [1, 0]) :

Figure 2

 2 Figure 2 -TDCB Geometry specimen

Figure 3 -

 3 Figure 3 -Force-displacement plot and crack length -displacement plot for different loading rates

Figure 4 -

 4 Figure 4 -From left to right a) Force-displacement and crack length vs displacement plot b) stress field in y-direction c) damage field (all contour plots showed for u = .02 mm)

  10 • C, 20 • C and 30 • C. The difference in temperature is accounted for in the model by shifting the retardation times τ i by the shift factor α T calculated using the WLF equation (Eq. (4)). Crack length vs displacement

Figure 5 -

 5 Figure 5 -Force-and crack length-displacement plots for different temperature at uimposed = .1 mm/s.

Table 1 -

 1 Dimensions for the TDCB geometry

	KV unit (i) 0	1	2	3	4	5	6	7	8	9	10	11
	E i (M P a) 31770 87398 123414 65830 62457 62661 7305 12500 418	1743 79	39
	τ