G Jodin
email: jodin@laplace.univ-tlse.fr

Y Bmegaptche Tekap

J M Saucray

J F Rouchon
email: rouchon@laplace.univ-tlse.fr

M Triantafyllou

M Braza

Optimized design of real-scale A320 morphing high-lift flap with shape memory alloys and innovative skin

Keywords: morphing, morphing skin, Shape Memory Alloys, optimal design, industrial aircraft

This article proposes a design approach based on multi-criteria optimization. Applied to the conceptual design of a cambered control flap of Airbus A320, different skin technologies and actuator topologies are compared. Amongst the possible shape memory alloy actuators, an agonist-antagonist solution is found the most suitable for the application. Elastic skins are compared to a bioinspired innovative skin and feather concept. The results are based on simplified but realistic aircraft specifications that consider environment and industrial concerns. Finally, the optimization results in a feasible low weight, low power consumption morphing wing design. This design is a basis for detailed incoming integrated aircraft system.

Introduction

Limiting energy consumption has become a central concern for reducing the aircraft operational costs. Therefore, improving aerodynamic performance is a way to reduce the fuel consumption during flight.

Current airfoil shapes are generally optimized for one working point, corresponding to nominal cruise conditions. During flight, the altitude, the weight and the speed are continuously changing. Hence this design is suboptimal for the whole aircraft mission. Changing the shape of the wing during a mission can save several percents of fuel for a regional passenger aircraft [START_REF] Lyu | Aerodynamic shape optimization of an adaptive morphing trailing edge wing[END_REF]. The concept of real time shape adaptation enabling multipoint optimization is called morphing. Wing morphing is of great interest to increase aerodynamic performance and to decrease noise throughout the different flight steps. Morphing requires a structure flexible enough to be easily deformed whilst being stiff enough to withstand the aerodynamic loads. This paradoxical compromise leads to current issues in skin and actuator design.

Within the framework of aircraft aerodynamic performance, morphing has been known for decades [START_REF] Barbarino | A review of morphing aircraft[END_REF]. It was shown by Lyu et al. [START_REF] Lyu | Aerodynamic shape optimization of an adaptive morphing trailing edge wing[END_REF] that camber control of the trailing edge of a wing is very efficient to improve airliner performance. Their work also exhibits that morphing applied to a limited part of the airfoil chord may feature effectiveness comparable to that of entirely morphed airfoils. Additionally, the required deformations result in a camber change by the order of 7% of the high-lift flap's chord.

Despite the potential gains, research concerning civil aircraft design has not led to commercial applications. One can explain this as technology is not mature enough: reliability, maintenance, mass and power consumption of the added devices deteriorate the assessment. The skin or the interface between the structure and the outside airflow is the current technological bottle neck. Studies simply do not deal with this issue (experiments at low velocity do not exhibit issue with elastomeric skins), others propose a degraded skin (i.e. corrugated skin for example that cause turbulent transition). As the skin is a part of the morphing system, it has to be taken into account. Thill et al. [START_REF] Thill | Morphing skins[END_REF] draw a review of morphing skin and its challenges. With this regard, relatively high Technology Readiness Level (TRL) projects, targeting current industrial airliners at true scale, have been undertaken. The European research program SARISTU 6 and the European research program CleanSky7 focus on operating cost reductions as well as on improving the aerodynamic performance. Work packages of these two programs focus on adaptive morphing trailing edge, [START_REF] Dimino | Smart Intelligent Aircraft Structures (SARISTU): Proceedings of the Final Project Conference[END_REF][START_REF] Pecora | Toward the bi-modal camber morphing of large aircraft wing flaps: the cleansky experience[END_REF]. This actuation provided by servomotors and hinges allows for camber control during the different flight phases/assets. Another trailing edge morphing concept, called Adaptive Compliant Trailing Edge, was developed by NASA in cooperation with FlexSys Inc. 8 . This concept features an adjustable structure which can be actively deformed. Endurance flight tests for this concept were performed and described by [START_REF] Kota | Flight testing of the flexfloil tm adaptive compliant trailing edge[END_REF].

However these new adaptive structures are actuated through conventional A c c e p t e d M a n u s c r i p t actuators like electromechanical or hydraulic servomotors. Recent advances in the field of smart materials show the potential to overcome difficulties to make a wing both stiff enough to withstand the loads, and flexible enough to be easily deformed [START_REF] Barbarino | A review on shape memory alloys with applications to morphing aircraft[END_REF]. The related research focuses mainly on low TRL, mainly applied to low scale Micro Air Vehicles. Among electroactive materials, Shape Memory Alloys (SMAs) are frequently used. SMAs are characterized by thermomechanical behaviors, and most applications use an electrical resistor or the resistance of the SMAs themselves to activate the transformation. Different morphing concepts were developed, an overview of which is presented by Barbarino et al. in [START_REF] Barbarino | A review of morphing aircraft[END_REF]. A counter example of the low TRL research is a recent industrial application of the SMAs within the Boeing CLEEN research program. A topic of the program consists in a flap actuated by a SMA twist tube, [START_REF] Calkins | Flight test of a shape memory alloy actuated adaptive trailing edge flap[END_REF]. Safety issues have been solved using a redundant hydraulic actuator. It also provides brake and damping functions to maintain a deflected flap without requiring actuation energy. It shows that certification issues due to the use of SMA can be addressed, as a flight test campaign has been successfully done, reaching TRL 7. A recent article [START_REF] Sun | Morphing aircraft based on smart materials and structures: A state-of-the-art review[END_REF] provides a wide overview of the current state of the art regarding aircraft morphing using smart materials. One of the conclusion recommendations is that structures with smart materials should be optimized to allow future commercial aircraft applications. This is a aim of the present study.

In the works previously mentioned the desired operation is generally accomplished using one specific smart material. The proposed approach reflects an airflow point of view. Large flow instabilities can be manipulated through camber control whereas smaller instabilities can be modified through active turbulence manipulation of the trailing edge wake's shear layers. For ten years of collaborative effort from two laboratories (IMFT and LAPLACE) this approach has been leaded to the electro-active hybrid morphing concept: a low frequency (< 1 Hz) camber control (∼10% of the chord) thanks to SMAs and a higher frequency vibrating trailing edge (fractions of millimeters up to 400 Hz). It was demonstrated that the flow dynamics are significantly affected by the trailing edge actuation. The wing's wake energy was reduced, leading to an improvement in aerodynamic performance, according to Sheller et al. in [START_REF] Scheller | Trailing-edge dynamics of a morphing {NACA0012} aileron at high reynolds number by high-speed {PIV}[END_REF]. Also the hybrid morphing has significant effects on lift, according to Jodin et al. in [START_REF] Jodin | Dynamics of a hybrid morphing wing with active open loop vibrating trailing edge by time-resolved piv and force measures[END_REF]. Following reduced scale studies, the collaboration now aims true scale applications.

This article deals with a design approach based on multi-criteria optimization. Applied to the preliminary design of a cambered control flap of Airbus A320, the first section describes the morphing flap concept with design requirements. Then, shape memory alloy actuator topologies are examined. Thirdly, an automatic design routine is used in the optimization algorithm. A fourth section compares two different skin technologies based on the optimized morphing flaps. This multi-objective optimization algorithm allows for the selection of the final design in a last section before conclusion. It is a basis for detailed incoming integrated aircraft system.

The design of the higher frequency vibrating trailing edge as well as precise design and certification of the resulting industrial system is out of the scope of the study.

The reader may note that the approach focuses on conceptual system design which is under-constrained regarding true industrial systems. This simplification allows for the use of new computational techniques to design a system considering multiple components with wide possibilities. This point of view on adaptive structure is the major contribution of this article.

Morphing wing concept and optimized design

Specifications

The proposed morphing concept is applied to a high-lift flap, but the assumptions are valid for a whole morphing wing. The flap profile has been adapted from Airbus specifications. The flap is 1 m chord and 2 m span. A chordwise loading is specified, equivalent to 1.5 tons of aerodynamic upward forces. The force distribution is illustrated on Figure 2 by a sketch describing the morphing concept.

The function assessed by the proposed morphing is to adapt the wing shape configuration. This allows to change the shape that correspond to the lowest drag for every flight step. The high lift function is not addressed so that safety is not critical for this function. A shape envelop is defined, the morphing flap must be able to smoothly vary between the two shapes presented in Figure 1. These shapes correspond to up or down displacement of the trailing edge by 10 cm i.e. 10% of the flap chord.

Articulated concept

The morphing flap is based on an articulated ribs where Shape Memory Alloy actuators control the rotations of the elements around the hinges. The proposed concept, presented in Figure 2 • hinges allow the rotation of the articulated ribs. Parts of forces are transmitted through these components without generating much parasitic force (or torque) when rotated.

• actuators are devices that transmit mechanical energy to the structure. The actuators are responsible for the shape control and have to counteract aerodynamic forces as well as internal forces coming from the other components.

• skins or covering devices guarantee the airtightness of the wing, transmit the aerodynamic forces to the structure and ensure a smooth shape during morphing. The skin must endure deformation without unexpected displacements like bumps or wrinkles.

Additionally, mechanical stops are provided to limit the rotations of the articulations, thus preventing overloads in actuators. The internal structures represented by the articulated ribs actually consist in an engineered mechanical structure composed of ribs and spars. The fine design of this structure is outside the scope of this study, it is accepted that the lower the force in the structure, the lighter the structure. This reasonable assumption could be justified by the aim of matter inside the wing. Excluding the skin and the actuators, the inner matter is only dedicated to hold the components and transmit the forces from the skin to the spar, through the actuators and hinges. This means that we neglect dynamic and thermal inertial effects of the inner structure. Considering a given structure topology (e.g. ribs, spar stiffeners), and because strength of material relies on the sections and volume of matter to withstand forces

Vertical force distribution

Leading edge Trailing edge A c c e p t e d M a n u s c r i p t (e.g. maximal stress), lowering the inner structure forces decreases the required matter so the weight. From an industrial point of view, the torsion resistance along the span direction is critical issue for commercial aircraft design. By the way, torsion resistance for flap is not a crucial function; the torsional stiffness of the wing comes from the main wing box, the flap is neglected. The flap's main solicitations come from the resistance to spanwise bending moment. For the developed application, the articulations axes are spanwise, therefore the spanwise bending stiffness is not compromise. Regarding the torsion resistance, the consider flap is cut in multiple wing boxes linked together by the actuators and the hinges. All the torsion forces are transmitted through the actuators, and these forces are parts of the force specifications.

Articulated ribs

Gliding bearings (plain bearings with steel-TEFLON contacts) have been selected for the hinge function. They are suitable for low rotation velocities. They are compact, lightweight and generate low friction torque.

The proposed actuators consist in cylinder like actuators. Composed of shape memory alloy wires, they are able to pull on the articulated ribs, thus imposing the rotations. More detailed are presented in section 3.

Different skin technologies can potentially address the previous specifications. A comparison between two technologies is presented in section 5, based on optimization.

The original purpose of this study is the comparison of the different technologies that can be applied.

Actuator technology

The proposed morphing concept relies on actuators using Shape Memory Alloy (SMA). This section first shows that SMAs are suitable for the actuators. Secondly a design model of such actuators is proposed. Thirdly, actuation topologies are compared.

Non linearity effects like material hysteresis, loose SMA wires or dependencies between actuators are not discussed in this article. These specific issues can be taken into account by suitable control laws that will be addressed in a future study. However, a few elements regarding the control are discussed at the end of the section 3.4.

Shape Memory Alloy behavior

The Shape Memory Alloys (SMAs) are used liked artificial muscles embedded on both sides of a hinge. This antagonist topology allows for the actuation by contracting one of the SMA. Amongst smart materials, SMA are metallic alloys that exhibit an impressive thermo-mechanical coupling, due to crystallographic phase changes at microscopic scale. They have been studied for decades, Lexcellent draws a complete handbook about SMAs, [START_REF] Lexcellent | Shape-memory alloys handbook[END_REF]. The most common are based on Nickel-Titanium alloys with small amount of chemical additives to tailor thermo-mechanical properties. Basically, cold SMA is martensitic which exhibits low stiffness and pseudo-plastic behavior. Hot SMA is austenitic, characterized by a higher stiffness and super-elastic behavior. This corresponds to a "two-material" behavior characterized by strain levels up to 7% and stress level up to 600 MPa, as presented by hot and cold stress-strain characteristics in Figure 3. The SMA's high specific actuation energy (about 1 kJ/kg), stainless property and growing maturity make this material of interest for aeronautic applications. SMA actuators are not designed in respect of the maximum stress and strain, but in respect of fatigue. A cycle life up to one million actuation cycles have been shown on small diameter SMA wires at limited stress and strain, [START_REF] Jani | A review of shape memory alloy research, applications and opportunities[END_REF]. Generally, the lower the diameter and strain-stress loading, the longer the cycle life. Areas drawn on Figure 3 indicate the expected maximum number of cycles. Therefore, a reliable design relies on a specification area where the working points are included. To ensure about 100,000 cycles with 1.5 mm diameter for the first demonstrator, limits are set at σ SM A max = 150 MPa stress and SM A max = 2% strain. This choice allows a wire section large enough limiting the number of wires to be integrated, simplifying the prototype assembly.

Actuator sizing 3.2.1. Specifications

As the morphing displacement is slow, quasi-static conditions are assumed. Thus the Newton's second law can be applied on a wing section comprising a rib section and all the following ones until the trailing edge. The resulting forces and moments expressed at the hinge location are: the actuator forces F acti and moments M acti , the hinge friction moment M hinge and reaction force F hinge , the 7 A c c e p t e d M a n u s c r i p t resulting aerodynamic force F aero and moment M aero applied on all the wing elements following the hinge.

The force balance is presented in equation 1.

F hinge = -(F aero + ΣF acti ) M act = -(M hinge + M aero ) (1)
Two effects sum up this balance: 1-the force transmitted to the previous ribs F hinge allows the sizing of the hinge; 2-the resulting moment specifies the sizing of the actuator moment M act , where the actuators' forces and lever arms are design variables.

Design algorithm

The SMA actuators are sized according to the specified moment M act . The design algorithm used in this study attempt to size the actuator with minimal mass and force, whilst respecting the available space and material limits. The algorithm outputs the actuator lever arm and length. The first approach is to use the maximum available length to maximize the lever arm, thus minimizing the actuator force for the given output torque. If the lever arm exceeds these limits, the length is decreased and the lever arm is set to the maximum. The required actuator force is evaluated, then the SMA section is calculated. This provides the amount of SMA mass (number and size of the wires) as well as the actuator mass and the energy consumption. Penalties are calculated if material or geometry limits are exceeded.

Actuation topology: counter spring versus agonist-antagonist concepts

This subsection deals with the actuation topology. As the SMA actuators can only pull, a device has to be implemented to recover the initial shape. Since the aerodynamic forces will not always act in the direction of desired actuation, they cannot be relied upon exclusively for recovery from unidirectional actuation. Two return mechanisms are compared: a passive counter-spring and an antagonistic SMA actuator. The representative comparison cases are presented in Figure 4.

Analytical models are used to compare the solutions. It is assumed that the forces applied on the articulated ribs are the actuator force F act and the counter-spring force F sp or respectively the antagonistic actuator force F ant . Force moments are calculated at the hinge's center. The lever arms are noted h 1 and h 2 . The usable actuator output torque is M u . The SMA actuator minimum force is F act cold and maximum force is F act hot . The counter spring is assumed as a linear spring, with its stiffness K sp and its relaxed length L sp0 .

Counter spring solution model

The counter spring aims to counter act the cold SMA force from the actuator to recover the initial shapeas presented in Figure 4a. The counter spring prestrain must apply a minimum moment equal to the opposite of the cold SMA moment. In order to limit the exceed in this force, the stiffness K sp have to 8 A c c e p t e d M a n u s c r i p t be minimum, this is possible by selecting a low stiffness spring with a high prestrain. For integration reason, the pre-strain is reasonably limited to 70% of the actuator length L SM A (δL = 70%L SM A ). The torque balance gives the expression of the output moment in equation 2. This equation also contains the resulting force F X in the hinge, where δL is the spring pre-strain. It appears that the output torque of the spring solution depends on the rotation angle.

L SMA h 2 1 2 M u h 1 θ K act L SMA h 2 1 2 h 1 θ act ant M u (a) Counter spring concept. L SMA h 2 1 2 M u h 1 θ K act L SMA h 2 1 2 h 1 θ act ant M u (b) Agonist-antagonist concept.
Maximum actuation:

M u = F act hot • h 1 -F act cold • h 1 • (1 + h1 • θ δL )
Minimum actuation:

M u = -F act cold • h 1 • h1 • θ δL F X = F act + F act cold h 1 h 2 (1 + h 1 • θ δL ) (2)

Antagonist solution model

Regarding the agonist-antagonist concept, the representation on Figure 4b distinguishes the actuator F act and the antagonistic actuator F ant which is design to be actuated to counter act the cold actuator. The torque balance as well as the resulting hinge force are presented in equation 3. Because of assuming that cold SMA exerts constant force, the output torque does not depend on the rotation angle.

Maximum actuation: M

u = F act hot • h 1 -F ant cold • h 2 Minimum actuation: M u = F act cold • h 1 -F ant hot • h 2 F X = F act + F ant (3)

Topology comparison

The sizing algorithm has been applied to an articulated hinge for a range of output torques. The articulations' parasitic torque is taken into account. The performance of the technologies are presented in Figure 5. The hinge force and the amount of SMA are drawn depending on the specified output torque. For information, a non-functional system without antagonistic actuator nor spring is also presented by dash lines on the figure. This "no-return system" can only be actuated in one way, thus it is not suitable for morphing applications. This system is not valid for comparison, but it allows the reader to understand that the addition of an extra antagonistic actuator does not increase the inner structure forces by a large factor. It also provides an idea of the performance limit of the technology, this system respects the one way actuation requirements, so two way actuation systems that require more devices are heavier. The conclusion is significant: counter spring solution adds forces that requires larger total SMA mass and larger hinges.

Finally, the agonist-antagonist actuator concept is more efficient than a counter spring. The solution is lighter, requires a lower SMA mass with lower forces in the structure, even if a more elaborate control system would be designed. From the integration point of view, the counter spring can be materialized by an elastic skin, as presented in the next Section 5. The retained concept is a bio-inspired analogy regarding animal's articulations where a simplified version of an elbow is actuated by two antagonistic muscles (biceps and triceps). This "naturally selected solution" is preferred than one large muscle with a counter spring materialized by tendons. actuator topology is selected for the morphing wing design.

Control of the actuated hinges

The following paragraph deals with the control strategy of the SMA actuated antagonistic hinges. Accurate control of SMA is well understood in the literature [START_REF] Ma | Position control of shape memory alloy actuators with internal electrical resistance feedback using neural networks[END_REF][START_REF] Kha | Position control of shape memory alloy actuators by using self tuning fuzzy pid controller[END_REF], including antagonistic configurations [START_REF] Gorbet | A novel differential shape memory alloy actuator for position control[END_REF]. Example of antagonistic actuator configurations for morphing wings with bending beam deformation mechanisms is developed in [START_REF] Jodin | Implementation of a hybrid electro-active actuated morphing wing in wind tunnel[END_REF][START_REF] Jodin | On the multidisciplinary control and sensing of a smart hybrid morphing wing[END_REF].

Every actuator is equipped with a force sensor a temperature sensors. Figure 6 illustrates the controller's structure. The control loops are nested. The outside loop is a position loop for each hinge. The output of the PID controller is a reference torque to be applied by the actuators. Depending on the sign of the torque, the antagonist or the agonist actuator is piloted whereas the other is regulated to a given constant force to ensure SMA wire are taut. These forces are regulated by inner force control loops using bang-bang controllers. Therefore the hinge stiffness as well as the position are controlled. Such a topology is close to impedance control used in robotic systems. Temperature sensors are present for safety and monitoring.

Design algorithm for morphing flap

Introduction of the algorithm's aims

In order to allow comparison and optimization of flap designs, an algorithm has to design autonomously the morphing flap from input parameters. This algorithm -presented in the following -requires: the previously detailed actuator technology, the skin technology and the hinges' positions. The actuator design algorithm and skin design algorithm are discussed in sections 3.2.2 and 5.

Once the number of articulations and the technologies for the actuators and skin are selected, the only optimization variables of the global algorithm -i.e. the parameters sufficient to define a flap -are just the x and y positions of the hinges. Then, the locations, the size and the pre-strain of the actuators, as well as the skin and the hinges are automatically set by inner deterministic design loops. To guarantee the controllability of each articulation angle, every hinge is actuated by a dedicated actuator. After optimizations of flaps for different number of hinges, a compromise between complexity and accuracy in shape control has been found for five actuated hinges.

As the internal structure is not precisely modeled in the sizing algorithm, the number of actuators spread along the span direction of the wing is not an important matter. A realistic compromise between force distribution and complexity is a system of four series of actuators along the span. The optimization is then based on four similar wing-boxes of 50 cm each, forming the 2 m span flap.

Algorithm working principle

The design algorithm is programmed to be a cost function of a global optimization problem. The design algorithm is illustrated in Figure 7. The first step of the algorithm requires the hinge positions as input parameters. Starting from a hinge distribution and sizing rules (i.e. SMA limits and geometrical parameters), the algorithm firstly computes the optimized rotation angles that best fit the specified shapes for both low and high cambered airfoils. A gradient-based local optimization method has run separately on both two objective shapes, in order to find the extreme rotation angles. Secondly, the forces induced by aerodynamic loads (i.e. moments at hinge positions) are evaluated. Afterwards, the forces due to the considered skin technology are evaluated, based on the aerodynamic pressure and the rotation angles. Thirdly, actuators and hinges are sized. During the first sizing iteration the unknown forces are assumed zero. Consequently, the first actuators' designs 12 Page 12 of 25 AUTHOR SUBMITTED MANUSCRIPT -SMS-106292.R2

A c c e p t e d M a n u s c r i p t

do not take into account the cold antagonistic actuators nor the hinge parasitic torque. Next the antagonistic actuators are computed based on the cold actuator forces. Finally the hinges are sized, based on the forces of all the actuators and antagonistic actuators. This first loop is followed by other sizing loops. Several iterations are done until the modifications in the hinge parasitic torque and the force from the cold antagonistic actuators do not change by more than 5%. This percentage is in accordance with he design requirements and provides fast computations. If the design does not converge, the sizing algorithm stops with a penalty cost. As a fourth step, a sized system is obtained where the actuators, the antagonistic actuators and the hinges can work together. The sizing algorithm finally computes weights, power consumptions, shape approximation error, lifespan and penalization:

• The total weight is the sum of the agonistic and the antagonistic actuators, the hinges and an assumed simplified rib structure.

• The power consumption comes from the SMA activation, linked to the SMA amount and a thermal modeling of the actuators. The maximum power is computed when heating all the agonist actuators. The maintaining power corresponds to the average power needed to maintain all the agonist actuators at maximum temperature.

• The shape approximation errors are assessed with the evaluation of S AV , the average absolute error and S max , the maximum absolute error.

• The lifespan concerns the actuators. The system expected lifespan is limited by the actuator with the lowest lifespan. This is estimated from maximum stress and strain levels of the sized actuator, as shown in Figure 3.

• Penalizations can be added, when sizing difficulties are encountered. Penalties are calculated when dimensions of dedicated actuator place are lower than the minimum lever arm, when there are not enough places for hinges and when the design is not converged.

Skin technologies

The skin is the device that covers the whole wing and defines the interface with the airflow. Two different technologies are compared for the articulated flap. The first proposal consists in an elastic skin covering the flap composed by parts glued on the surface and unfixed ones to allow flap deformations 5.1. The second proposal is a promising innovative bio-inspired feather-like skin, investigated in section 5.2. For both technologies, the specifications are to withstand the aerodynamic pressure with less than 1 mm deformation.

Elastic skin 5.1.1. Working principle

The elastic skin concept consists in a taut membrane that is hung on the articulated ribs. The skin is free between the ribs at articulations levels, as presented in Figure 8. The free skin elements can be elongated when the articulations moves while the dynamic pressure tends to deform it.

An important issue of flexible skins are the possibility of wrinkling and bulging due to pressure and flap motion. In the following, the effect of the pressure is considered for both extreme flap deformations. The minimal length case with maximal aeroload is the worse case. A minimal tension force is determine to prevent both wrinkling and bulging. The skin is always taut to prevent wrinkling and the minimal tension force limits the pressure induced bulging to the 1 mm specification. By the determination of the length and tension force of the skin, this problem is addressed by design and does not require additional system. Mutlifunctional materials used as active skin with adaptive stiffness could potentially be elements to overdrive the current issues in flexible skin for commercial morphing aircrafts. For example, [START_REF] Sun | Mechanical properties of shape memory polymer composites enhanced by elastic fibers and their application in variable stiffness morphing skins[END_REF] is a study on the design, making and characterization of an electroactive composite material based on Shape Memory Polymer.

Model for optimal design

To model the skin, one free element is focused on the detailed view in Figure 8. The free length is note L, T is the skin tension force (considering a span length Span) and H is the normal deformation due the pressure P (supposed constant on the free length). The skin is modeled as a membrane; but reduced to two dimensions, the skin is assumed to be modeled like a taut cable. The relation between the normal deformation and the other parameters is depicted in Equation 4. In a case corresponding to constant length L, a choice of the couple length, tension (L, T ) is made to ensure a small normal deformation H < H max = 1 mm.

Due to the movements of the articulations, the skin length changes. Additionally the elasticity of the skin cause changes in the tension. Therefore,

14

A c c e p t e d M a n u s c r i p t

Equation 5 models the skin as a pre-stressed elastic material; where t is the skin thickness, E its Young modulus and δL the skin extension due to articulation movements. The skin pre-stress T pre corresponds to the minimum tension in skin when morphing movement impose minimum skin length. This tension is calculated from Equation 4to ensure a small normal skin deformation.

The tension T changes with the morphing movements. This tension has to be withstand by the structure and the actuators, so that this force has to be minimized for an efficient overall design. It can be shown that the choice of an optimal free length L opt minimizes the maximum skin tension T max . The analytic expressions of L opt and T max are presented in Equation 6 and 7.

H = P • Span • L 2 8 • T (4)
T = T pre + E • t • δL L (5)
L opt = 3 4E • t • δL • H max P (6)
T max = P • Span 8H max L 2 opt + E • t • Span δL L opt (7)
Finally, these equations are integrated in the design algorithm. This routine takes the skin displacement and the available room as inputs. It computes the free skin lengths and pre-stresses that fit the available room -optimum length are selected if possible. Outputs are the tensions that are transmitted by the hinges, and the maximum torques that actuators have to counter act, due to skin elastic forces when stretched during morphing.

Innovative feather and skin concept

The proposed innovative concept is bio-inspired from birds. Birds' bodies are covered by heterogeneous feathers. When a bird moves its wings, the wing profile stays smooth, as the feathers glide on each others. The feathers follow the aerodynamic shape. Birds have skin under the feathers. Inspired from this twofold structure, the proposed concept consists of internal airtight skins and external gliding devices. The internal skin does not have specific shape constraints but only allowing the motion and being airtight to prevent internal airflow going outside. It has to carry the average aerodynamic loads. The external slipping devices, like artificial feathers, slide on the airfoil shape and are thin enough to respect the 1 mm shape accuracy. Figure 9 presents the concept applied to one articulation.

The sizing of the internal skin does not require it to be extended, therefore the forces due to the internal skin are neglected. Particular attention is paid to ensure enough available space between the feathers and the actuators. The sizing of the feathers is decomposed threefold: geometric constraints, friction induced parasitic torque , pressure induced deformation. Geometric constraints concern the feathers' thicknesses that must be under 1 mm. Also, the feathers' lengths L must be larger than the gliding length to ensure the airfoil profile is always closed.

The friction induced torque comes from the calculation of the force exerted between the feather and the rib. Due to the rotation, the feather tip displacement y is evaluated. Assuming the feather is a cantilever beam in aluminum, the length have to be higher than a minimum value depending on the elasticity modulus and elasticity limit, not to damage the feathers. The maximum displacement y allows evaluation of the rib contact force F y . In the worth case corresponding to maximum deformation and force, the parasitic friction torque is evaluated.

The pressure induced deformations are static and dynamic. The feather tip deformation is calculated assuming the aerodynamic pressure is applied on only one side of the feather. This assumption is conservative because as the feathers are not airtight, the mean pressure between the feather and the skin is close to the outside mean pressure. The tip calculated displacement is very low (about 50 µm, for L = 15 mm, t = 0.5 mm), even within the conservative assumption. The dynamics aspects are important, feather flutter must not occur. Therefore, the first resonance frequencies of the feathers are estimated to be high enough to avoid fluid-mechanic resonance. This condition is respected for all the feathers.

Comparison of optimized designs

This section compares the two skin solutions. An optimization algorithm has run two times, the first designs a flap with elastic skin and the second, a flap with skin and feathers.

16

Page 16 of 25 AUTHOR SUBMITTED MANUSCRIPT -SMS-106292.R2

A c c e p t e d M a n u s c r i p t

Optimization algorithm

The study aims at designing a feasible and realistic flap by a standard optimization algorithm, not the development of new optimization algorithm. The computations are based on the genetic optimization algorithm from Mathworks' Matlab software [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF]. Optimizations have run on a laptop with 8 core CPU at 2.3 GHz. The cost function evaluation -i.e. computation of a completely sized flap with specified shape, actuators, antagonistic actuators and skin, as presented in Section 4 -takes about 0.6 s. Each optimization lasts about 10 min approximately and is proceeded with a population of N = 150 individuals (initially randomly generated) and a cross-over fraction of 0.9.

As the optimization problem is multi-objective and requires only one criterion, cost weightings are used. The optimization is processed for multiple cost weightings to ensure the suitability of the results.

Optimized designs 5.4. Technological impacts on the objectives

This section compares the results from two technological choices: elastic skins with elastic beam hinges and feathered skins with gliding hinges. We consider a 0.5 m span flap section. This corresponds to reasonable compromise of distributing the actuators every 0.5 m. The estimated mass comes from the actuators (SMA and approximation of structures and anchors), the skin and a simplified structure model. A dedicated thermal study of the actuators provided recommendations on thermal insulation, heating power to actuated and heating power to maintain the actuation temperature. For example, independently of the topology, a 60 s dynamic leads to a maintaining heating power per SMA mass of 293 W/kg. This detailed study is not detailed here for brevity purpose.

Flap with elastic skin with elastic beam hinges

The first considered flap uses the elastic skins modeled in Section 5.1. For a 0.5 m span flap section, the total mass is estimated to 4.9 kg and a maximum heating power of 8 kW. The internal forces are quite large; for example, the first agonistic actuator is designed to apply a force equivalent to 4.3 tons. This high force value on the actuator is due to the aerodynamic forces and internal forces, as shown in table 1. More than half the actuator specified torque is dedicated to counter act the skin elasticity. The part of the force required for the hinge elasticity is below 5 %. This indicates that the actuators can be downsized by a factor 2 if the skin is removed.

Flap with skin and feathers

The second flap uses the skin and feather concept modeled in Section 5.2. Instead of elasticity, this technology parasitic torque is due to friction. For the 0.5 m span flap section, the total mass is estimated to 2.6 kg and a maximum heating power of 2 kW. The internal forces are lower than the solution with elastic skins. For example, the first agonist actuator is designed to apply a force equivalent to 2.0 tons. This force in the actuator is mainly due to the aerodynamic forces, as visible in table 2. The part of the force due to friction from hinges and feathers is below 1.3 %.

Technological choice

For the elastic skin technology, more than half of the actuators' output torques are dedicated to counteract the skin elasticity. For the feathered flap instead, the actuators are designed to only compensate the aerodynamic forces and the antagonistic actuators. In addition, the feathers demonstrate here their ability for a low power consumption of the morphing flap with low internal forces. Therefore this technology is selected.

Final choice, a multi-objective optimized design

The previous optimization results are highly depend on cost weightings. These weightings are empirically tailored and the impact from the shape approximation accuracy or the actuator integration are not distinguished. To discriminate the impact of the different costs, a multi-objective approach based on Pareto front method [START_REF] Deb | Multi-objective optimization using evolutionary algorithms[END_REF] is used in the following.

Multi-objective optimum

Considering a population of N flaps, named X i , i ∈ [[1, N ]], the costs C sha (X i ) -the suffix sha indicating the shape and int the actuators' integration -(C int (X i ) respectively ) are related to the shape and to the integration of actuators respectively. The flap X j is Pareto optimal if there is no flap with a better cost, i.e. ∀i ∈ [

[1, N ]], C sha (X i ) > C sha (X j ) or C int (X i ) > C int (X j ).

18

The cost definitions for this multi-objective optimization are written in Equation 8. It is remembered that: S maxU P /S maxDOW N are the maximum skin error between the reference and the effective deformed shape for up and down deformations; k is the hinge/actuator/antagonistic actuator number of the considered flap, M aero k is the aerodynamic torque on the k th hinge; θ max k is the angle range of the k th hinge, h max k and L max k are the maximum available lever arm and length of the considered agonistic/antagonistic actuator.

C sha (X i ) = S maxU P (X i ) + S maxDOW N (X i ) C int (X i ) = k M aero k θ max k • h max k L max k (8)
To justify the integration cost, consider that the k th actuator is designed with the lowest force. To apply the required torque with the lowest force, the lever arm must be maximum (h max k ). Then the stroke of an actuator is θ max k • h max k . The actuator's length is assumed maximum. The relation between the stroke and the length is L max k = θ max k •h max k / SM A . Depending on the geometric parameters, the SMA strain SM A can be expressed by SM A =

θ max k •h max k L max k
. The actuator is feasible if the SMA strain is low enough to ensure the expected fatigue life. Thus the minimization of cost C int (X i ) tends to minimize the forces in the structure, whilst maximizing the space for actuator integration. A physical interpretation of this cost comes from the fact that is the sum for all the actuators of the product of the actuation work by the lifespan.

Optimization results

Considering these two costs, the genetic algorithm has spread the population of the flap on the Pareto front presented in Figure 10. Here the costs have been normalized, i.e. linear transformed to be between 0 and 1. Two main areas of the Pareto front are identified: an area where the shape approximation costs are low (< 0.3) for a wide range of integration costs, and an area where the integration costs are low (< 0.1) and decreases with the increase of the shape costs. At the extremity of these two areas are located the best shape approximation flap and the best actuator integration capability flap. These two designs are drawn on the figure with red highlighted drawbacks. For the flap with the best shape cost, the integration is not carried out for one actuator. For the flap with the best integration cost respectively, the shape approximation is very poor in the highlighted area. This observation indicates that an optimized flap which has both the best shape approximation and the best actuator integration capability does not exist. A compromise must be done.

Corresponding to the Pareto front flaps, the locus of the hinge positions is represented on Figure 11. The locus of hinges #5 does not move, locus for hinges #1, #3 and #4 are quite continuous whereas the positions of hinges #2 are spread on two spots. The Figures 11b and11c indicates that the best hinge positions corresponding to the best shape approximation or the best actuator integrability are not identical. A compromise must be done. The design of the two flaps discussed on Figure 10 have different articulations' positions. Finally, a fine design of the actuators for both considered flaps shows that both are feasible with comparable performance in terms of weigh, lifespan and power consumption. Therefore the choice of the ultimate flap design is the flap with the best shape approximation. The articulation positions corresponding to the final choice are indicated on Figure 11.

Design choice

The justification of the selected design is presented in the previous Section 6. The hinge positions are presented in Figure 11. The design estimates a 9 kg flap section with an average consumption of 0.515 kW. One can notice that the expected fatigue life is limited to 40,000 cycles, which is very low. This estimation comes from the lever arm limitation in the region near the trailing edge. This issue is alleviated during detailed design by inclining the actuators: a new calculation of the active lengths and lever arms is carried out for such inclined actuators.

Additionally, all the actuators have been inclined to increase the lever arm and then to decrease the forces in the structure. Taking into account the available space in the span direction, the different actuators are distributed in the span. The result is a macro-actuator composed of the different actuators, the articulated ribs and the hinges. The flap is then made by the assembly of four macro-actuators plus spars, stiffeners, skin and feathers. Figure 12 presents a view of the manufacturing CAD model of the macro-actuator, without the skins.

Finally, one can remark the estimated 9 kg are not realistic. It corresponds to the weight of the complete actuators, the hinges and rough estimations of skins and structures. The precise weights of the structure, the skins and devices like electronics are not known, even if they could be non-negligible. Nevertheless, these neglected elements elements do not affect the worthiness of the design choice. Indeed, the comparison of the relative the added internal forces, added power consumption and added weight are enough to compare the technologies and make a choice.

Conclusion

The present article treats a morphing wing design enabled by optimization. The target is a true scale flap of a civil aircraft (e.g. Airbus A320 type). It is noticeable that the designed flap can carry realistic aerodynamic loading of several tons. A camber control concept based on articulated ribs is proposed. This concept can be used with different hinge technologies and different skin technologies. The actuators, based on innovative smart materials like shape memory alloys, are also investigated. Comparing different actuation topologies, it has been shown that in terms of weight and internal forces, the agonist-antagonist topology is better than a solution with one actuator and a counter spring. An algorithm that automatically designs a working morphing flap from input parameters is proposed. This algorithm has been used within an optimization process to compare the different technologies. As a result, the elastic skin causes

22

A c c e p t e d M a n u s c r i p t large issues according to the internal forces and is responsible for over-sizing by more than 100% the actuator forces. Therefore, an innovative bio-inspired skin and feather concept has been selected, taking advantage of the specifications and constraints. As the optimization objectives are plentiful (weight, shape approximation, energy consumption, fatigue life), a multi-objective optimization based on Pareto front is used to choose the ultimate flap design. The resulting design is proven feasible and a detailed realization of the flap has been achieved. This is possible because of the collaborative work with industrials to select aircraft-compatible technologies. Future works will focus on the design, integration, and control of the macro-actuators in a full scale flap prototype. Then electromechanical and aerodynamic testing will be accomplished on dedicated static test bench and wind-tunnel platform.

3 A

 3

Figure 1 :

 1 Figure 1: The two objective shapes superimposed with the original non-deformed profile.

4 Page 4

 44 , is decomposed fourfold. All the following elements are worth designing wisely : of 25 AUTHOR SUBMITTED MANUSCRIPT -SMS-106292.R2A c c e p t e d M a n u s c r i p t• articulated ribs define the geometry and carry the other components. They have to withstand the internal and external (i.e. aerodynamic) forces, whilst being low weight.

Figure 2 : 5 Page 5

 255 Figure 2: 2D illustrative sketch of the proposed concept. Articulated ribs are placed between fix leading edge and trailing edge. Actuators and specific skins are located within the morphing flap. The aerodynamic force distribution is indicated by the blue arrows.

6 Page 6 Figure 3 :

 663 Figure 3: Tensile test of SMA wires at two temperature. n represents the life cycle of the SMA actuators. Areas n < 20k, n > 20k and n > 50k respectively indicate that for actuators working in these areas, less than 20, 000, more than 20, 000 and more than 50, 000 cycles can be expected. The specified working area is lower than 200 MPa and 2%.

Figure 4 :

 4 Figure 4: Sketch of the two compared solutions.

9 A

 9 kg) Actuated torque (N.m) No -return system Counter spring Agonist-antagonist

Figure 5 :

 5 Figure 5: Performance of the counter spring and the antagonistic solutions compared to the no return system.

  Thus the agonist-antagonist 10 Page 10 of 25 AUTHOR SUBMITTED MANUSCRIPT -SMS-106292.R2 A c c e p t e d M a n u s c r i p t

Figure 6 :

 6 Figure 6: Nested control loops for antagonistic SMA actuators.

11 AFigure 7 :

 117 Figure 7: Schematic representation of the flap sizing algorithm on the left. The top flap pattern presents the input data: the hinge positions. The middle drawings present the shape approximation evaluation. The bottom flap pattern shows the actuator implementation, assuming an agonist-antagonist topology.

13

 13

Figure 8 :

 8 Figure 8: Morphing flap with elastic skin. Skin is hung when represented in green, whereas it is free and taut when represented in orange. The right part presents a detailed view of a parametrized taut skin element.

15

 15

Figure 9 :

 9 Figure 9: Skin and feather concept applied to one articulation. The flaps presented correspond to neutral and maximum cambered shape. A schematic view of the mechanism is also presented on right.

19

 19

Figure 10 :

 10 Figure 10: Pareto front normalized Shape cost versus normalized Integration cost. The population is represented by blue points and the Pareto optimal individuals by red circles. The two selected best compromise flaps are represented.

Figure 11 :

 11 Figure 11: Hinge locus. The positions of the hinges of the Pareto front's flaps are plotted threefold: 11a color shades represent the hinge number, 11b color shades represent the actuator integration cost and 11c color shades represent the shape approximation cost.

21 Page 21

 2121 of 25 AUTHOR SUBMITTED MANUSCRIPT -SMS-106292.R2 A c c e p t e d M a n u s c r i p t (a) 3D CAD view of the detailed design of the macro-actuator of the flap. Skins and feathers are not represented (b) Side view of the detailed design of the macro-actuator of the flap.

Figure 12 :

 12 Figure 12: Detailed design views of the macro-actuator of the flap.

Table 1 :

 1 Design moments for the actuators at hinge locations (not including the torques from agonistic and antagonistic actuators)

	Page 17 of 25	AUTHOR SUBMITTED MANUSCRIPT -SMS-106292.R2

Table 2 :

 2 Design moments for the actuators at hinge locations (without the torques from agonistic and antagonistic actuators)

http://www.saristu.eu

http://www.cleansky.eu/

http://www.flxsys.com/ 2

Acknowledgements

The authors are grateful to Airbus that provided partly funding for this research. The authors are thankful to Dominique HARRIBEY for his valuable contribution, as well as Oussama FILALI and Maxime RUFFEL for their useful advice and help within this work. The authors thank the technical services of IMFT and LAPLACE for their assistance in the realization of this study.

This work is part of the Smart Morphing and Sensing for aeronautical configurations. This project has received funding from the European Unions H2020 program for research, technological development and demonstration under grant agreement no 723402.


Othmane Zerhouni

Maria-Gabriella Tarantino

Kostas Danas

Influence of the internal geometry on the elastic properties of materials using 3D-printing of computer-generated random microstructures

Understanding elastic properties of rocks is a scientific challenge due to the complexity of their microstructures. This study combines the numerical tools to generate models for internal geometry of pores with the 3D printing technology in order to control the shape and size of pores as well as the distribution of pore's network inside the sample. Accuracy of the printing has been assessed by optical microscopy. The numerical and experimental tests conducted on generated microstructures show that the elastic properties are independent of the size of the pore. In turn, the shape of the pore has a strong effect on the elastic properties. This study shows that in case of multiple pore types, the ones with small aspect ratios have a strong impact though of minor volume fraction. The methodology that is developped can be extended to investigate the influence of more parameters in case of connected porosity and tests theories that have been proposed to link physical properties of reservoir rocks to their internal geometry.

INTRODUCTION

An important challenge in the developpement of reservoir models with high predictive capacity of production and the discovery of new wells is to understand and quantify the relationship between the internal structure of the rock and their elastic properties. The porous space of sedimentary rocks is characterized by heterogeneous and complex structures that spread over multiple scales. This multiscale geometry is difficult to capture with imaging techniques that are limited by the sample scales that can be observed and therefore need upscalling approach to get the overall properties of the field, i.e. need the determination of the equivalent permeability and bulk properties of each scale.

Laboratory analysis of rock properties, such as core flooding or mercury porosimetry, usually involves destructive testing which prevent from re-using the same rock in multiple tests. In turn, the numerical simulation offers a faster and nondestructive way to compute the elastic properties. However, the accurate characterization depends on the careful selection of the representative volume to overcome the dependence of the properties on the boundary conditions and the comparison between the experimental and computational testing does not involve the same scale of porosity due to the different size of samples that are used for each method.

The development of 3D printing manufacturing process has led to the development of new tools in geoscience to advance knowledge of the link between properties and microstructure. Such process has been used to derive the influence of the pro-cesses in sedimentary rock on the their experimental properties [START_REF] Ishutov | Using resin-based 3d printing to build geometrically accurate proxies of porous sedimentary rocks: Groundwater[END_REF][START_REF] Head | Effects of changes in rock microstructures on permeability: 3-d printing investigation[END_REF]. Their approach is based on image manipulation of the porous network of real rocks in order to highlight the microstructure resulting from a specific process such as compaction or dissolution. Their resulted images are then printed and tested in order to compare between experiments on printed samples and numerical simulations on images. The preliminary results shows good quantitative agreement for permeability related to the calculated porosity which shows the potential of additive manufacturing in developping methods to characterize the porous network influence on rock's properties and compare the experimental tests and the numerical simulation with the same samples. Another potential use of 3D printing consists of manufacturing the mathematical models of fractured rocks and porous network in order to validate the theoritical studies and the numerical methods [START_REF] Suzuki | fracture network, flow experiment, 3d printer, ct scan images: Stanford Geothermal Workshop[END_REF]. By controlling the fracture and pore geometry in the mathematical models which is in turn transferred into the experimental samples, this approach provide quantitative results of the influence of several aspects of the geometry such as the aperture of the interconnected porosity and the size of fractures on the flow properties.

The present study investigates the influence of the pore geometry on elastic properties of rocks. To that aim, we create random porous microstructures with controlled porosity and well-defined size and shape of pores. These microstructures are then experimentally and numerically tested to measure the linear elastic properties of the porous specimens. Using this methodology, we separately investigate the influence of each of the porosity, the size and the shape of voids on the elastic response of porous microstructures.

MATERIAL AND METHOD

We use a random sequential adsorption (RSA) algorithm [START_REF] Rintoul | Reconstruction of the structure of dispersions[END_REF][START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Lopez-Pamies | The nonlinear elastic response of suspensions of rigid inclusions in rubber: Ii-a simple explicit approximation for finiteconcentration suspensions[END_REF] to construct the virtual microstructures of the present study. We specify to the algorithm the total porosity c and parameters of each pore type i which are its size a i , its aspect ratio α i and its volume fraction c i to create ellipsoidal voids randomly distributed in a cubic cell of size L 3 . For illustration purposes, Fig. 1 shows three such RVEs comprising monodisperse spherical voids of increasing porosity c = {10, 20, 30}%.

We build an experimental protocol to determine the elastic properties of the 3D-printed controlled microstructures. To that aim, we design virtual test specimens which ensure uniaxial stress conditions in the gage section by assembling lengthwise five representative cubic cells which are enclosed between the heads allowing us to mount the specimen onto the uni- axial machine. These specimens are then 3D -printed using an acrylic photopolymer available in our EDEN 260VS 3Dprinter purchased by Stratasys. We determine the effective Young's modulus and Poisson's ratio of the 3D-printed random porous microstructures using uniaxial tensile loading while measing the axial and transverse strains by the extensometers as shown in Fig. 2. The elastic properties of the printing material are isotropic with Young modulus E = 1.4 GPa and Poisson's ratio ν = 0.42 which correspond using equivalence between the P-wave velocity V p and S-wave velocity V s with the static bulk modulus κ and shear modulus µ in the case of isotropic material :

V p = κ + 4/3µ Φ , V s = µ Φ (1)
The printing material wave velocities are around V p = 1750m.s -1 and V s = 660m.s -1 and density Φ = 1170kg.m -3 .

In turn, we use finite-element methods with a conforming quadratic mesh to obtain numerically the overall elastic properties of generated microstructures. Subsequently, the microstructures are discretized using the automatic mesh generator NETGEN [START_REF] Schöberl | Netgen an advancing front 2d/3d-mesh generator based on abstract rules[END_REF] and simulated using ABAQUS software with periodic boundary conditions. For each set of parameters, we create four different samples to address the statistical deviation of the RSA process.

The host material is modeled by an isotropic linear elastic constitutive behavior with Young's Modulus and Poisson's ratio corresponding to those of the virgin matrix phase used during the printing.

RESULTS AND DISCUSSION

Effect of porosity on the elastic properties

We first focus our experimental and 3D printing efforts to isotropic random microstructures consisting of single sized (i.e. monodisperse) spherical voids embedded in a homogeneous solid matrix in order to study the effect on elastic properties only due of the variation of the total porosity. We compare of their elastic effective properties with the HS isotropic bounds [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF][START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] given by

κ HS = 4(1 -c)µ m κ m 4µ m + 3c κ m , µ HS = (1 -c)(8µ m + 9κ m )µ m 4µ m (2 + 3c) + 3κ m (3 + 2c) (2)
where κ m and µ m denote the bulk and shear moduli of the matrix phase, respectively and c the volume fraction of the inclusion phase or porosity in the present context. Results in Figs 3a and3b show that increasing the porosity creates larger void regions which implies significant interactions between voids and softens the material.

Pore size effect

Pore size is pointed out as an important parameter in the variation of P-wave and S-wave velocities during pore-type characterization that links between microscopic images from thin section and experimental tests over the core sample. However, none of the existing analytical estimates take into account the size of the pore in their predictions but only the shape of pore and the porosity. In order to determine if pore size affects the elastic properties, we generate microstructures that have the same porosity c = {5%, 10%, 15%, 20%, 25%} and spherical pores of decreasing size. In order to link between the pore-type distribution and the elastic properties of a porous material, one approach by [START_REF] Weger | Quantification of pore structure and its effect on sonic velocity and permeability in carbonates[END_REF][START_REF] Fournier | Elastic properties of microporous cemented grainstones[END_REF] consists in deriving two parameters (α κ and α µ ) corresponding to the aspect ratio of ellispodal pores in the different effective medium (DEM) model [START_REF] Berryman | Single-scattering approximations for coefficients in biot's equations of poroelasticity[END_REF][START_REF] Mavko | The rock physics handbook: Tools for seismic analysis of porous media[END_REF]) that has similar P-wave and S-wave velocities to the well-log acoustic response of the microstructure. The accuracy of pore-type classification based on such approach is closely dependent on understanding the most influent pore type among the pore types withdrawn, from thin section observations. In order to quantify the uncertainties related to this choice, we generated microstructures with controlled porosity c ∈ [1%, 25%].

We use aspect ratios taken from petroacoustic parameter α κ obtained by studying carbonates rocks [START_REF] Fournier | The equivalent pore aspect ratio as a tool for pore type prediction in carbonate reservoirs[END_REF] to generate microstructures with different pore shapes and volume fraction of each. Realizations of such microstructures are shown in Fig. 6 to highlight the different internal geometry that can be found for the same porosity c = 20%. We use numeri- 7b highlights the effect of pore shape on the elastic properties of the material. However, the results shown are assuming isotropic effective behavior to retrieve bulk and shear modulus and subsequently P-wave and Swave velocities. This assumption is far from being consisting with the elastic properties of the microstructures that exhibits high anisotropy in case of small AR.

To highlight the importance of this study, we compute the deviation from the effective result if considering that the most present pore shape is the most influencing on the elastic properties. Results in Fig. 8 show that small aspect ratio are more impacting the elastic properties and that the acoustic signature of the microstructure is highly related to the presence of these highly elongated pores.

ACKNOWLEDGMENTS

Authors acknowledge support from TOTAL S.A. The computational part of this work was also supported by the ANR under contract number ANR-10-EQPX-37. Figure 8: Dominant phase Error in the elastic properties when considering a phase dominant over another in two pore type coexisting in the same microstructure.

Figure 1 :

 1 Figure 1: Representative cells of 12.8 mm (correponding to 0.5in) with N randomly distributed spherical particles of monodisperse sizes in a total porosity (a) c = 10% and N = 160 (corresponding to a = 0.68 mm), (b) c = 20% and N = 275 (corresponding to a = 0.73mm), (c) c = 30% and N = 400 (corresponding to a = 0.72). (d-f) Representative meshes corresponding to the undeformed configuration of the representative cubic cells.

Figure 2 :

 2 Figure 2: Experimental setup consisting of sample gripped by self-tightening grips mounted on a uni axial tensile machine which is connected to a LabView software to acquire the macroscopic measures of the axial and transverse strain.

  Figure 4 shows different realizations of these microstructures for c = 20% with pore size being the only changing parameter. FE results in Fig 5 are very similar for all computed sizes showing no effect of size on the elastic properties of the material.The variations that can be observed are due to statistical sampling.

Figure 3 :

 3 Figure 3: Effective (a) P-wave velocity V p , (b) S-wave velocity V s of microstructures with spherical single-size voids as a function of the porosity c. Comparison between the analytical HS bounds, FE monodisperse numerical estimates and experimental results.

Figure 6 :

 6 Figure 6: Microstructures with porosity c = 20% composed of pores of different aspect ratios (AR) : a) 100% of ellipsoidal pores of AR= 0.3, b) 100% of ellipsoidal pores of AR= 0.13, c) 100% of ellipsoidal pores of AR= 0.1, d) 90% of ellipsoidal pores of AR= 0.3 and 10% of ellipsoidal pores of AR= 0.1

Figure 7 :

 7 Figure 7: Results of :a) P-wave velocity b) S-wave velocity, variations due to change in the aspect ratio of the ellipsoidal voids.


O Zerhouni
email: othmane.zerhouni@polytechnique.edu

M G Tarantino

K Danas

C N R S Lms

Experimental homogenized elastic properties of computer-generated 3D-printed random porous material

Keywords: Additive manufacturing, Hashin-Shtrikman bounds, Random Sequential Adsorption, Random composites, Homogenization, Porous Materials, Effective properties

The present study introduces a methodology that allows for the combination of 3D-printing, experimental testing, numerical and analytical analysis of random porous materials with controlled homogenized elastic properties. These microstructures are computer-generated based on a random sequential addition algorithm with statistically controlled morphological properties such as volume fraction, shape and size of voids as well as isotropic distribution of their centers. We first focus on porous material containing single-size (monodisperse) spherical voids. The porous specimens are fabricated by 3D printing with polymer jet technology and then microstructures are aposteriori investigated by optical microscopy and SEM. The influence of the 3D printing process parameters is also experimentally assessed. An experimental protocol is developed in order to obtain the effective elastic properties of the porous specimens.

Introduction

Porous materials are present in nature (e.g., rocks and geomaterials) but can also be designed to allow for controlled stiffness of lightweight structures. Recent development of additive manufacturing processes has led to new methods for generating complex microstructures that can achieve high-performance lightweight properties. The porous materials fabricated by 3D printing find applications in high-stiffness lightweight structures [START_REF] Zok | A protocol for characterizing the structural performance of metallic sandwich panels: application to pyramidal truss cores[END_REF][START_REF] Berger | Mechanical metamaterials at the theoretical limit of isotropic elastic stiffness[END_REF], acoustic and vibration dampers [START_REF] Göransson | Acoustic and vibrational damping in porous solids[END_REF][START_REF] Banhart | Damping properties of aluminium foams[END_REF][START_REF] Ma | A nonlinear auxetic structural vibration damper with metal rubber particles[END_REF], impact energy absorbers [START_REF] Davies | Metallic foams: their production, properties and applications[END_REF], high electric capacitors [START_REF] Wang | 3d aperiodic hierarchical porous graphitic carbon material for high rate electrochemical capacitive energy storage[END_REF] and filtration [START_REF] Alderson | An auxetic filter: a tuneable filter displaying enhanced size selectivity or defouling properties[END_REF] among others.

In this study, we are interested in transferring theoretical and virtual microstructures that are very close to the Hashin-Shtrikman bounds numerically to fabrication and experimental testing.IThis aims at introducing a methodology that allows to combine 3D printing, experimental testing, numerical and analytical modeling in order to create random close-cell porous materials with effective elastic properties that approach the Hashin-Shtrikman bounds. In this study, we focus on isotropic random microstructures consisting of single sized (i.e. monodisperse) spherical voids embedded in a homogeneous solid matrix. Those microstructures are shown to have almost identical effective elastic properties with multiple size (i.e. polydisperse) microstructures for volume fractions up to 30%.

Generation of random isotropic porous materials

RSA generation of isotropic random spherical inclusions

We use the random sequential adsorption algorithm [START_REF] Rintoul | Reconstruction of the structure of dispersions[END_REF][START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF][START_REF] Lopez-Pamies | The nonlinear elastic response of suspensions of rigid inclusions in rubber: Ii-a simple explicit approximation for finite-concentration suspensions[END_REF] to construct the virtual microstructures of the present study. We use as input parameters the porosity c and the number of pores N to create spherical voids that are distributed randomly in a cubic cell of volume L 3 . In thiss study, we set the minimal distance between two voids at 150µm. We note in passing that this distance can be further reduced to allow larger volume fractions.

However, a minimum spacing between non connected voids is needed so that they can be properly manufactured by our 3D printer. For illustration purposes, Fig. 1 shows three such RVEs comprising monodisperse spherical voids of increasing porosity c = {10, 20, 30}%.

RVE size and related boundary conditions

The choice of N is made in order to have the ratio D/L, with D and L denoting the diameter of the voids and the edge length of the cubic cell, large enough to obtain numerical effective properties that are independent on boundary conditions. [START_REF] Suquet | Homogenization Techniques for Composite Media: Lectures Delivered at the CISM International Center for Mechanical Sciences[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] These are computed from the average stress and strain fields using the overall constitutive equation [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF]

σ(x) = C : ε(x) (1)
where ε(x) = ε is the macroscopic imposed strain, C denotes the overall stiffness tensor and the notation

. is used to denotes volume averages. In practice, adopting Voigt notation, each column of the stiffness tensor C is computed by imposing a strain field in a specific direction [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]. For example, imposing the overall strain field ε ≡ (ε 11 , 0, 0, 0, 0, 0), we obtain six linear relations: 11,22,33,12,23,31. The computation of all six columns of C is done by application of six independent loads ε as described in [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF].

C kl11 = σ kl /ε 11 with kl =
The matrix phase is modeled by an isotropic linear elastic constitutive behavior with Young's Modulus E = 1.4 GPa and Poisson's ratio ν = 0.42 corresponding to those of the matrix phase used in the 3Drinting. Results of the experiments on microstructures from 1% to 30% porosity are shown in Fig. 2. 3

Additive manufacturing and experimental methods

In order to construct the virtual geometry of the porous samples, we adopt the protocol shown in Fig. 3.

We first assemble length-wise five representative cubic RVE cells of length L = 12mm in order to build the reduced uniform section of the test specimens. The latter is enclosed between the heads of the sample which, in turn, have a solid section thus allowing us to mount the specimen onto the uniaxial machine.

We use commercial software NETGEN to transform the 3D virtual model of our test specimens into a stereolitography format (i.e. STL) to 3D-printing. We assess the dimensional accuracy of the 3D-printed porous microstructures by use of an optical microscope. Overall, Fig. 4 reveals that pores with circular section can be manufactured with very good dimensional accuracy within the range of size required by our numerical RVE analysis for L = 12mm. Finally, it is worth mentioning that no additional micron-size porosity has been observed in the matrix phase as a result of the printing process. The main aim of the experiments is then the measurement of the effective Young's modulus and Poisson's ratio of the 3D-printed random porous microstructures. Following the experimental methodology presented in [START_REF] Hossain | Experimental study and numerical modelling of vhb 4910 polymer[END_REF], we measure the elastic parameters by means of a multi-step relaxation testing, as shown in Fig. 5a, to overcome any effect of the viscosity in the material. To assess reproducibility of the data, we use four specimens for each test.

The effective Young's modulus, E (or E m for the pure polymer), is evaluated as the slope of the line connecting the equilibrium (fully relaxed) stress points in the nominal stress-strain curve as shown in Fig. 5b and the calculation of the Poisson's ratio is found by the slope of transverse strain-axial strain curve.

a) b) c) The experimental results for the matrix and support material following this procedure provide experimental evidence that the matrix material (VeroWhite) is statistically homogeneous and exhibits an elastically isotropic behavior with a Young's modulus E m = 1400 ± 120MPa and a Poisson's ratio of ν m = 0.42 ± 0.02.

These moduli correspond to a bulk modulus of κ m = 2920MPa and a shear modulus of µ m = 493MPa. In turn, it indicates that the support material, used by 3D printer to fabricate the spherical pores, has negligible elastic properties when compared with those of the matrix material (VerowhitePlus) thereby having no effect on the effective porous material response.

Results and Discussion

For isotropic porous materials, the Hashin-Shtrikman bounds [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF][START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] are given by

κ HS = 4(1 -c)µ m κ m 4µ m + 3c κ m , µ HS = (1 -c)(8µ m + 9κ m )µ m 4µ m (2 + 3c) + 3κ m (3 + 2c) (2)
where κ m and µ m denote the bulk and shear moduli of the matrix phase, respectively and c the volume

Concluding remarks

In this work, we have investigated the influence of some geometrical features of the microstructure on the effective elastic propertiesof random isotropic close-cell porous materials by means of a numericallyaided additive manufacturing procedure. Results show that the size of the pores does not affect the overall experimental elastic properties of the 3D printed material and that strain localization seems to cause the effective properties of the material to depart from the Hashin-Shtrikman bound.

The present study can be extended to larger porosities, at least up to c = 70% (or relative density ρ = 1 -c = 0.3) using polydisperse distributions of particles.

Figure 1 :

 1 Figure 1: RVEs of unit volume L 3 with N randomly distributed spherical particles of monodisperse sizes for a total porosity (a) c = 10% and N = 160, (b) c = 20% and N = 275, (c) c = 30% and N = 400. (d-f) Representative meshes corresponding to the undeformed configuration of the representative cubic cells.

Figure 2 :

 2 Figure 2: (Right axis) Converged pore to cube size ratio D/L and (left axis) pore diameter (in mm) for a cube side length L = 12mm as a function of the porosity c.

Figure 3 :Figure 4 :

 34 Figure 3: (Top) Virtual testing sample generated by assembling length-wise 5 cubic porous RVEs and by adding gripping heads of solid material. (Left) A zoom of the RVE which defines our gage section. (Bottom) 3D printed testing sample after cleaning support material from open pores and boundary surface. (Right) A zoom of the 3D printed RVE showing the size of the actual axial and transverse gage sections in the experiment.

Figure 5 :

 5 Figure 5: Overall applied strain rate ε = 10 -5 s -1 : (a) Stress-Strain curves obtained during multi-step relaxation steps. The Young's modulus is obtained by the slope of the curve that connects the fully relaxed stress states excluding the first point. (b) Axial versus transverse strain curve during the mult-step relaxation test. Slopes during loading are found to be fairly independent of the level of the axial strain.

Figure 6 :

 6 Figure 6: Normalized with the matrix properties effective (a) Young's modulus E/Em, (b) Poisson's ratio ν as a function of the porosity c. Comparison between the analytical HS bounds, FE monodisperse numerical estimates and experimental results.