J Serra
email: joelserra@isae-supaero.fr

A Trellu

C Bouvet

S Rivallant

B Castanié

L Ratsifandrihana

Combined loadings after medium velocity impact on large CFRP laminated plates: Discrete ply model simulations

Keywords: Composite structures, Structural testing, Medium velocity impact, Combined loading, Numerical modeling

To decrease development costs of aeronautical composite structures, industrials and academic researchers are turning to "virtual testing" methods. To help achieve this objective, a new methodology has been developed at the Institut Clément Ader. It is based on a complex loading test rig for technological specimens (more representative than the usual level of coupons) and the Discrete Ply Model (DPM), a model allowing the main damage occurring in composite structures to be simulated. This research focuses on the combined loading of CFRP plates after impact. In a previous work, plates were experimentally impacted and loaded after the impact in the VERTEX test rig and a methodology to transfer boundary conditions to the numerical model was validated using a simple model without damage. In the work presented here, the DPM was used to simulate both the impact and the complex loading after impact tests. Numerical results and experimental results showed satisfactory correlation. The research undertaken demonstrates that the impact damage has very little influence when the tests are performed at larger scales than coupons. It is therefore paramount to use a pyramidal approach to design and test composite structures, in order to avoid oversizing.

Introduction

Composites present high mechanical properties and are now widely used in aeronautical structures but they are sensitive to low velocity impacts. These impacts, due to unexpected dropping of tools during aircraft ground operations or the presence of small debris on the landing area [START_REF] Morteau | Composites @ Airbus, Damage Tolerance Methodology. FAA Workshop for Composite Damage Tolerance and Maintenance[END_REF], for example, can significantly reduce the residual strength of the structure even if the damage is not visible (Figure 1). The residual compression strength after impact is the most critically affected due to the buckling of delaminated plies [START_REF] Abrate | Impact on Composite Structures[END_REF][START_REF] Perugini | Influence of delamination growth and contact phenomena on the compressive behavior of composite panels[END_REF] and can fall by 50% [START_REF] Tropis | Certification of the composite outer wing of the ATR72[END_REF][START_REF] Caprino | The significance of indentation in the inspection of carbon fibrereinforced plastic panels damaged by low-velocity impact[END_REF]. This impact damage issue imposes a damage tolerance approach in the design of composite structures in order to ensure that the structure will withstand in-service loads even if the damage is undetectable. This design philosophy has led to standards for composite coupons under low velocity impact [START_REF] Astm D ; M | Standard Test Method for Measuring the Damage Resistance of a Fiber-Reinforced Polymer Matrix Composite to a Drop-Weight Impact Event[END_REF] and compression after impact [START_REF] Astm D ; M | Standard Test Method for Compressive Residual Strength Properties of Damaged Polymer Matrix Composite Plates[END_REF].

Simulation of failure in composites

The complex failure phenomena affecting composites depend on several parameters, such as the stacking sequence, the nature of the fibers, the ply thickness, etc., and thus prohibit the use of analytical or semi-analytical models [START_REF] Abrate | Impact on Composite Structures[END_REF][START_REF] Zhou | The use of experimentally-determined impact force as damage measure in impact damage resistance and tolerance of composite structures[END_REF][START_REF] Suemasu | Compressive behavior of multiply delaminated composite laminates part 1: experiment and analytical development[END_REF][START_REF] Suemasu | Compressive behavior of multiply delaminated composite laminates part 2: finite element analysis[END_REF][START_REF] Olsson | Analytical prediction of large mass impact damage in composite laminates[END_REF][START_REF] Huang | Analytical modeling of impact resistance and damage tolerance of laminated composite plates[END_REF] alone in the design of structures subject to low velocity impact. These methods are coupled with experimental campaigns following the pyramid of tests design philosophy [START_REF] Rouchon | Certification of Large Aircraft Composite Structures, Recent Progress and New Trends in Compliance Philosophy[END_REF]. To replace these long and costly experimental campaigns (tens of thousands of tests per aircraft), numerical simulations based on the Finite Element (FE) method are being developed and can simulate damage occurring during low velocity impact and compression after impact. This is the "Virtual Testing" approach [START_REF] Yang | Virtual Testing for Advanced Aerospace Composites: Advances and Future Needs[END_REF].

The main damage occurring within the ply (intra-laminar) during a low velocity impact is matrix cracking and fiber breakage. In addition, delamination can develop between two plies (inter-laminar) [START_REF] Lopes | Low-Velocity impact damage on dispersed stacking sequence laminates. Part II: Numerical simulations[END_REF][START_REF] Lopes | Physically-sound simulation of low-velocity impact on fiber reinforced laminates[END_REF][START_REF] Wisnom | Modelling discrete failures in composites with interface elements[END_REF][START_REF] Menna | Numerical simulation of impact tests on GFRP composite laminates[END_REF]. Fiber breakage and delamination are the most critical types of damage for laminate failure. These two failures decrease the laminate stiffness [START_REF] Abrate | Impact on Composite Structures[END_REF][START_REF] Zafer Kutlu | Modeling Compression failure of laminated Composites Containing multiple through-the-width delaminations[END_REF] as broken fibers cannot absorb forces and delamination principally reduces compression and bending resistance. Even if very low velocity impacts do not induce fiber failures, the delamination can lead to strongly decreased structural efficiency [START_REF] Aslan | Buckling behaviour and compressive failure of composite laminates containing multiple large delaminations[END_REF][START_REF] Dang | A numerical study on impact and compression after impact behavior of variable angle tow laminates[END_REF]. Matrix cracking is not critical for the safety of a structure but is the first damage occurring during an impact [START_REF] Abrate | Impact on Composite Structures[END_REF][START_REF] Choi | A model for predicting damage in graphite/Epoxy laminated composites resulting from low-velocity point impact[END_REF] and plays an important role in the impact damage scenario. Bending of the global structure induces tension in bottom plies, which leads to matrix cracking. These cracks propagate into the ply thickness and create delamination when they meet a ply with a different orientation. It is therefore important to represent these three main types of damage (fiber failure, matrix cracking and delamination) and their coupling efficiently.

Various methods to model such interlaminar and intralaminar damage have been described.

 Interlaminar (delamination): The literature principally mentions three different methods to simulate delamination. The Virtual Crack Closure Technique (VCCT) was introduced by Krueger [START_REF] Krueger | The virtual crack closure technique: history, approach and applications[END_REF] and, in particular, used by Li et al. [START_REF] Li | Modelling damage of multiple delaminations and transverse matrix cracking in laminated composites due to low velocity lateral impact[END_REF] to propagate an existing crack node by node, comparing the energy required to close the crack to the material tenacity. However, the VCCT suffers from some difficulties: it can only simulate the crack propagation and requires topological information from nodes ahead and behind the crack front to calculate the energy required to close the crack. This need for nodal information involves a remeshing for crack propagation. The Cohesive Zone Model (CZM), introduced by Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt [27], has been used to create cohesive interface elements located in a crack zone. The use of CZM involves knowing the crack location, which is not a problem with laminates because delamination appears between plies with different orientations. Cohesive interface elements can predict both crack initiation and propagation. The CZM is widely used in the literature to simulate delamination occurring during low velocity impact [START_REF] Lopes | Physically-sound simulation of low-velocity impact on fiber reinforced laminates[END_REF][START_REF] Faggiani | Predicting low-velocity impact damage on a stiffened composite panel[END_REF][START_REF] Shi | Modelling damage evolution in composite laminate subjected to low velocity impact[END_REF][START_REF] Aymerich | Prediction of impact-induced delamination in cross-ply composite laminates using cohesive interface elements[END_REF][START_REF] Sun | Interaction of inter-and intralaminar damage in scaled quasi-static indentation tests: Part 2 -Numerical simulation[END_REF]. Studies of cohesive parameters, such as interface stiffness and mesh size effects, to ensure good delamination behavior are proposed in the literature [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF][START_REF] Harper | A study on the influence of cohesive zone interface element strength parameters on mixed mode behavior[END_REF] and lead to the conclusion that a short cohesive zone enables better results. Menna et al. have employed 3D elements to represent laminate plies and model delamination by using a contact "surface-to-surface tiebreak" between any two plies [START_REF] Menna | Numerical simulation of impact tests on GFRP composite laminates[END_REF]. The contact is driven by interface strength parameters. This third method is used but is not as common as cohesive interfaces. It has the advantage of not needing coincident meshes for plies.

 Intralaminar (matrix cracking and fiber breakage): intralaminar damage modeling is most often based on continuum damage theories. Damage is not represented geometrically and kinematically; only its effect on the mechanical properties is taken into account. When an initiation criterion is reached, the laminate stiffness is degraded in the direction considered until failure occurs. Several researchers [START_REF] Lopes | Low-Velocity impact damage on dispersed stacking sequence laminates. Part II: Numerical simulations[END_REF][START_REF] Lopes | Physically-sound simulation of low-velocity impact on fiber reinforced laminates[END_REF][START_REF] Ji | A progressive damage model for predicting permanent indentation and impact damage in composite laminates[END_REF][START_REF] Donadon | A progressive failure Model for composite laminates subjected to low velocity impact damage[END_REF][START_REF] Maimi | A continuum damage model for composite laminates: Part I -Constitutive model[END_REF] have used the continuum damage model approach to simulate both matrix cracking and fiber breakage in 3D elements with a non-oriented mesh. However the limit of this method concerning damage locations should be underlined: the model response is strongly dependent on the mesh size. To overcome this issue, Pinho et al. [START_REF] Pinho | Fracture toughness of tensile and compressive fibre failure modes in laminated composites[END_REF] and Bouvet et al. [START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF] simulated the fiber failure by dissipating a surface energy that is independent of the mesh size, as introduced by Bazant & Oh [START_REF] Bazant | Crack band theory for fracture of concrete[END_REF]. Another way to represent the matrix cracking is to use cohesive elements in the ply thickness between 3D elements oriented in the fiber direction. This method has been used by several researchers [START_REF] Lopes | Physically-sound simulation of low-velocity impact on fiber reinforced laminates[END_REF][START_REF] De Moura | Modelling the interaction between matrix cracking and delamination in carbon-epoxy laminates under low velocity impact[END_REF][START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[END_REF] who assume that the matrix cracking location is known. Moura et al. [START_REF] De Moura | Modelling the interaction between matrix cracking and delamination in carbon-epoxy laminates under low velocity impact[END_REF] place cohesive elements where matrix cracking is observed experimentally and Leopold & Harder study the relation between ply thickness and the distance between two instances of matrix cracking [START_REF] Leopold | Influence of 90° layer thickness on damage initiation and propagation in CFRP cross ply laminates[END_REF].

Because matrix cracking plays an important role in the occurrence of delamination, the importance of the coupling between interlaminar and intralaminar damage in composite laminates is well accepted in the literature [START_REF] Lopes | Physically-sound simulation of low-velocity impact on fiber reinforced laminates[END_REF][START_REF] Dang | A numerical study on impact and compression after impact behavior of variable angle tow laminates[END_REF][START_REF] Lammerant | Modelling of the interaction between matrix cracks and delaminations during impact of composite plates[END_REF][START_REF] Iarve | Mesh-independent matrix cracking and delamination modeling in laminated composites[END_REF][START_REF] Hallett | Modelling the interaction between matrix cracks and delamination damage in scaled quasi-isotropic specimens[END_REF]. The use of cohesive elements in the ply thickness and between plies allows the coupling between these two types of damage to be captured well and thus gives good results in low velocity impact simulation.

To estimate the residual strength in compression after impact and evaluate the damage tolerance of the laminate, some models simulate both low velocity impact and compression after impact on Abaqus/Explicit [START_REF] Dang | A numerical study on impact and compression after impact behavior of variable angle tow laminates[END_REF][START_REF] Gonzalez | Simulation of dropweight impact and compression after impact tests on composite laminates[END_REF][START_REF] Tan | Predicting low velocity impact damage and Compression-After-Impact (CAI) behavior of composite laminates[END_REF][START_REF] Panettieri | Delamination growth in compression after impact test simulations: Influence of cohesive elements parameters on numerical results[END_REF][START_REF] Caputo | Numerical study of the structural behaviour of impacted composite laminates subjected to compression load[END_REF] and on Abaqus/Implicit [START_REF] Abir | Modelling damage growth in composites subjected to impact and compression after impact[END_REF].

About Multiaxial testing

As detailed by Trellu et al. [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF], only a few test rigs in the world allow structural details to be tested with multiaxial loading. The state of the art concerning this kind of multiaxial test rigs can be found in [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF][START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project. Part 1: Experimental Setup[END_REF]. Two families of machines are described. In the first, actuators directly reproduce the compression/shear loads on a square meter specimen [START_REF] Romeo | Buckling and Post-buckling Behavior of Anisotropic Plates under Combined Biaxial Compression and Shear Loads[END_REF][START_REF] Wolf | An Efficient Test Method for the Experimental Investigation of the Post-buckling Behavior of Curved Shear Panels[END_REF][START_REF] Rouse | Structural Stability of a Stiffened Aluminum Fuselage Panel Subjected to Combined Mechanical and Internal Pressure Loads[END_REF][START_REF] Bergan | Full-scale test and analysis of a PRSEUS fuselage panel to assess damage-containment features[END_REF]. These test rigs look technically complex, with the use of many actuators to obtain compression/shear loading. The second family uses rectangular boxes to apply complex loading to the specimen. In 1984, Peters [START_REF] Peters | Buckling tests of flat rectangular plates under combined shear and longitudinal compression[END_REF] obtained an experimental curve of compression/shear buckling with a square section box stressed in torsion and bending. Klein [START_REF] Klein | General about buckling tests with thin-walled shells[END_REF] used a box structure embedded at one end and loaded with two actuators at the other end. Very recently, Zucco et al. [START_REF] Zucco | Static test of a thermoplastic composite wingbox under shear and bending moment[END_REF] designed and manufactured a testing fixture capable of introducing a prescribed shear force and bending moment at one end of a variable thermoplastic composite wingbox and reacting the load at the other end. This test rig is selfequilibrated and so does not need any joints with the ground (foundations).

Castanié et al. [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF][START_REF] Castanié | Theoretical and experimental analysis of asymmetric sandwich structures[END_REF] developed a multiaxial test rig for asymmetric sandwiches with a box structure. The specimen is bolted to the test rig center on its upper face. Specimen overall dimensions are 248 x 306 mm 2 and the area of interest is 200 x 200 mm 2 . The two actuators at the box extremities allow the structure to be loaded in four-point bending and the specimen in traction or compression. The two actuators in the box center allow the structure to be loaded in torsion and the specimen in shear. This test rig was successfully used to show that asymmetric sandwich structures are extremely strong in compression [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF]. However Castanié et al. pointed out some drawbacks. As for the other test rigs of this type, it is very difficult to estimate the stress flows directly entering in the specimen because of the numerous structural redundancies and the non-linear behavior of asymmetric sandwich structures. Then, the use of small specimens induces Saint-Venant effects, which perturb the strain field. This is why a new experimental test rig called the "VERTEX test rig" was developed [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project. Part 1: Experimental Setup[END_REF][START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project. Part 2: Load envelopes for the assessment of panels with large notches[END_REF]. It is based on the same philosophy as used in Castanié's test and is described in Figure 2. The specimen is now larger (Figure 3-a Testing specimens at the structural detail scale under combined loading is more complex than performing standard tests but detects mechanical phenomena that are much closer to those found in real-life structures. Standards recommended by industrials for low velocity impacts [START_REF] Astm D ; M | Standard Test Method for Measuring the Damage Resistance of a Fiber-Reinforced Polymer Matrix Composite to a Drop-Weight Impact Event[END_REF] and for compression after impact [START_REF] Astm D ; M | Standard Test Method for Compressive Residual Strength Properties of Damaged Polymer Matrix Composite Plates[END_REF] on composite coupons are too conservative. For example, boundary conditions of the compression after impact standard impose a pure compression failure without postbuckling and the specimen size results in rapid propagation of the impact damage to the free edges. In real life, impacted panels are subject to postbuckling and it is interesting to study interactions between impact damage and postbuckling.

Objectives

The Discrete Ply Model (DPM) developed at the Institut Clément Ader by Bouvet et al. makes it possible to simulate low velocity impact and compression after impact, capturing matrix cracking, fiber breakage and delamination [START_REF] Dubary | Damage tolerance of an impacted composite laminate[END_REF][START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF][START_REF] Rivallant | Failure analysis of CFRP laminates subjected to compression after impact: FE simulation using discrete interface elements[END_REF]. The DPM was validated for different stacking sequences with an experiment/model dialog in terms of delamination areas and force displacement curves. It has also shown good results with residual compression strength [START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF][START_REF] Rivallant | Failure analysis of CFRP laminates subjected to compression after impact: FE simulation using discrete interface elements[END_REF] and has been used to optimize a composite laminate design according to low velocity impact damage tolerance [START_REF] Dubary | Damage tolerance of an impacted composite laminate[END_REF]. More recently, simulation of notched behavior in tension has also been achieved with several scales of specimen [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF][START_REF] Serra | Experimental and Numerical Analysis of Carbon Fiber Reinforced Polymer Notched Coupons under Tensile Loading[END_REF].

In the present work, the DPM is used in a "Virtual Testing" approach to simulate medium velocity impacts (50 -110 m/s) and complex loading after impact on large composite laminate plates. This medium speed range was chosen to simulate a new application case with the DPM. To adopt an original position, the laminate scale was chosen at the second level of the pyramid of tests [START_REF] Rouchon | Certification of Large Aircraft Composite Structures, Recent Progress and New Trends in Compliance Philosophy[END_REF]. As presented in the experimental part associated with this research [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF], impacted plates were loaded under complex stresses in order to go further than the recommended standards [START_REF] Astm D ; M | Standard Test Method for Compressive Residual Strength Properties of Damaged Polymer Matrix Composite Plates[END_REF], which seem conservative, and to be able to study interactions between impact damage and postbuckling. The complex loading after impact was obtained with the VERTEX test rig at the Institut Clément Ader [66]. The context of testing composite panels under complex loading with structural test rigs is discussed in detail by Castanié et al. in [START_REF] Castanié | Combined shear/compression structural testing of asymmetric sandwich structures[END_REF][START_REF] Castanié | Theoretical and experimental analysis of asymmetric sandwich structures[END_REF] and the first tests achieved with the VERTEX test rig are presented in [START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project. Part 1: Experimental Setup[END_REF][START_REF] Serra | Validation of aeronautical composite structures under multiaxial loading: the VERTEX Project. Part 2: Load envelopes for the assessment of panels with large notches[END_REF].

A brief reminder is given of the DPM in Section 2, with a description of the numerical modeling used and the experimental campaign of impact and complex loading after impact.

Numerical and experimental results are compared: for medium velocity impact in terms of delaminated area and force displacement curves and for complex loading after impact in terms of loading paths. The aim is to validate the use of the DPM as a tool to study interactions between impact damage and postbuckling.

Numerical model and experimental campaign

The Discrete Ply Model "DPM" is extensively presented in [START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF][START_REF] Rivallant | Failure analysis of CFRP laminates subjected to compression after impact: FE simulation using discrete interface elements[END_REF]. The DPM is run with Abaqus/Explicit and a VUMAT. The model is recalled only briefly here as further details have already been published [START_REF] Serra | Scaling effect in notched composites: the Discrete Ply Model approach[END_REF]. In this section, damage modeling philosophy as well as the material and boundary conditions used, are described. Then the experimental campaign of impact and complex loading after impact is presented.  Delamination is simulated with classical interface elements between two consecutive plies, each ply being modeled with one 3D element in the thickness. This damage is driven using fracture mechanics. Once delamination is initiated, the propagation is modeled by releasing stresses until dissipation of the critical energy release rate.

Discrete Ply Model

 Matrix cracking is modeled with interface elements located between volume elements in the ply direction. This damage is driven using Hashin's Criterion (Eq. 1) calculated in neighboring volume elements.

( ) ( ) (Eq. 1)
with the positive value of the transverse stress, the transverse failure stress, and the shear stresses, and the shear failure stress. When the criterion is reached, element stiffness is set to zero, and the two volume elements become independent. The use of interface elements in the ply thickness imposes a complex mesh but allows the coupling between intra-and inter-laminar damage to be obtained naturally.

 Fiber failure is taken into account using conventional continuum damage and failure mechanics. When the damage initiation strain in traction or in compression is reached, a damage variable corresponding to a linear decrease of the stress is calculated until the final damage strain is obtained, and stresses are determined from the damaged stiffness matrix. is calculated to ensure that the fiber critical energy release rate in mode I, in traction or in compression, is dissipated.

The following equation (Eq. 1) describes the fiber failure criterion.

∫ (∫ ) (Eq. 2)
where V and S are the volume and the section of the element and and the longitudinal stress and strain.

Material and boundary conditions

Medium velocity impact and complex loading after impact were performed on 558 × 536 × and its elastic properties (Table 1) are associated in these shell elements. The use of these elements allowed us to introduce exactly the same boundary conditions as were used with the implicit model [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF] and to decrease computation time. These two different parts were tied with a node to surface constraint [START_REF]Abaqus 2016 Analysis user's manual[END_REF].

The specimen was impacted at its center by a 28 g hemispherical impactor 19 mm in diameter (Figure 5). The boundary conditions during the impact test were given by the contact between the plate and a fixed rigid body, representing the impact window (400 × 400 mm 2 ). After the impact step, a step of relaxation was simulated to eliminate vibrations and waves induced by the impact. During the step of loading after impact, displacements were imposed from DIC (Digital Image Correlation), as presented in a previous article [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF], in order to recreate stress flows entering the specimen (Figure 5). The model was run with Abaqus/Explicit for each step. The step time was set to 30 ms, while the real test of loading after impact lasts about 3 min. Therefore, the DPM capability to simulate quasi-static tests is also assessed.

Experimental campaign

An experimental campaign was performed in order to validate the use of the DPM as a numerical tool, to predict damage occurring during medium velocity impact, and to evaluate the residual strength of the impacted structure under multiaxial loading. Eight specimens were draped. Seven of them (specimens A to G) were impacted at different medium velocities from 54 m/s to 110 m/s and one was kept non-impacted (specimen H). Impacts were performed with the gas launcher LG40 (a 40 mm diameter gas gun) of the impact platform STIMPACT [68]. After impact, all these specimens underwent different loading in the VERTEX test rig. During the impact phase, a speckled pattern was made on the impactor and its trajectory was tracked with a high-speed camera. Images recorded with the high speed camera were processed with a program developed by Passieux et al. [START_REF] Pierré | Finite Element Stereo Digital Image Correlation: framework and mechanical regularization[END_REF] to evaluate the impactor displacement and thus the impact force. The delaminated area after impact was evaluated with ultrasonic investigations. Curves of impact force versus impactor displacement and C-scan for each specimen are compared to numerical results (Figure 6 and The impact velocities and type of post-impact loading selected for each specimen are summarized in Table 2. Experimental and numerical results in terms of delaminated area and curves of impact force versus impactor displacement are shown in Figure 6 (specimens A to C) and Figure 7 (specimens D to G). Concerning the delaminated area, except for specimen A, the extension of the delamination area at 45° at the lower interface is well propagated with the DPM. For other interfaces, the numerical delaminated area is a little greater than the experimental one. This explains why the global delaminated area, except for specimen A, is slightly overestimated even though the correlation between experimental and numerical results remains efficient. Experimental and numerical curves of impact force versus impactor displacement are also similar, except for specimen A, where the numerical impactor displacement is underestimated. This small value of displacement can explain the poor numerical delamination obtained at the lower interface for specimen A. In fact, the delamination propagation at the lower interface is due to the overall bending of the structure. These results demonstrate the capabilities of the DPM to simulate a medium velocity impact on a large composite plate.

Complex loading after impact

In this work, the numerical and experimental results of combined loading after impact were compared for specimens A, B, C, F and G. Therefore, every type of loading was simulated: compression, shear, compression+shear and tension+shear. The pure tension case was not tested because the tension flux needed to break the impacted plate would have been too

close to the limits of the VERTEX multiaxial rig (3000 N/mm). Boundary conditions were applied as displacement on the top surface, as described in Figure 5. For each case, numerical stress/strain curves obtained with the DPM (with and without the simulation of the impact) were compared to stress/strain curves obtained with the experiment. DPM stress/strain curves were obtained with the same method as the experimental curves, as presented in detail in [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF]. The methodology to determine average stresses and global strains from displacement obtained with the help of stereo digital image correlation is briefly recalled in Figure 8. A comparison between numerical and experimental results, in terms of fiber failure scenario,delamination area evolution and displacements was also performed. Compression/shear test (specimen A)

Specimen A, impacted at 54 m/s, was loaded in compression/shear. Stress/strain curves are shown in Figure 9. In compression/shear, curves follow the theoretical and experimental stiffness until the buckling appears. Numerical failure appears earlier than the experimental one when the shear behavior is examined. No difference is visible whether the impact is simulated or not. This is consistent with the very small delaminated areas obtained after impact (numerically and experimentally).

The DPM simulation makes it possible to analyze the failure scenario. The fiber failure is initiated in the bottom right corner and propagated through the impact point to reach the opposite corner (Figure 10). This is very consistent with thermal images recorded during the real test [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF], which show that heat was first emitted in the bottom right corner. The upper ply is oriented at 45°. It was loaded in compression during impact and loaded again in compression during shear loading. It is therefore not surprising to observe it failing first in compression fiber failure mode.

The DPM shows that the delaminated area after impact is still the same as when the first failure occurred and does not evolve during the failure propagation (Figure 11). Only a delamination propagation according to fiber failures in the top ply appears. failure is not well reproduced by the models and therefore the numerical failure stress is higher than expected. As for specimen A, no difference is visible whether the impact is simulated or not. This is consistent with the very small delaminated areas obtained after impact (numerically and experimentally). Concerning the failure profile (Figure 14) with the DPM simulation, the failure is initiated in the bottom right corner and propagated through the impact point to reach the opposite corner. This is in contradiction with images recorded during the real test, which seem to show that failure is initiated in the impact zone, and then propagates to the bottom right and the top left corners. This may be explained by a poor quality of the boundary conditions associated with a poor displacement correlation due to several specular reflections of light on the speckle pattern (Figure 16). The discrepancy may also stem from an insufficient simulated fiber failure on the impacted side, since this probably initiates the global failure.

The DPM shows that the delaminated area after impact is still the same as when the first failure occurred and does not evolve during the failure propagation (Figure 15). Only a delamination propagation according to fiber failures in the top ply appears. Figure 16 compares displacement fields between experimental results and numerical results. The displacement fields are correctly described but the amplitudes are slightly overestimated (grey and black zones). In other words, the post-buckling behavior is overestimated. This corroborates the observations made on Figure 13. Specimen C (impacted at 75 m/s) was subjected to tension/shear loading. The stress/strain curves (Figure 17) allow the tensile and shear behaviors to be analyzed. In tension and in shear, the numerical stiffnesses are very similar to both the theoretical and the experimental values until the buckling point is reached. The correlation of the buckling point is very good for the shear behavior but there is a small difference in the tension behavior: the buckling point is underestimated. After the buckling, the global strain value decreases as the experimental value and the stress continues to increase. The final failure occurs around the same stress value but the numerical failure longitudinal strain (-500 µdef) is notably lower than the real failure strain (-200 µdef). This difference is due to underestimation of the buckling point by the numerical models: the curves are shifted to the left of the figure but follow the same trend. Of course, this discrepancy should be considered with caution, given the very low value of the imposed longitudinal strain. In shear after buckling, the behavior is well simulated by the numerical models. Simulating the impact does not have a substantial influence here; only the final failure appears earlier (tensile behavior) when it is taken into account.

The DPM final failure profile (Figure 18) is close to the experimental failure. The numerical failure starts from the initial delaminated area due to impact and then propagates to both the top left and bottom right corners. The experimental data do not make it possible to confirm this scenario precisely because the damage propagates very suddenly. The infrared camera was not used at a high enough frame rate to capture its propagation in both directions. As for specimen B (in pure shear), the delaminated area after impact is not modified by the tension/shear loading (Figure 19). Only delamination according to the fiber failure in the top ply propagates. Figure 20 shows adequate correlation between DPM and experimental displacement fields before failure. The displacement fields are correctly described but the amplitudes are slightly overestimated (grey and black zones). In other words, the post-buckling behavior is overestimated. This corroborates the observations made on Figure 17. Specimen F, impacted at 100 m/s, was loaded in tension/shear. Stress/strain curves in tension and shear are shown in Figure 21. The oscillations observed on the experimental curves are simply due to the use of a higher number of data points during the testing of specimen F, which exhibit the servo-control of the multiaxial test-rig. In tension, curves follow the theoretical and experimental stiffness until the buckling appears. Numerically determined failure stress is lower than that found by experiment. Regarding the shear behavior, the initial stiffness is slightly overestimated by the simulations. The buckling point, and therefore the failure stress, is markedly underestimated by the numerical models. The behavior in tension for this specimen is better than for specimen C, which was also loaded in tension/shear but impacted at 75 m/s. The opposite is, however, true for the behavior in shear. No difference is detected whether the impact is simulated or not.

The numerical final failure profile is similar to the experimental failure profile. The fiber failure initiates in the bottom right corner (Figure 22) as observations made with the thermal camera in [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF] suggested. Afterwards, fiber failure occurs in the top left corner and later in the central damaged area due to preliminary impact. Eventually, fiber failures coalesce in a global -45° orientated line. As for specimen C, also loaded in tension/shear, the failure scenario of specimen F is influenced by the initial damage. As for specimens A, B and C, the delaminated area after impact follows the fiber failure scenario exactly (Figure 23).

In Figure 24, DPM numerical and experimental displacement fields are compared.

Numerically obtained displacements are comparable to those extracted from the experiment. The specular reflections due to a poor choice of paint for the speckle pattern (glossy instead of mat) led to zones of the images where the correlation was impossible to estimate (bottom left corner). Therefore, incomplete displacement data led to incomplete boundary conditions and wrong simulations of the displacement in this area. Specimen G, impacted at 110 m/s, was loaded in compression and the corresponding stress/strain curves are shown in Figure 25. The curves follow the theoretical and experimental stiffness until buckling appears. The first and second buckling modes are correctly predicted by the models (Figure 25). After the second buckling, the numerically obtained behavior diverges slightly form the experimental one. This might have to do with the complex evolutions of the displacement fields due to two buckling modes and the way averaged stresses are obtained [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF]. A small difference in the buckling behavior can translate into large difference in the averaged stress. No difference is seen whether the impact is simulated or not. Concerning the failure profile (Figure 26), the DPM simulation is consistent with the experiment since the failure is initiated in the bottom right corner and propagates towards the center of the plate. The numerically determined orientation is, however, different from that observed experimentally.

The DPM shows that the delaminated area after impact is still the same as when the first failure occurred and does not evolve during the failure propagation (Figure 27). Only a delamination propagation according to fiber failures in the top ply appears. A comparison of displacement fields between experimental and numerical results is shown in Figure 28. The displacement fields are correctly described. With the exception of specimens C and F in the global strains graph, the loading path of each specimen up to final failure is very well represented with the DPM. The DPM is efficient to capture the buckling and the associated buckling modes but sometimes leads to an earlier failure. The DPM shows its capabilities as far as studying the interaction between the delaminated area after impact and the structure postbuckling is concerned. As estimated in [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF], damage induced by the impact has little influence on the post-buckling failure except in specimen F. Final failure profiles are also well modeled and the DPM proves to be useful to study the failure initiation and propagation.

Conclusion

In this paper, the DPM was used to simulate VERTEX tests described in [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF]. The objective of DPM simulations is to be able to answer a few questions regarding the location of the failure initiation and the interaction between impact damage and postbuckling behavior.

The first step was to validate DPM capabilities to simulate medium velocity impact, from 54 m/s to 110 m/s, on large composite laminate plates. Thanks to a numerical/experimental comparison of force versus impactor displacement curves and delaminated areas obtained with C-scan, this first step was successfully validated.

The second step was to simulate complex loading after the impact. Boundary conditions recovered from DIC were applied to the DPM with the method validated using an implicit model without any damage representation [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF]. In order to implement these boundary conditions, the plate was meshed in three parts: one with damageable volume elements used with DPM behavior laws, one with volume elements (1 per ply) used as a transition zone and one external part with one non-damageable thick shell element in the plate thickness. This configuration allowed the authors to use exactly the same boundary conditions as those used in the simplified mode in [START_REF] Trellu | Combined loadings after medium velocity impact on large CFRP laminate plates: Tests and enhanced computation/testing dialogue[END_REF]. Five specimens were studied in this work: specimen A in compression/shear loading (impacted at 54m/s), specimen B in pure shear loading (impacted at 70 m/s), specimens C and F in tension/shear loading (respectively impacted at 75 m/s and 100 m/s) and specimen G in pure compression (impacted at 110 m/s).

DPM simulations gave good results in all types of loading. Tension/shear behaviors could be improved since the buckling points appeared too early in the simulations. The post-buckling behaviors were therefore influenced and diverged slightly from the experimental ones. The discrepancy between numerical and experimental results can probably be explained by the quality of the boundary conditions extracted from the experimental data (from stereocorrelation) and implemented in the model. For the present study, the shape of the plate once bolted to the test bench was not taken into account in the numerical model. In other words, the experimental boundary conditions were applied to a perfectly flat plate whereas the experimental observations show that, after this bolting operation, some loading is already applied to the specimen, which deforms accordingly. This initial small curvature should have some influence on the global behavior, especially when loading (compression or shear) triggers buckling. Additionally, it should be mentioned that displacements used as boundary conditions were applied to both the top and the bottom surface of the plate. To be efficient, it would be necessary obtain the rotation of the plate through different displacements applied to the top and the bottom surfaces. Research is still on-going to assess the influence that the type of boundary conditions used in a numerical model has on the correlation with experimental data.

In most cases, the DPM confirms the location of the failure occurrence observed experimentally. Except for specimen B, the failure appears first in the plate bottom right corner and then propagates to the top left corner through the impact mark. This also makes it possible to show that the delaminated area after impact was not propagated during the loading after impact. This behavior seems to indicate that the impact damage does not play a predominant role in the final failure.

The DPM gave encouraging results for simulating the medium velocity impact and the complex loading after impact in the same model. The implementation of boundary conditions in an explicit scheme must be investigated if a better correlation between numerical and experimental results is to be obtained. Finally, computation times associated with such simulations are too long (50 h with 36 CPUs). This last point needs to be improved before industrial use of the DPM can be envisaged.

Damage from the impact had little influence on the post-buckling failure except for specimen B. This result is very interesting because it contradicts the classical results of the literature [START_REF] Tropis | Certification of the composite outer wing of the ATR72[END_REF][START_REF] Caprino | The significance of indentation in the inspection of carbon fibrereinforced plastic panels damaged by low-velocity impact[END_REF]. It is noteworthy that the energy involved in the present study (from 40.8 J to 169.4 J) is very high compared with the classical ones (around 50 J). These conflicting results should be explained by the size of the specimen. The standard tests of impact and compression after impact [START_REF] Astm D ; M | Standard Test Method for Measuring the Damage Resistance of a Fiber-Reinforced Polymer Matrix Composite to a Drop-Weight Impact Event[END_REF][START_REF] Astm D ; M | Standard Test Method for Compressive Residual Strength Properties of Damaged Polymer Matrix Composite Plates[END_REF] are performed on small coupons of 100 × 150 mm 2 ; which induces greater impact damage (for a given impact energy) and, at the same time, a larger zone of impact damage relatively to the coupon size. The present study shows that, for larger specimens, the impact damage is less prevalent and decreases the residual strength less. In fact, the buckling mode clearly affects the final failure more than the impact damage does. This result stresses the interest of the pyramid of tests design philosophy [START_REF] Rouchon | Certification of Large Aircraft Composite Structures, Recent Progress and New Trends in Compliance Philosophy[END_REF], even though its benefit is already very clear in the industry.

In conclusion, the pyramid of tests design is clearly a crucial point of composite structure design, and test rigs, such as the VERTEX machine [66], able to evaluate the residual strength after impact on a real structure, or at least at the level 2 of the pyramid of tests [START_REF] Rouchon | Certification of Large Aircraft Composite Structures, Recent Progress and New Trends in Compliance Philosophy[END_REF] are, more than ever, needed to avoid oversizing composite structures.

Figure 1 :

 1 Figure 1: Damage tolerance concept with detectable and undetectable damage [2]

  ) and is bolted to the structure with 128 fasteners, 8 and 10 mm in diameter (Figure3-b). The VERTEX test rig was used in the present study to carry out complex loading tests.

Figure 2 :Figure 3 :

 23 Figure 2: VERTEX test rig

Figure 4

 4 Figure 4 presents the concept of the DPM. The main failure modes occurring in an impacted

Figure 4 :

 4 Figure 4: DPM concept: element types and associated damage [64]

Figure 5

 5 Figure 5 shows the mesh and numerical boundary conditions used in the model. In this

Figure 5 :

 5 Figure 5: Boundary conditions of the finite element model

Figure 7 .

 7 Figure 7. For post-impact loading in the VERTEX test rig, a pair of 5 Mpx cameras recorded

Figure 6 :Figure 7 :

 67 Figure 6: Delaminated areas and curves of impact force versus impactor displacement for

Figure 8 :

 8 Figure 8: Method to determine average stresses and global strains [51]

Figure 12 comparesFigure 9 :Figure 10 :Figure 11 :Figure 12 :

 129101112 Figure 12 compares displacement fields between experimental results and numerical results

Figure 13 :

 13 Figure 13: Stress/strain curves for specimen B in shear loading

Figure 14 :Figure 15 :

 1415 Figure 14: Comparison between DPM and experimental failure profile for specimen B

Figure 17 :Figure 18 :Figure 19 :Figure 20 :

 17181920 Figure 17: Stress/strain curves for specimen C in tension/shear loading

Figure 21 :Figure 23 :Figure 24 :

 212324 Figure 21: Stress/strain curves for specimen F in tension/shear loading

Figure 25 :

 25 Figure 25: Stress/Strain curves for specimen G in Compression loading

Figure 26 :Figure 27 :SynthesisFigure 29

 262729 Figure 26: Comparison between DPM and experimental failure profile for specimen G

  Figure 29: Comparison

  3.5 mm 3 laminate plates of T700/M21 UD carbon/epoxy composite. The stacking sequence chosen for this model was quasi-isotropic, with 14 plies [45 2 /-45 2 /0 2 /90 2 /0 2 /-45 2 /45 2 ]. The ply thickness was 0.25 mm. Material properties used in the DPM are presented in Table 1.

Table 1 . T700/M21 properties used in the model

 1

	Elastic Properties
	Tensile Young's modulus in fiber direction	130 GPa
	Compressive Young's modulus in fiber direction	100 GPa
	Transverse Young's modulus	7.7 GPa
	Poisson ratio	0.3
	Shear modulus	4.75 GPa
	Delamination & matrix cracking
	Transverse tensile strength	60 MPa
	In-plane shear strength	110 MPa
	Interface fracture toughness for opening mode (I)	0.5 N/mm
	Interface fracture toughness for shear mode (II &	1.6 N/mm
	III)
	Fiber Failure
	Fiber tensile strain at damage initiation	1.70 %
	Fiber compressive strain at damage initiation	1.40 %
	Fiber tensile fracture toughness	90 N/mm
	Fiber compression fracture toughness	30 N/mm

Table 2 . Overview of impact velocities and post-impact loadings

 2

		Velocity	Impact Energy
	Specimen			Loading after impact
		(m/s)	(J)
	A	54	40.8	Compression/shear
	B	70	68.6	Shear
	C	75	78.7	Shear/tension
	D	90	113.4	Compression/shear
	E	98	134.4	Compression
	F	100	140	Shear/tension
	G	110	169.4	Compression
	H	0	0	Shear/tension

Acknowledgements

The authors gratefully acknowledge CALMIP (CALcul en Midi-Pyrénées, (https://calmip.univtoulouse.fr) for access to the HPC resources under the allocation p1026.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:


Joud Loïc
email: loic.joud@danielson-eng.fr

Daniela Chrenko
email: daniela.chrenko@utbm.fr

Alan Keromnes
email: alan.keromnes@u-bourgogne.fr

Rui Da

Luis Le Moyne
email: luis.le-moyne@u-bourgogne.fr

Predictive Energy Management for Hybrid Vehicles based on Driving Cycle Recognition

Keywords: Plug-In Hybrid Vehicle, Series Hybrid Vehicle, Energy management, Cycle Recognitions, Dynamic Programming

This article presents a novel approach to design an optimum energy management for a series plug-in hybrid electric vehicle based on driving cycle recognition. Thanks to an intensive study of driver habits, an important part of cycles is predictable. Optimized energy management can be developed for representative cycles via DP. Those results can be applied to the same cycle with slight adaptations thereafter with quality factors above 98%.

I. INTRODUCTION

Facing limitation of resources and global warming, users are forced to find more environmentally sensible, safe and cost effective solutions for their personal transport. Electric vehicles do not yet meet the demands due to limited battery lifetime, long recharge times, lacking recharge infrastructure [START_REF] Devloo | Component Improvements in the Electrification of Passenger Vehicles Drivetrains[END_REF][2] and high cost. However, Plug-in Hybrid Electric Vehicles (PHEV) can be the missing link from conventional to electric vehicles as they can be recharged by both liquid fuel and electricity and their electric autonomy of 30 to 50 km allows to cover at least 60% of driving cycles [START_REF] Wirasingha | Plug-in hybrid electric vehicle developments in the US: Trends, barriers, and economic feasibility[END_REF][4] [START_REF] Ericsson | Variability in urban driving patterns[END_REF] [START_REF] Khan | Predicting the market potential of plug-in electric vehicles using multiday GPS data[END_REF]. In order to benefit from their full potential, it is however required that the PHEV arrives at the destination just at the moment when its battery is empty and that electrical energy from renewable sources is available at destination [START_REF] Mehar | EV-planning: Electric vehicle itinerary planning[END_REF] [START_REF] Gong | Trip Based Power Management of Plug-in Hybrid Electric Vehicle with Two-Scale Dynamic Programming[END_REF]. If thermal liquid fuel is used, whereas the cycle could have been covered by electric energy, nonnecessary pollution is created; if the battery is empty before arriving at destination, the internal combustion engine has to cover all working points, reducing the global efficiency [START_REF] Wirasingha | Classification and Review of Control Strategies for Plug-In Hybrid Electric Vehicles[END_REF] [START_REF] Larsson | Commuter Route Optimized Energy Management of Hybrid Electric Vehicles[END_REF].

Wirasingha and Emadi [START_REF] Wirasingha | Classification and Review of Control Strategies for Plug-In Hybrid Electric Vehicles[END_REF] present an exhaustive analysis of control strategies for PHEVs, but it does not cover aspects of prediction in detail.

The purpose of the present work is to know if a particular effort on the energy management on usual cycles and driver's habits (journeys, speed profile) is meaningful. If it is the case, the goal is to develop a tool in order to optimise the energy flux inside a series hybrid powertrain. This is done by suggesting a management, which identifies usual cycles and which is able to adapt the strategy according to the results of the usual cycle recognition with the goal to reach the best energy distribution along the journey. The work will give rise to an innovative predictive management. Different approaches have been presented [START_REF] Filev | Contextual On-Board Learning and Prediction of Vehicle Destinations[END_REF], [START_REF] Ye | Predict Personal Continuous Route[END_REF].

The article is structured as follows: An overview of vehicle use in France is provided in Section II, an analysis of vehicle use in the Nièvre Department is presented in Section III. Section IV presents the energy management including cycle prediction. Conclusions and Perspectives close the article in Section V.

II. OVERVIEW OF PERSONAL VEHICLE USE IN FRANCE

In order to define the best strategy, it is first necessary to verify that the use of the vehicle is driven by habits. Figure 1 presents the main displacement reasons in France in 1994 and 2008 [START_REF]Commissariat Général au Développement Durable[END_REF]. It can be noticed that more than 50% of the journeys are constrained (work, shopping and study). This distribution is stable over the last decades. Therefore, it might be possible to identify the destinations during working days. Moreover, the average distance of each journey is short (less than 30 km, Table 1) and could be covered electrically by the PHEV. However, these average values do not give any information about the influence of driver behavior (speed profile, acceleration), which depends also on the car's characteristics (power, weight,…) and will affect the energy management. Even if it is possible to create realistic driving cycles thanks to modeling, it can be interesting to study the real use of vehicles in order to identify if a limited number of cycles is regularly reproduced. Therefore, in order to determine these recurrent cycles and the influence of driving behavior, as well as to verify the previous values in the rural context of the Nièvre

: Mobility Evolution of Local displacement reasons on working days in

France for people older than 6 years old) [START_REF]Commissariat Général au Développement Durable[END_REF] department in France, it has been decided to perform a mobility study. The approach is presented in the nest section.

III. ANALYSIS OF VEHICLE USE IN NIEVRE DEPARTMENT

A. Measurement Approach to Evaluate Vehicle Use

In order to be able to get significant results which correspond to previous study, it has been necessary to define a panel, which corresponds to the quota method (age, sex) inside the working population of a rural area. In Nièvre, the working population is composed by 47% of women and 53% of men and the age distribution is presented in Figure 2. It can be noticed that the working population is equally distributed according to the age criterion between 20 and 60. Only a few percentage of the workers are older than 60 due to retirement (legal retire age: 62 years). In order to fulfill the criterion a panel of 40 volunteers has been determined and the test were performed over 4 months.

In order to perform the mobility study, ten identical autonomous data acquisition systems have been developed based on Arduino systems. These systems are placed inside the vehicle during one month and do not require any action from the driver. They are composed by a 9 axis inertial unit, a GPS (speed and position) and a 32 Go micro SD card in order to store the data. The power unit is based on 18 Lithium-ion batteries. Data is acquired at a frequency of 10 Hz.

B. Results

Table 1 : Comparison between Nièvre against mean of France [START_REF]Commissariat Général au Développement Durable[END_REF] During the mobility study of 4 months including 40 drivers, a total of 1625 journeys have been performed, which corresponds to over 13,000 km. The results are compared to previous study from [START_REF]Commissariat Général au Développement Durable[END_REF] in Table 1 on a daily basis in rural areas. It can be noticed that the numbers of journeys and the average daily distance are very similar. However, in the present study the average time is shorter due to a higher average speed. The difference is due to a motorway, which is often used by drivers. Detailed results are presented in Figure 3. As mentioned previously, it can be noticed that most drivers travel between 20 and 30 km per day. However, some drivers are located further away from work and travel up 55 km per day. In order to optimize the energy distribution, the first step is to identify if some usual cycles can be determined for each driver. So, a cross-correlation factor is calculated for each driver considering daily repetition. Based on this, it is possible to quantify the number of occurrences of the different trips.

For each starting position, the number of identical journeys is compared to the total number of journeys starting from this point. The identical journeys are determined thanks to a normalized cross-correlation factor, σ (equation 1) and the journeys are considered to be identical if σ is greater than 0.9. The predictability is defined as the ratio between the cycles with σ greater than 0.9 and the total number of cycles starting at this position. The results (Figure 4) indicate that, on working days, more than 75% of the trips are performed more than five times. Therefore, it is possible to identify similarities during working days, as presented in Table 2 for all the database and a selection of drivers, since the correlation factor is greater than 82% during working days. However, drivers' habits are more difficult to identify during week-ends since the overall number of trips is reduced and the number of repetition of the same trip is reduced: 80% of the trips are performed less than five times. Moreover, most journeys are performed during working days since they represent more than 80% of the journeys (1304 trips for a total trip number of 1625). However, some drivers, such as the woman between 40 and 50 years in Table 2, have a limited journey range and a more important number of journeys with a higher number of destinations. In this case the correlation factor is reduced to 60%. A further analysis of this driver indicates that most journeys are performed in urban areas.

= ∑ , , , ̅ , ∑ , , , ∑ , ̅ , (1)

IV. ENERGY MANAGEMENT INCLUDING CYCLE PREDICTION

The energy management is based on cycle recognition including speed profile. The selected speed cycle is then used inside a vehicle model to determine the energy and power needs during the cycle and optimize the energy management.

A. Cycle prediction

Since it is possible to identify a repetitive scheme, it is necessary to store and sort the journeys as a function of destination. This will allow the identification of the potential destination as a function of origin and day as well as the determination of the best energy distribution over the journey. Thus, it is necessary to determine a representative speed cycle for a journey. The use of an average speed cycle is not possible since it would completely erase the driver behavior and the speed fluctuations. Therefore, it has been decided to determine among the stored speed cycles for each journey the most representative one based on a cross-correlation criteria calculated from the speed profile. The equation is similar to equation 1. This speed cycle is used to calculate the most efficient energy distribution. Each time a journey is performed, the new speed cycle is stored and compared to previously stored cycles using the cross-correlation. The energy management system only keeps the ten cycles with the highest crosscorrelation criteria. This reduces the impact of chaotic traffic issues (traffic jams, accident and road work) on the speed cycles.

B. Vehicle model

The vehicle model is based on Energetic Macroscopic Representation (EMR) and its Inversion Based Control (IBC) of a plug-in series hybrid vehicle [START_REF] Jin | The control strategy and cost analysis for series Plug-in hybrid electric vehicle[END_REF]. The model has been previously validated over the NOAO racing car which was developed and tested in Magny-Cours [START_REF] Asus | Model and Control Strategy Simulation of a Racing Series Hybrid Car[END_REF]. This model includes all the different thermal and electrical components of the powertrain and is able to calculate the mass flow rate of fuel and the state-of-charge (SOC) of the battery. The energy and power needs are calculated based on the effort (aerodynamics, inertia, friction) applied on the vehicle and the respect of the speed cycle. Based on the rules used for the energy management in the NOAO, the model reproduces the power profile from the electric machine and the fuel consumption of the internal combustion engine with a very good accuracy. Finally, the maximum difference between the experimental and numerical SOC is below 1%

C. Optimisation strategy

It is commonly admitted that Dynamic Programing (DP) is able to provide the best optimization strategy of a system [START_REF] Wirasingha | Classification and Review of Control Strategies for Plug-In Hybrid Electric Vehicles[END_REF], [START_REF] Chrenko | Novel Classification of Control Strategies for Hybrid Electric Vehicles[END_REF]. However, the calculation time makes it impossible to use on-line. Therefore, Dynamic Programming is used off-line at the end of the journey in order to determine which would have been the best energy management along the journey taking into account the state-of charge of the battery and the fuel consumption. In order to verify the potential gain thanks to DP, it has been used in order to calculate the best energy management and the resulting fuel consumption for each cycle. The results are presented for a selection of five similar cycles of the same journey in Table 3. For those similar cycles, the energy management determined using DP for the most representative cycle (i.e. cycle 5 determined by crosscorrelation) has been applied for the other occurrences of this cycle. This results in an increase of the fuel consumption by up to 7% according to the cycle when compared to the DP calculation of the cycle. Therefore, a blind application of DP is not the best strategy. A second optimization strategy, based on least square minimization, has been implemented in order to adapt the energy management online during the journey. The method minimizes the difference between (i) the optimal consumption calculated thanks to DP on the most representative cycle and (ii) the consumption obtained with the energy repartition chosen at time t. The algorithm also calculates the energy management over a foreseeable window. An optimality factor has been determined in order to compare the fuel consumption determined by the strategy with the optimum fuel consumption calculated off-line thanks to dynamic programming. This factor corresponds to the fuel consumption predicted by DP divided by the consumption obtained thanks to the optimization. As an example, when applying this new optimization strategy to cycle number 1, it has been possible to reach a fuel consumption of 0,742 L. This represents an optimality factor of 0.983. The extension of this strategy to the experimental database is presented in Figure 5. It is possible to reach an optimality factor greater than 0.95 for all the speed cycles of the database recorded during week days (Figure 5).

V. C ONCLUSION AND PERSPECTIVES

The present study aims at developing an optimization tool in order to obtain a predictive energy management inside a plug-in series hybrid vehicle. An experimental investigation of mobility revealed the more than 82% of the displacement can be predicted during working days due to constrained mobility. Therefore, the optimization strategy could be applied in order to reduce the fuel consumption thanks to the prediction of the journey based on a cross correlation calculation.

The optimization strategy has been compared to the energy management calculated thanks to dynamic programming using an energetic macroscopic representation of the vehicle. The developed optimization strategy allows to reach at least 95% of the optimality predicted by dynamic programing.

In the future, the prediction of the speed cycle will be made more precise since it will take into account the starting time. Moreover, in order to improve the prediction, the instantaneous location will be compared to stored starting position in order to calculate the optimization strategy based on a larger database and include more speed cycle. This will also allow a better prediction and a reduced fuel consumption when the journey is performed for the first time.

Figure 1

 1 Figure 1: Mobility Evolution of Local displacement reasons on working days inFrance for people older than 6 years old)[START_REF]Commissariat Général au Développement Durable[END_REF]

Figure 3 :

 3 Figure 3: Average daily distance for the drivers in the present study

Figure 4 :

 4 Figure 4: Repetition of trips

Figure 2 :

 2 Figure 2: Selection of tested drivers by age

Figure 5 :

 5 Figure 5: Distribution of the optimality reached thanks to the optimization

Table 2

 2

	: Comparison of driving habits of 4 drivers during working days and
		week-ends
		Number of Journeys	Average Distance	Predictability
		Working days
	W 20/30	53	15 km	87%
	M 30/40	70	9,9 km	90%
	W 40/50	80	7,5 km	60%
	M 50/60	35	5,6 km	80%
	overall	1304	8,6 km	82%

Table 3 :

 3 Impact of DP strategy on fuel consumption for 9.2 km journey with an average speed of 44 km/h for a starting SOC of 0.54 and finishing at 0.37

	Cycle	1	2	3	4	5
	Correlation σ (speed)	0,983 0,979 0,976 0,981	1
	DP Fuel consumption (L)	0,73	0,64	0,61	0,73	0,70
	Cycle 5 strategy (L)	0,79	0,67	0,65	0,77	0,70
	Gap (%)	7,5	4,4	6,1	5,2	/

ACKNOWLEDGEMENTS

The authors want to acknowledge all the volunteers who participate to the mobility investigation and the Burgundy Council for the financial support via ID-Motion laboratory.