Daniela Yassuda Yamashita

Ionel Vechiu

Jean Paul Gaubert

Samuel Jupin

Autonomous observer of hydrogen storage to enhance a model predictive control structure for building microgrids

Keywords:

Hydrogen energy storage has emerged as a promising technology to improve the integration of renewable energy sources in building microgrids. However, inaccuracies in the modelling of fuel cells and electrolysers reduce the performance of building microgrids' energy management system. To improve the flexibility of building microgrids, this paper proposes to associate a two-level hierarchical model predictive controller empowered with an Autonomous Observer of Hydrogen Storage (AOHS). This novel observer evaluates the hydrogen production and consumption rates, storing little past data and needing no tuning of the parameters. Relying only on instantaneous data measurement, the algorithm can estimate the tank's level of hydrogen with an average relative error inferior to 2 %, even under measurement noise. A case-study based on a building microgrid currently under construction serves as the basis for all simulations. The performance of the AOHS is evaluated by comparing the self-consumption rates of the case-study when governed by two-level energy management system: one level using a fixed parameters model and the other one equipped with the proposed AOHS algorithm. Results show that the microgrid associated to the AOHS has better self-consumption compared to the microgrid with fixed parameters, as well as a better robustness regarding the measurement noise and modelling error. Furthermore, this algorithm demonstrates a planning function as it facilitates the energy planning from the aggregator's point of view and the external grid management.

Introduction

Even though energy is essential for society, the energy sector is built on an unsustainable system that meets about 70 % of the primary energy demand with fossil fuel and nuclear power plants [START_REF]Perspective for the Clean Energy Transition 2019 -The Critical Role of Buildings[END_REF]. To foster the decarbonisation of the electricity systems, the fast energy transitioning scenario aims to integrate more than 5 thousand TWh of Photovoltaic (PV) panels into the grid by 2050 [START_REF]Perspective for the Clean Energy Transition 2019 -The Critical Role of Buildings[END_REF]. In this context, building microgrids, such as the ones described in [START_REF] Hu | Neighborhood-level coordination and negotiation techniques for managing demand-side flexibility in residential microgrids[END_REF][START_REF] Wu | Towards collective energy community: potential roles of microgrid and blockchain to go beyond P2P energy trading[END_REF][START_REF] Fontenot | Modeling and control of building-integrated microgrids for optimal energy managementa review[END_REF][START_REF] Farinis | Integrated energy management system for microgrids of building prosumers[END_REF], have emerged as a promising grid topology to enable the massive installation of renewable energy sources directly in the distribution sector.

To avoid radical changes in the electricity system when incorporating building microgrids into the current power grid, the simplification of the role of market operators is paramount. This simplification is performed by introducing the concept of self-consumption of PV electricity [START_REF] Masson | Comprehensive Renewable Energy[END_REF]. This grid regulation aims to minimise the building's microgrid daily power imbalance, allowing grid operators to supply the energy demand as efficiently as possible [START_REF] Gomez | European Union electricity markets: current practice and future view[END_REF].

To attain high levels of self-consumption, energy storage systems (ESSs) are key elements to shift the building local energy surplus toward periods of energy deficit [START_REF] Bartels | Influence of hydrogen on grid investments for smart microgrids[END_REF]. Particularly, due to the seasonality throughout the year, the combination of short-and long-term ESSs is increasingly envisaged for building microgrid implementation [START_REF] Ren | Optimal operation of a grid-connected hybrid PV/fuel cell/battery energy system for residential applications[END_REF]. In recent years, hydrogen-based storage systems have become important seasonal ESSs, thanks to their high energy capacity (300-1200 Wh/kg) and their nearly zero self-discharging rates [START_REF] Lokar | The potential for integration of hydrogen for complete energy self-sufficiency in residential buildings with photovoltaic and battery storage systems[END_REF][START_REF] Maestre | Challenges and prospects of renewable hydrogenbased strategies for full decarbonization of stationary power applications[END_REF][START_REF] Argyrou | Energy storage for electricity generation and related processes: technologies appraisal and grid scale applications[END_REF].

To take full advantages of hybrid ESSs while maximising PV selfconsumption, it is necessary to design a proper building energy management system (EMS) [START_REF] Fontenot | Modeling and control of building-integrated microgrids for optimal energy managementa review[END_REF]. Nonetheless, buildings are not homogeneous and require solutions tailored to their specific conditions. Therefore, this paper aims to propose a flexible solution to adapt EMSs to any hydrogenbased building microgrid automatically and independently of the electrical equipment size.

Among the existent EMS algorithms [START_REF] Yamashita | A review of hierarchical control for building microgrids[END_REF][START_REF] Zia | Microgrids energy management systems: a critical review on methods, solutions, and prospects[END_REF], model predictive control (MPC) has proved its robustness against environmental disturbances [START_REF] Valverde | Integration of fuel cell technologies in renewableenergy-based microgrids optimizing operational costs and durability[END_REF][START_REF] Hu | Model predictive control of microgridsan overview[END_REF][START_REF] Kamal | Model predictive control and optimization of networked microgrids[END_REF]. The capacity to consider prediction data and periodic optimisations over a sliding window are the main strengths of the MPC * structure, making it highly appreciated for industrial applications. However, the performance of the MPC structure strongly depends on the quality of the predictions and on the models used.

In the context of long-term power flow optimisation, as it occurs in energy management for building microgrids, MPC requires precise ESS models. However, in the literature, MPC architectures are usually composed of time-invariant parameters derived from the manufacturer's technical specifications, therefore holding inaccuracies concerning the real system [START_REF] Garcia-Torres | Optimal economical schedule of hydrogen-based microgrids with hybrid storage using model predictive control[END_REF]. To improve the model precision, a model can be devised from the linearisation of non-linear experimental curves [START_REF] Valverde | Integration of fuel cell technologies in renewableenergy-based microgrids optimizing operational costs and durability[END_REF] or by fitting sample measurements to a logarithmic equation [START_REF] García-Valverde | Simple PEM water electrolyser model and experimental validation[END_REF]. Alternatively, advanced techniques to estimate the intrinsic parameters of proton exchange membrane (PEM) fuel cell and PEM electrolyser result in a precise estimation of the tank level-of-hydrogen. For instance, models based on the physical structure of the PEM fuel cell lead to an accurate model [START_REF] Abdin | PEM fuel cell model and simulation in Matlab-Simulink based on physical parameters[END_REF]. Similarly, Kalman filters and sliding mode observers [START_REF] Yuan | Model-based observers for internal states estimation and control of proton exchange membrane fuel cell system: a review[END_REF] anticipate fuel cell model changes such as oxygen and hydrogen partial pressure and stack temperature. However, these models are too complex to be embedded in MPC prediction as it requires many parameters that are usually hard to measure. Furthermore, these methods are usually time-invariant which may be inaccurate for longterm operation, especially due to the ageing of components and unexpected changes in the operating environment.

Aiming to reduce the effects of both imprecision in model-based state estimation and stochasticity in the building microgrid power imbalance without raising the complexity of the MPC algorithm, this paper proposes the association of a two-level hierarchical MPC structure with an Autonomous Observer of Hydrogen Storage (AOHS). The primary objective of the proposed EMS is to maximise the annual selfconsumption rate by optimising the building microgrid power flow. The innovation of the proposed EMS resides in the use of the AOHS and its direct inclusion in the MPC algorithm. This way, no parameter needs to be tuned to set up the MPC framework. Thanks to online data processing embedded in the AOHS, the EMS is capable of identifying the hydrogen ESS model parameters even under noise measurements. Moreover, the proposed AOHS algorithm achieves a highly precise modelling performance without storing past measurements. Consequently, the identification of the hydrogen ESS model is implemented only using a few stored variables.

The remainder of this paper is organised as follows. Section 2 presents the context of the building microgrid and this paper's case-study. Section 3 describes the hierarchical MPC structure by highlighting its cost function and constraints. Section 4 details the hydrogen ESS autonomous observer. Section 5 compares the performance of an MPC architecture associated with the AOHS algorithm with a classical MPC architecture without an AOHS. Finally, Section 6 concludes this paper.

The context of the building microgrid

The building microgrid selected as a case study is the future grid-

Abbreviations

AOHS

i ∈ {1, 2, 3} (A W -1 ) M k
Measure point at instant k G zone i x X-coordinate of the median of the three-zone power lines σ PEM, k Approximate angular coefficient linking the current and power consumed/generated by the electrolyser or fuel cell connected microgrid that will be installed in the existing building of ESTIA Institute of Technology in Bidart, France. This medium-sized public building will be equipped with PV solar panels, Li-ion batteries and a hydrogen energy storage system with technical specifications presented in Table 1. The resulting microgrid is represented in Fig. 4, where the different elements and their connections are displayed, as well as the control structure, discussed further in Section 3.

Since the objective of the building microgrid is to supply as much of its energy demand as possible with renewable energy, the PV panels will be controlled via the maximum power point tracking algorithm [START_REF] Bendib | A survey of the most used MPPT methods: conventional and advanced algorithms applied for photovoltaic systems[END_REF]. The PV power generation was modelled using real profiles of solar irradiation and temperature in Bidart [START_REF]JRC Photovoltaic Geographical Information System (PVGIS) -European Commission[END_REF], resulting in the annual power generation shown in Fig. 1. On the other hand, the building energy consumptionshown in Fig. 2 follows the energy demand profile of the ESTIA building, in which weekends and holidays were considered for conceiving the dataset. It is worth noting that the building demand is about 25 % more elevated during the winter than in summer days due to heating devices.

Due to low PV energy capacity, the energy exchange between the grid and the building microgrid is moderated by a community aggregator [START_REF] Moret | Energy collectives: a community and fairness based approach to future electricity markets[END_REF]. In this configuration, prosumers send their day-ahead forecast load consumption and energy injection to the aggregator. After that, the aggregator buys (or sells) electricity from the market operator at wholesale market prices and sells to (or purchases from) prosumers in the microgrid community at retail market prices. Consequently, the building EMS must provide a one-hour resolution day-ahead grid energy exchange plan to the community aggregators.

One important restriction imposed by French grid regulations is the limitation concerning the annual energy injected into the grid [START_REF]et l'exploitation d'Installations de production d'électricité à partir d'énergies renouvelables en autoconsommation et situées en métropole continentale[END_REF]. The current grid code allows building microgrids to inject only the renewable energy generated locally and restricts building microgrids from discharging their ESSs to inject into the grid. Moreover, the building microgrid must be well-sized to minimise the annual PV energy generation injected into the public grid.

Therefore, to reduce the energy injection and maximise the annual PV self-consumption rate, the building microgrid has to optimise its power flow by assigning the right setting points for its hybrid ESS. Particularly, to handle the seasonality throughout the year, the building microgrid under study will rely on a hydrogen energy storage which can reserve energy long-term, which enables shifting the summer energy surplus to the periods of energy deficit during the winter. In particular, the building microgrid is equipped with PEM technology for electrolysis and fuel cells. This kind of technology can be coupled with RESs, thanks to its faster response time and lower degradation rate when subjected to an intermittent power rate, compared to other technologies [START_REF] Rakousky | Polymer electrolyte membrane water electrolysis: restraining degradation in the presence of fluctuating power[END_REF].

The power-to-power hydrogen ESS includes three main parts: a hydrogen producer device (fuel cell), a hydrogen consumer device (electrolyser), and the hydrogen reservoir (tank), as shown in Fig. 4. The electrolysis reaction happening in the PEM electrolyser cells converts water and electricity into heat, oxygen and hydrogen gases, as described in Eqs. ( 1)-( 3) [START_REF] Abdin | Modelling and simulation of a proton exchange membrane (PEM) electrolyser cell[END_REF]. As shown in Fig. 3a, the electrons pass through the electrical circuit attached to the electrolyser, while protons pass through the membrane dividing the anode and cathode. Subsequently, the protons arriving through the PEME membrane at the cathode combine with the electrons arriving through the ancillary electrical circuit, generating hydrogen in the gaseous form.

Complete electrolysischemical reaction :

H 2 O→H 2 + 1 2 O 2 (1)
Anode :

H 2 O→2H + + 1 2 O 2 + 2e - (2)
Cathode :

2H + + 2e -→H 2 (3)
The hydrogen produced in the cathode side of the PEM electrolyser is pushed into the hydrogen tank by using a hydrogen compressor [START_REF] Wang | Modelling and control of hybrid wind/photovoltaic/fuel cell distributed generation systems[END_REF]. The compressor adapts the difference of pressure between the PEME cathode and the hydrogen reservoir pressure, enabling the hydrogen to flow inside the tank. The stored hydrogen can be used later to generate electric power according to the building needs by regulating the operating pointing of PEM fuel cells and the outlet hydrogen flow.

On the other hand, the chemical energy of the combustion of hydrogen with oxygen happening in the PEM fuel cells, produces water, electricity and heat, as specified in Eqs. ( 4)-( 6) [START_REF] Abdin | Modelling and simulation of a proton exchange membrane (PEM) electrolyser cell[END_REF]. As illustrated in Fig. 3b, the electro-oxidation of the hydrogen happens in the anode side, while the electro-reduction of oxygen occurs in the cathode side. The regulation of the outlet hydrogen throughput is assured by setting up the position of a mechanical valve located in between the tank and the anode of PEM fuel cell, whereas the operating point of the PEMFC is assured by controlling its DC-DC power converter. One important remark to assure high PEM fuel cell electrochemical efficiency is that the PEM fuel cell hydrogen utilisation rate must be regulated to make the PEMFC consume all the hydrogen flowing throughout the cells. In this way, all hydrogen crossing the PEMFC is used to generated power, avoiding resources being wasted.

Complete combustion chemical reaction :

H 2 + 1 2 O 2 →H 2 O ( 4
)
Anode :

H 2 →2H + + 2e - (5)
Cathode :

1 2 O 2 + 2H + + 2e -→H 2 O (6)

Hierarchical model predictive control structure

To adapt ESTIA's building microgrid into the electrical grid structure building microgrid. The upper EMPC optimises the flow of energy, giving references to the lower MPC, in charge of managing the power flow. The AOHS measures the current behaviour of the hydrogen storage system and communicates dynamic parameters to both MPC algorithms.

Taking into account the stochasticity in the net power imbalance forecast (E pv -E cons ), the upper MPCnamed Economic MPC (EMPC)determines the grid energy planning (E grid, k import and E grid, k inject ) for each discrete period k, where one period is an hour long (T s = 1 h), the batteries state of charge (SoC ref, k ), and the tank level of hydrogen (LoH ref, k ). The objective of the EMPC is to minimise both the day-ahead energy imported (E grid, k import ) and the energy injected into the grid (E grid, k inject ). This is guaranteed through the optimisation of the cost function expressed in Eq. ( 7), solved using Mixed-Integer Linear Programming of the CPLEX.

min θ ∑ 48 k=1 E import grid,k + E inject grid,k where θ = { SoC ref , LoH ref ,E import grid , E inject grid } (7)
The optimisation of Eq. ( 7) is implemented every day, at midnight or whenever the error between the planned energy exchange and the real imported and injected energy is higher than 7 %. This mechanism aims to reduce the number of unnecessary optimisations of the EMPC cost function.

Notably, the EMPC sends references for the battery state-of-charge (SoC) and the level of hydrogen (LoH) of the tank to the lower MPC, whereas the first 24 h of the grid energy planning (E grid import and E grid injected ) are sent to the community aggregator. A horizon of 48 h was chosen to better estimate the SoC and the LoH at the end of the day. A prediction horizon twice longer than the daily optimisation prevents the controller from fully discharging its ESS unnecessarily, therefore keeping the optimality only for the current day but also for the next one.

Meanwhile, the lower MPC aims to fulfil the grid exchange plan sent to the community aggregator by running the hybrid ESS at adequate power set points. Therefore, every hour it implements a tracking MPC (TMPC) to minimise the quadratic errors between the LoH and SoC references transmitted by the EMPC, as expressed in the cost function defined in Eq. [START_REF] Bartels | Influence of hydrogen on grid investments for smart microgrids[END_REF]. Its optimisation horizon is eight times shorter than the one of the upper MPC, so as to reduce the need for high computation resources. In this equation, each quadratic error is normalized regarding its respective maximum values and weighted to give more importance to instantaneous references than to farther references in the TMPC horizon. Consequently, TMPC will determine the power references for fuel cells (P fc ), electrolysers (P els ) and batteries (P bat ) that follow as much as possible the references of the upper MPC considering the updated ESS states.

min Pfc ,Pels,Pbat ∑ 6 k=1 ( N TMPC h -k -1 SoC max ) 2 ( SoC ref ,k -ŜoC k ) 2 + ( N TMPC h -k -1 LoH max ) 2 ( LoH ref ,k -LoH k ) 2 (8)
One of the most important goals of the hierarchical MPC is to keep the energy balance between load (E cons ) and generation (E pv ) by considering the equality constraint [START_REF] Ren | Optimal operation of a grid-connected hybrid PV/fuel cell/battery energy system for residential applications[END_REF]. In this equation, the power consumed by the hydrogen compressor (P comp ) to push the hydrogen into the tank is considered whenever the electrolyser is operating (δ els = 1). Furthermore, the hierarchical MPC optimisation is constrained to maintain the safe operation of each ESS by limiting their power rate according to its maximum and minimum values, as expressed in Eqs. [START_REF] Lokar | The potential for integration of hydrogen for complete energy self-sufficiency in residential buildings with photovoltaic and battery storage systems[END_REF] for the hydrogen storage and (11) for batteries. Notably, δ els , δ fc , δ bat ch and δ bat dis are integer variables that are worth one when the ESS is active;

and zero, otherwise. To avoid charging and discharging batteries and prevent the controller from filling in and filling out the hydrogen tank simultaneously, the inequality constraints in Eq. ( 12) must be embedded into the two-level MPC.

E pv,k -E cons,k + T s ( P bat,k + P els,k + P fc,k + P comp • δ els,k ) + E import grid,k + E injected grid,k = 0 (9) -P max els,k • δ els,k ≤ P els,k ≤ 0; 0 ≤ P fc,k ≤ P max fc,k • δ fc,k (10)
P min bat,k • δ ch bat,k ≤ P ch bat,k ≤ 0; 0 ≤ P dis bat,k ≤ P max bat,k • δ dis bat,k (11)
0 ≤ δ els,k + δ fc,k ≤ 1; 0 ≤ δ ch bat,k + δ dis bat,k ≤ 1 (12)
The hierarchical MPC also considers the grid code for small prosumers in France [START_REF]Enedis l'électricité en réseau, Conditions de raccordement des Installations de stockage[END_REF] by limiting the charging of ESSs with the internal power surplus (P surplus ) and constraining their discharging to only supply the building internal power deficit (P deficit ), as defined in Eqs. ( 13) and ( 14), respectively.

-

⃒
) ⃒ ⃒ ⃒ (14)
The hierarchical MPC structure decides to use the batteries or the hydrogen storage according to their round-trip energy efficiency, which are embedded into the MPC internal models through the inequality constraints defined in Eqs. ( 15) and [START_REF] Hu | Model predictive control of microgridsan overview[END_REF]. In these two equations, the values of SoC and LoH for the next period (SoC k+1 and LoH k+1 ) are computed via the linear functions f bat and f H2 that are calculated through autonomous observers.

In the case of batteries, an iterative algorithm that is further detailed in [START_REF] Yamashita | Real-time Parameters Identification of Lithium-ion Batteries Model to Improve the Hierarchical Model Predictive Control of Building MicroGrids[END_REF] determines the parameters of f bat . On the other hand, the parameters of f H2 are updated daily by the AOHS algorithm detailed in the next section. The AOHS was briefly introduced in the authors' previous work [START_REF] Yamashita | Two-level hierarchical model predictive control with an optimised cost function for energy management in building microgrids[END_REF] to explain the context of another study focusing on a technical-economic analysis of installation of hybrid ESS in building microgrids. In this paper, the AOHS is mathematically detailed, and its operation performance is analysed in depth. Additionally, contrary to the work presented in [START_REF] Yamashita | Two-level hierarchical model predictive control with an optimised cost function for energy management in building microgrids[END_REF], in this paper, the AOHS robustness against measurement noise is assessed to verify its impact on the annual selfconsumption rate.

SoC min ≤ SoC k+1 = SoC k + f bat ( P ch bat,k , P dis bat,k ) ≤ SoC max ( 15
)
LoH min ≤ LoH k+1 = LoH k + f H2 ( P els,k , P fc,k ) ≤ LoH max (16)

Algorithm of the AOHS

As shown in Fig. 4, in a closed loop with the hierarchical MPC structure and the hydrogen storage, the AOHS aims to continuously improve the precision of the LoH estimation. The proposed AOHS updates the parameters θ els and θ fc composing f H2 defined in Eq. ( 17). This operation is performed with the same sample time as the smallest controller's discretisation time (T s = 1 h) and relies solely on local measurements. The θ els and θ fc parameters are determined either through physical values provided by the manufacturers (ς els ,ς fc and V tank defined in Table 1) or through the process implemented by the AOHS.

f H2 ( P els,k , P fc,k ) = T s ς els ⋅V tank ⋅1000 ⏟̅̅̅̅̅̅̅̅̅ ̅⏞⏞̅̅̅̅̅̅̅̅̅ ̅⏟ θels,k P els,k - ς fc ⋅T s V tank ⋅1000 ⏟̅̅̅̅̅̅⏞⏞̅̅̅̅̅̅⏟ θfc,k P fc,k (17)
The continuous identification of θ els and θ fc is based on two iterative correlation processes and on a step of model polishing. Each of these operations will be detailed in the next subsections.

Step 1: correlation between current and pressure variation

The relationship between the current and the tank pressure is inferred from the combination of Faraday's law of electrolysis with the ideal gas law. Defined by Eq. ( 18), Faraday's law of electrolysis links the current in the PEM stacks (i PEM in A) to the hydrogen flowing ( ṅPEM H2 in mol/s) [START_REF] Njoya | A generic fuel cell model for the simulation of fuel cell vehicles[END_REF]. In this equation, N cells PEM is the number of cells in the PEME or the PEMFC, and F is the Faraday constant in C/mol. Notably, the notation PEM refers to the set PEM = {fc, els}, referring to PEMFC or PEME, respectively.

ṅPEM H2 = N PEM cell ⋅i PEM 2⋅F (18)
Considering that the temperature variation in the hydrogen tank is much slower than the model update time (T s ), the number of moles stored in the tank (n H2 tank ) is proportional to the tank pressure (P tank ), following the ideal gas law defined by Eq. ( 19), where R is the ideal gas constant in m 3 •Pa K•mol . By computing the derivative, the hydrogen flow ( ṅPEM H2,k ) can be estimated from the variation of tank pressure. As a result, Faraday's equation in Eq. ( 18) can be rewritten in function of the tank pressure, instead of the number of moles, as shown in Eq. [START_REF] Abdin | PEM fuel cell model and simulation in Matlab-Simulink based on physical parameters[END_REF].

n tank H2,k P tank,k = V tank R⋅T tank → d dt ṅtank H2 = V tank R⋅T tank ⋅ P tank,k -P tank,k-1 T s ( 19
)
ΔP tank,k = R⋅T tank V tank ⋅ N PEM cells 2⋅F ⋅T s ⏟̅̅̅̅̅̅̅̅̅̅̅⏞⏞̅̅̅̅̅̅̅̅̅̅̅⏟ β PEM Faraday ⋅i PEM ( 20
)
Instead of attributing static values to compute β Faraday PEM in Eq. ( 20), the AOHS algorithm determines it automatically thanks to a dynamic update process. To identify β Faraday, k PEM , the pressure and current measurements are used. The sample measurement composed by the couple M k = (ΔP tank, k meas , i PEM, k meas ) are used to calculate β Faraday, k meas, PEM , by using Eq. ( 21), where PEM = {els, fc} refers to either electrolyser or fuel cells.

β meas,PEM Faraday,k = ΔP meas tank,k i meas PEM,k , if i meas PEM,k ∕ = 0 (21)
With the instantaneous value of β Faraday, k meas, PEM

, the estimated value of β Faraday, k PEM is updated through a weighted average between its last estimation and the current measurement, as specified in Eq. ( 22). As shown in Fig. 5 and defined in Eqs. ( 23) and ( 24), the variable ω PEM, k is a confidence weight that grows with the acquisition of new measurements and is linearly reduced by a constant variable τ PEM decay . The value of τ PEM decay should be tuned according to the desired model robustness against noise measurement. Higher values of τ PEM decay will lead the model to adapt to changes in the external environment more easily, but it will reduce the model's robustness. Notably, the static parameter Δω PEM -also known as the learning rateis the variation of the confidence weight between two consecutive iterations. Therefore, the weight ω PEM, k increases as soon as a new measurement M k is acquired, following Eq. [START_REF]JRC Photovoltaic Geographical Information System (PVGIS) -European Commission[END_REF]. Subsequently, at each iteration (T s = 1 h), all confidence weights are reduced according to a pace defined by Eq. [START_REF] Moret | Energy collectives: a community and fairness based approach to future electricity markets[END_REF]. Fig. 5 shows that the weights ω PEM={els,fc} decrease when the PEMFC or PEME are not operating.

β PEM Faraday,k+1 = ω PEM,k ⋅β PEM Faraday,k + Δω PEM ⋅β meas,PEM Faraday,k ω PEM,k + Δω PEM (22) ω PEM,k = ω PEM,k + Δω PEM , if i PEM,k ∕ = 0 (23) ω PEM,k = ω PEM,k -τ decay PEM (24)

Step 2: correlation between current and power

The relation between the current and the power is non-linear. Consequently, representing the hydrogen flow through the tank's inlet and outlet as a single linear function, as suggested in most studies [START_REF] Valverde | Integration of fuel cell technologies in renewableenergy-based microgrids optimizing operational costs and durability[END_REF][START_REF] Garcia-Torres | Optimal economical schedule of hydrogen-based microgrids with hybrid storage using model predictive control[END_REF], may result in modelling inaccuracies. Considering this problem, a new intermediary linear model based on the intensity of power is proposed. This novel linear model divides the power references assigned to electrolysers (P els ) and fuel cells (P fc ) into three operating zones, which leads to the definition of three powers for each PEM technology, namely P PEM zone1 , P PEM zone2 and P PEM zone3

. The objective is to enhance the model accuracy by linearising the non-linear power-current curve shown in Fig. 6 into three linear functions.

Without storing any past measurements, the identification of these three lines relies only on the present measurements of the currents (i fc and i els ) and of the voltages (v fc and v els ) of PEMFC and PEME, and on the coordinates of three points, namely A, B and C. These three points are the vertices of the intersection between the power zones and the linear model, as shown in Fig. 6a. Therefore, as summarised in Table 2, their ycoordinates are fixed and dependent on the maximum power rate supported by each PEM technology (P PEM max ), whereas their x-coordinates are defined by angular coefficients of each of the three lines, i.e. α PEM zone1 , α PEM zone2 and α PEM zone3

. Consequently, the linear model is determined by the AOHS through the identification of these three angular coefficients.

At the first iteration, the vertices A, B, and C are initialised so that the three vertices are aligned with the first non-null measure M k , as shown in Fig. 6a. Therefore, if the first point M k = (i PEM, k meas , v PEM, k meas ), in which i PEM, k meas ∕ = 0, then the three angular coefficient are worth α els, k zone1 = α els, k zone2 = α els, k zone3

= 1/v PEM, k meas .
Subsequently, for the next iterations, the vertices are updated according to the location of the measurement M k . If the point M k belongs to a specific power zone, then one of the vertices limiting this power zone will be modified. For instance, if the M k is in zone3, as shown in Fig. 6b or c, either the vertices B or C will be modified. The upper vertex (point C) will be updated, if M k is nearer to the vertex C than vertex B, as shown in Fig. 6b. On the opposite case, the lower vertex (point B) will be modified if the point M k is nearer to vertex B than vertex C, as illustrated in Fig. 6c.

To update the upper vertex C, the angular coefficient of the third zone, i.e. α PEM zone3

, is calculated. Implemented through Eq. ( 25), this updating process corresponds to a weighted dynamic averaging operation of the last estimation of the angular coefficient of the line BC (i.e. 25). To deal with the ageing of the PEMFC and PEME, at each sampling time (T s = 1 h), the confident weights linearly decrease according to the variable τ decay , as defined in Eq. [START_REF] Wang | Modelling and control of hybrid wind/photovoltaic/fuel cell distributed generation systems[END_REF].

α zone3 PEM,meas,k = C PEM y,k -v meas PEM,k C PEM x,k -i meas PEM,k (26) ω zone3 PEM,k = ω zone3 PEM,k + Δω PEM (27) ω zone3 PEM,k = ω zone3 PEM,k -τ decay ( 28
)
In the opposite scenario, if the point M k belongs to zone3 and is nearer to vertex B than vertex C, the x-coordinate of the vertex B is modified. Firstly, as shown in Fig. 6c, an intermediate point

B ′ = (B x ′ , B y
′ ) is determined. This point B ′ corresponds to the intersection between

CM k and the line dividing zone 2 and zone 3. Subsequently, the point B is updated similarly to the upper vertex. Firstly, the angular coefficient of line AB ′ , namely α fc, meas zone2 , is calculated through Eq. ( 29). Thereafter, the same update process of the upper vertex is applied to the lower vertex.

α zone2 PEM,meas,k = B ′ y -A fc y,k B ′ x -A fc x,k (29)
To better understand the second step of the algorithm, Fig. 7 shows the temporal evolution of the confident weights. As expected, the weights ω fc zone1 , ω fc zone1 and ω fc zone1 grow when the PEME or PEMFC are operating and decrease when they are in an idle mode. After having identified the three angular coefficients (α PEM zone1 , α PEM zone2 and α PEM zone3

), the final step consists of determining a unique correlation parameter, namely σ PEM , linking power to current, as expressed in Eq. [START_REF] Yamashita | Real-time Parameters Identification of Lithium-ion Batteries Model to Improve the Hierarchical Model Predictive Control of Building MicroGrids[END_REF]. Instead of building a model with multiple variablesfor instance, power at zone1, zone2 and zone3 -the idea is to facilitate the MILP optimisation process by including only two control variables, i.e. power references for PEMFC (P fc ) and PEME (P els ).

P PEM,k = σ PEM,k ⋅i PEM,k (30)
Therefore, from the values of the three angular coefficients, it is possible to determine σ PEM, k through the minimisation of the weighted

Table 2

Coordinates of the vertices points A, B, and C. , as illustrated in Fig. 8) and the approximative linear model defined by Eq. ( 30). This approximation is assured by the minimisation of the cost function defined by Eq. [START_REF] Yamashita | Two-level hierarchical model predictive control with an optimised cost function for energy management in building microgrids[END_REF].

min σPEM ∑ 3 i=1 ω zonei PEM ⋅ ( G x zonei - (2⋅i -1)⋅P max PEM 6⋅σ PEM ) 2 (31)

Step 3: approximation and definition of the linear model

The final step of the AOHS algorithm consists of building the final model by identifying the direct correlation between power (P els and P fc ) and variation of the LoH (ΔLoH = f H2 (P els , P fc )). The real LoH can be estimated from the current tank pressure (P tank, k meas ) and the maximum tank pressure (P tank max ), as shown in Eq. [START_REF] Njoya | A generic fuel cell model for the simulation of fuel cell vehicles[END_REF]. By combining this equation to the correlation between current and pressure variation determined in the first step (β Faraday PEM of Eq. ( 22)) and the correlation between current and power (σ PEM of Eq. ( 31)) calculated in the second step, it is possible to determine the direct link between power and variation of the LoH. Consequently, the function f H2 in Eq. ( 17) can be calculated using Eq. (33).

LoH

Case studies

To assess the performance of the MPC empowered with the AOHS algorithm, the building microgrid of sizing defined in Table 1 was simulated in MATLAB Simulink® using the PEMFC SimPower-Systems model and a PEME model validated through experimental data reported in [START_REF] García-Valverde | Simple PEM water electrolyser model and experimental validation[END_REF]. The analysis of the proposed EMS is divided into two main parts. The first part is focused on the robustness of the proposed linear model against measurement noise. Meanwhile, the second part consists of a benchmark with the common linear model [START_REF] Garcia-Torres | Optimal economical schedule of hydrogen-based microgrids with hybrid storage using model predictive control[END_REF] defined by Eq. ( 17).

Precision and robustness against measurement noise

To verify the precision and the robustness of the proposed AOHS against measurement noise, the real voltage, current, and tank pressure measurements were disturbed by a normalized Gaussian noise of variance of 1 %, 5 %, and 10 % of their nominal values. To assess the error between the day-ahead LoH calculated by the EMPC and the real one, the TMPC in these scenarios was considered as a perfect router. Consequently, instead of optimising Eq. ( 8), it implements the control variables determined by the EMPC. Therefore, the effect of the TMPC is decoupled, highlighting the impact of errors in the estimated LoH.

After simulations of 50 days of operation of the PEMFC and PEME under these three measurement noise intensities, the temporal evolution of the current-voltage (σ fc and σ els ) and current-pressure (β Faraday els and β Faraday fc ) correlation factors are shown in Fig. 10. It is noteworthy that in these simulations, Δω PEM and τ decay were set up as 0.1 and 0.001, respectively. In addition,ω PEM , ω PEM zone1 , ω PEM zone2 and ω PEM zone3 were initialised as 0.5. Table 3 indicates that the AOHS can identify the current-power coefficients with an average error below 2 % regarding the theoretical model, whereas the pressure-current is determined with an average error inferior to 11 %. Notably, the theoretical model was obtained from the complete linear regression of the real PEME and PEMFC models. Despite including a measurement noise with a variance of 10 %, the graphs in Fig. 10 highlight that the final estimated hydrogen storage parameters are close to the theoretical model, demonstrating that the proposed AOHS algorithm is robust, especially against voltage and current noise measurements. Nonetheless, since the algorithm is based on the variation of the pressure (Eq. ( 21)), the pressure measurement noise is amplified, harming the precision of pressure-current parameter identification (β Faraday els and β Faraday fc ), as shown in Fig. 9. Besides its robustness, the precision of the proposed algorithm is comparable to the theoretical model. According to Table 3, the algorithm can assure almost the same average error under the three intensities of noise measurement, with an average difference below a 0.6 percentage point (p.p.) regarding the theoretical model. In particular, the AOHS demonstrated highly robust against noise measurement with variance inferior to 1 %, since the maximum and the average LoH errors were kept almost the same with 0 % and 1 % of noise variance.

Notwithstanding, as shown in Table 3, the maximum error in the estimation of the LoH when including noise measurements with a variance of 5 % and 10 % is higher than the theoretical model. As discussed previously and shown in Fig. 10, this is due to the inaccuracies in the estimation of the pressure-current correlation factors created by the computation of the pressure variation. Therefore, to keep the same precision as the theoretical model, a pressure measurement noise below 5 % must be guaranteed. This can be assured through the inclusion of moving average filters or using more precise pressure sensors.

Comparison with the time-invariant model

The proposed MPC with the AOHS is compared to an MPC with the time-invariant model specified in Eq. ( 17). This comparison aims to verify the impact of parametric errors on the annual building microgrid self-consumption and coverage rates, as well as the accuracy of dayahead grid energy planning transmitted to community aggregators.

In this view, the building microgrid and its HMPC (EMPC and TMPC) were simulated in MATLAB Simulink® for 365 days under four scenarios. The first three scenarios are when the HMPC is equipped with a time-invariant hydrogen storage model holding +20 %, 0 %, and -20 % parametric errors in θ fc and θ els defined in Eq. [START_REF] Kamal | Model predictive control and optimization of networked microgrids[END_REF]. Finally, the fourth scenario is when the HMPC is endowed with the AOHS algorithm with a noise measurement of variance of 5 %. Since no parameter needs to be tuned when employing the AOHS algorithm, it is not necessary to analyse the impact of parametric errors.

τ sc = 1 - ∑ 8760 k=1 E injected grid,k ∑ 8760 k=1 E pv,k ; τ c = 1 - ∑ 8760 k=1 E import grid,k ∑ 8760 k=1 E cons,k (34)
By comparing the annual self-consumption and annual coverage rates defined by Eq. (34) summarised in Table 4, it is possible to observe that the HMPC with the AOHS algorithm assures practically the same coverage and self-consumption rates as the time-invariant model with no parametric errors. Moreover, by comparing the annual coverage rate of these four simulation scenarios, it is possible to note that it is almost insensitive to parametric errors. In contrast, the self-consumption rate is reduced by 3 p.p. with a parametric error of -20 %, and it was unchanged with a parametric error of +20 %.

The reason behind this result is that the power reference of the PEME is almost the same (around 30 kWh) in the scenarios with the AOHS algorithm and the time-invariant models with +20 % and 0 %, as shown in Fig. 11. According to Fig. 11, the EMPC minimises the grid energy exchange by using its ESS as much as possible. Consequently, to maximise the self-consumption rate τ sc , the EMPC absorbs the PV energy surplus by storing hydrogen through the operation of the PEME.

As shown in Fig. 11, as in all scenarios the hydrogen tank is not full (the LoH is below 90 %), the EMPC can operate the PEME most of the time at maximum power rate (around 30 kWh) to minimise the grid energy injection. Since there are low parametric errors when employing the AOHS and the time-invariant models with 0 % (Fig. 11b andd), the TMPC can follow the LoH references strictly, resulting in an operation of the PEME close to the maximum power rate, as expected from the EMPC. On the contrary, even though the LoH reference is not followed by the TMPC when θ els is over estimated (Fig. 11c), the power set point of the PEME is equal to the EMPC estimated values. This is because the overestimation of θ els leads the LoH references to be always higher than the real capacities of the PEME. Therefore, this will induce the PEME to be operated over its estimated capacity, which is bounded by its maximum power rate. Consequently, the power references of both the TMPC and the EMPC are equivalent, assuring similar self-consumption rates.

On the other hand, when θ els is underestimated (Fig. 11a), the errors in the LoH estimation calculated by the EMPC will lead the PEME to be underused. In this scenario, the EMPC tries to operate the PEME at maximum power rate (around 30 kWh), but due to parametric errors in θ els it will transmit to the TMPC the references of the LoH that are smaller than the real hydrogen production rate. Therefore, the PEME will be underexploited, resulting in the reduction of the self-consumption rate.

It is important to note that a similar phenomenon happens with the operation of fuel cells. Nonetheless, due to the low round-trip efficiency (below 40 %) of hydrogen conversion, the use of fuel cells is minimal. According to graphs shown in Fig. 12, even though the annual use of fuel cells can vary up to 17 % (1.3 MWh) regarding the time-invariant model with no errors, this difference is small compared to the annual building consumption, which amounts to 242 MWh/year. Consequently, the coverage rate (τ c ) defined by Eq. ( 34) is almost unaffected by the inaccuracies in the PEMFC models.

Although the imprecisions in the LoH estimation do not always impact the annual self-consumption and coverage rates, it affects the grid energy planning sent to community aggregators. As shown in Fig. 13, compared to time-invariant models with no parametric errors, the error in the day-ahead grid energy exchange planning sent to the community aggregator increases by respectively 55 % and 10 % with time-invariant (T.I.) models holding -20 % and +20 % of parametric imprecision. Conversely, the proposed AOHS algorithm can assure an average error of 334 Wh, which is 52 % lower than the ideal timeinvariant. This result indicates that the AOHS algorithm can facilitate grid energy management from the point of view of the external grid by providing more reliable day-ahead grid energy planning to community aggregators.

Conclusion

In this paper, a two-level hierarchical model predictive controller (HMPC) with an Autonomous Observer of Hydrogen Storage was proposed to reduce the HMPC internal model inaccuracy. This modelling inaccuracy issue is particularly important when considering long operating times. The proposed solution incorporates an external observer to the HMPC structure, requiring little to no modification of the MPC algorithms while offering better performance. Even subjected to strong measurement noise, the designed energy management proved capable of identifying the real hydrogen production and consumption rates. This operation was performed while relying only on a few storage variables and requiring no tuning of the parameters. This facilitates its implementation in embedded systems. It could guarantee an average precision inferior to 2 % points regarding the theoretical model that holds the full knowledge of the system.

Table 4

Comparison of the coverage and self-consumption rates when using timeinvariant models and the AOHS algorithm. The flexibility of the algorithm in adapting automatically to external environments reduces the negative effects of model parameter inaccuracies provided by manufacturers. Through simulations, it was noticed that a 20 % underestimation of hydrogen production reduces the building's annual self-consumption by about 3 % when using a rudimentary model approach. On the other hand, through a more accurate estimation of the tank level hydrogen provided by the novel algorithm, the self-consumption and coverage rates are similar to the ideal model that holds the full system knowledge. Additionally, the algorithm enables the HMPC to provide reliable day-ahead grid energy planning to community aggregators, facilitating the external grid energy management and reducing the errors by about 52 % compared with the conventional modelling approach.

Since the ESTIA's building microgrid studied in this paper is under construction, the future work consists of testing the algorithms and solutions described in this paper using hardware in the loop approach. The fact that the three elements of the control structure described in this paper, EMPC, TMPC and AOHS, are operating with different time scales will make possible to implement them using different machines or different cores of the same real-time computer for each control element.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Photovoltaic power generation during an entire year, from January to December.

Fig. 3 .

 3 Fig. 3. Summary of chemical reaction in cathode and anode sides of (a) proton exchange membrane electrolysis cell and (b) proton exchange membrane fuel cells.

Fig. 4 .

 4 Fig. 4. Energy management system to optimise the power flow of a building microgrid.

  new measured α fc, meas zone3 . The variable α fc, meas zone3 is the angular coefficient of the line BM k that is calculated through Eq. (26), whereas ω fc, k zone3 is the confidence weight of the estimation of α fc zone3 . The term Δω fc also called learning rateindicates the importance of the new measurement regarding the current estimation of α fc zone3 . Both ω fc zone3 and Δω fc are initialised to a small value, typically in the order of 0.1. The weight ω fc zone3 grows with the acquisition of measurements close to the vertex C, following Eq. (

Fig. 5 .

 5 Fig. 5. Temporal evolution of pressure-current confident weights.

  ⋅α zone3 PEM,k + Δω PEM ⋅α zone3 PEM,meas,k ω zone3 PEM,k + Δω PEM[START_REF]et l'exploitation d'Installations de production d'électricité à partir d'énergies renouvelables en autoconsommation et situées en métropole continentale[END_REF]

Fig. 6 .

 6 Fig. 6. Detail of the fuel cell RTMI. (a) Initialisation step. (b) Update step when the measurement point is near to the upper vertex. (c) Update step when the measurement point is near to the lower vertex.

Fig. 7 .

 7 Fig. 7. Temporal evolution of power-current confident weights.

Fig. 8 .

 8 Fig. 8. Visualisation of the definition of the approximate model.

Fig. 9 .

 9 Fig. 9. Estimation of current-power (σ fc and σ els ) and current-pressure (β Faraday els

Fig. 10 .

 10 Fig. 10. Comparison between the final estimated model correlation factors (σ PEM and β Faraday PEM ) and their corresponding theoretical values in the case with noise measurement variance of 10 %.

Fig. 11 .

 11 Fig. 11. Comparison between the EMPC and TMPC control variables with the AOHS algorithm subjected to measurement noise of 5 % and with time-invariant models subjected to parametric errors of -20 %, 0 %, and +20 %.

Fig. 12 .

 12 Fig. 12. Annual use of energy storage systems when using time-invariant models with parametric errors and the AOHS algorithm with noise measurement of 5 %.

Fig. 13 .

 13 Fig. 13. Average (red squares) and standard deviation (blue bars) of the error in the day-ahead grid energy exchange planning sent to community aggregator. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1

 1 Investigated building microgrid electrical components.

	Equipment	Technical specification
	Photovoltaic panels	Peak power at 1000 W/m 2 : 107 kWc
		Annual energy generation: 131 MWh
	Building load	Annual energy consumption: 242 MWh
	Li-ion batteries	Nominal capacity: 167 Ah
		Nominal voltage: 720 V
		Nominal discharge current: 70A
		Maximum power rate: 60 kW
	Hydrogen tank	Maximum pressure: 30 bars @ 80

• C Maximum hydrogen mass: 18 kg Normalized volume (V tank ): 224 Nm 3 PEM electrolyser Nominal power: 25 kW Hydrogen production (ς els ): 4.18 kWh/Nm 3 Maximum power rate: 30 kW PEM fuel cells Nominal power: 20 kW Hydrogen consumption (ς fc ): 0.63 Nm 3 /kWh Maximum power rate: 48 kW Compressor Nominal power consumption (P comp ): 1 kW

Table 3

 3 Average precision of the level of hydrogen estimation according to the intensity of noise measurement.

	LoH error	Noise measurement intensity
		Theoretical	0 %	1 %	5 %	10 %
	Average (p.p.)	0.57	0.77	0.75	1.2	0.68
	Maximum (p.p.)	4.3	4.9	4.9	6.8	7.2

Acknowledgement

The authors would like to thank the New Aquitaine Region (AAPFO424953), ESTIA Institute of Advanced Industrial Technologies and the University of Poitiers for their financial support.

CRediT authorship contribution statement

Paper revision and discussion over the review comments.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Daniela Yassuda Yamashita reports financial support was provided by Region Nouvelle-Aquitaine. Jean-Paul Gaubert reports a relationship with University of Poitiers that includes: employment. Ionel Vechiu reports a relationship with ESTIA-RECHERCHE that includes: employment. Daniela Yassda Yamashita reports a relationship with ESTIA-RECHERCHE that includes: employment.


Joel Serra

Prajwal Christophe Bouvet
email: christophe.bouvet@isae-supaero.fr

Karinja Haridas

Léon Ratsifandrihana

P Karinja Haridas

Numerical simulations of combined size effects acting on an open-hole laminated composite plate under tension

Keywords: Size effect, Open-hole tensile test, Thickness effect, Blocked plies, Discrete Ply Modeling, Finite Element Model Size effect, Open-hole tensile test, Thickness effect, Blocked plies, Discrete Ply Modeling

Numerical and experimental research programs have been carried out to investigate the effect of scaling on the tensile strength of notched composites. This paper presents a computational study of scaled openhole tensile tests using the Discrete Ply Modeling (DPM) method. This finite element model is discrete, and only a small number of parameters are required from experimental characterization tests. Experimental and numerical strength values are compared here, and reveal that DPM simulations tend to slightly overestimate strength values, with an average discrepancy of 9.7%. However, DPM Results show that such modeling simulates both the reduction in strength when specimen size is increased for sublaminate level scaled specimens, where failure is fiber dominated, and the increase in strength when specimen size is increased for ply level scaled specimens, where failure is delamination dominated. In all cases, increasing the total thickness of the specimen leads to a decrease in strength and this effect is dominant over the effect of increasing hole diameter. As well as the variation in strength, three distinct failure mechanisms are observed: fiber failure with extensive matrix damage (pull-out failure), fiber failure with little or no matrix damage (brittle failure) and delamination failure. Comparisons with experiments demonstrate that tensile strengths, damage propagation scenarios and failure patterns are predicted with acceptable accuracy.

Aircraft that use composite structures have many parts containing holes or different cut-outs that induce stress concentrations and reduce the strength of the part. When they are subjected to tensile loading, damage and failure mechanisms around these holes/cut-outs are magnified because of stress concentrations. Notched strength is therefore one of the most important design drivers for composite structures [START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF].The mechanical performance of open-hole composite laminates depends on many factors, such as the size and thickness of the laminate; the hole size and its location; the plies orientations and the machining quality. Before real composite structures are tested, the effects of holes are studied on coupons 10-100 times smaller. Understanding the associated scaling effects is therefore of the utmost importance. Since the discovery of these scaling phenomena by Leonardo da Vinci in the early 1500s [START_REF] Bazant | Scaling of Structural Strength[END_REF], substantial amounts of research have been devoted to them but deep understanding is still lacking.

To explore the physical causes of scaling effect, numerous experiments have been conducted [START_REF] Bazant | Size Effect and Fracture Characteristics of Composite Laminates[END_REF]- [START_REF] Xu | An experimental investigation into size effects in quasi-isotropic carbon/epoxy laminates with sharp and blunt notches[END_REF]. Size effect can occur at different levels [START_REF] Camanho | Prediction of size effects in notched laminates using continuum damage mechanics[END_REF]. At the structural level, when thin plies are used, it has been proven that the strength of notched composite laminates decreases when the notch size increases. At the material level, the "in situ" effect has been detected and is influenced by the ply thickness. Wisnom et al. [START_REF] Wisnom | Scaling Effects in Notched Composites[END_REF] have analyzed this phenomenon on quasi-isotropic specimens .

The "hole size effect" has been studied at the structural level. In this configuration, increasing hole diameter leads to a decrease in the strength of a laminate. The presence of non-critical ply-level damage, such as matrix failure, delamination, fiber splitting and fiber fracture in the vicinity of the hole blunts the stress concentration [START_REF] Hallett | An experimental and numerical investigation into the damage mechanisms in notched composites[END_REF]. The strength difference between small and large specimens [START_REF] Lin | Fracture of filamentary composite materials[END_REF] is explained by the extent of the "fracture process zone", which is relatively independent of the size of the specimen.

Using a sublaminate level scaling method, Harris and Morris [START_REF] Harris | Role of Delamination and Damage Development on the Strength of Thick Notched Laminates[END_REF] investigated the effect of an increase in the total thickness of a laminate with a constant ply thickness. They observed that increasing thickness in quasi-isotropic laminates led to a decrease in strength because the damage prior to failure was confined to a boundary close to the surface, so thicker laminate experienced less stress redistribution throughout the laminate .

Using a ply-level scaling method, Vaidya et al. [START_REF] Vaidya | Effect of Ply Thickness on Fracture of Notched Composite Laminates[END_REF] increased the thickness of their laminates.

They observed similar fracture toughness values and similar amounts of sub-critical damage during loading for quasi-isotropic [ / / / ] s laminates with n = 1 and n = 2. However, with 0 90 45 45 n = 4, the failure mechanism changed from fiber failure to delamination failure along the 0/90 interface, resulting in the increase of sub-critical damage, accompanied by a marked decrease in fracture toughness. The influence of stacking sequence and hole diameter on the notched strength of a laminate was also investigated by Lagace [START_REF] Lagace | Notch Sensitivity and Stacking Sequence of Laminated Composites[END_REF]. He used a constant specimen width that gave varying stress concentration across the width. He found that increasing hole diameter led to a change in fracture mechanism, from a matrix-dominated to a fiber-dominated failure, for [ ] s laminates.This was 0/90 2 caused by the fact that increasing the hole-radius-to-laminate-thickness ratio led to a decrease in interlaminar stresses in the region around the hole boundary. Laminates with little or no delamination showed linear stress-strain curves up to failure whereas those with large amounts of delamination prior to failure exhibited large non-linearities. Chang et al. [START_REF] Chang | Damage Tolerance of Laminated Composites Containing an Open Hole and Subjected to Tensile Loadings[END_REF] found that damage propagates from a circular hole via matrix failure in the off-axis plies.

Preliminary sizing solutions for notched composites are mainly based on semi-empirical analytical models for predicting notched strength, such as the point stress or the average stress models [START_REF] Whitney | Stress Fracture Criteria for Laminated Composites Containing Stress Concentrations[END_REF], [START_REF] Wang | Experimental and numerical investigation of the tension and compression strength of un-notched and notched quasi-isotropic laminates[END_REF]. An extensive review of these technics has been provided by Awerbuch and Madhukar [START_REF] Awerbuch | Notched Strength of Composite Laminates: Predictions and Experiments-A Review[END_REF].Whitney and Nuismer [START_REF] Whitney | Stress Fracture Criteria for Laminated Composites Containing Stress Concentrations[END_REF] developed two related criteria based on normal stress distribution to predict the residual strength of laminated composites containing a through-thethickness discontinuity. They presented specific examples for a circular hole in quasi-isotropic material and for a notch in [ ] s glass/epoxy laminated plates. In both cases, comparison of 0/ ± 45 2 experimental data and theory showed good correlation. Various researchers have developed stressbased models, with varying levels of complexity and of success [START_REF] Lin | Fracture of filamentary composite materials[END_REF], [START_REF] Mar | Fracture mechanics correlation for tensile failure of filamentary composites with holes[END_REF]- [START_REF] Camanho | A finite fracture mechanics model for the prediction of the open-hole strength of composite laminates[END_REF]. Fast semi-analytical and numerical methods [START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavities -size effects[END_REF]- [START_REF] Barenblatt | The Mathematical Theory of Equilibrium Cracks in Brittle Fracture[END_REF] have given good agreement between experimental data and numerical predictions for quasi isotropic laminates. However, they do not take the influence of the stacking sequence into account since they consider that strain is constant within the laminate thickness [START_REF] Camanho | Numerical Modelling of Failure in Advanced Composite Materials[END_REF]. A more complete numerical model is therefore needed to predict the physics associated with failure phenomena (delamination and transverse cracks).

Many models have been proposed to fit experimental data [START_REF] Wisnom | Scaling Effects in Notched Composites[END_REF], [START_REF] Chen | Numerical analysis of size effects on open-hole tensile composite laminates[END_REF]- [START_REF] Wisnom | Modelling discrete failures in composites with interface elements[END_REF]. They can be classified according to their level of discretization [START_REF] Hongkarnjanakul | Modélisation numérique pour la tolérance aux dommages[END_REF], from the whole continuum model without any interfaces to the most discretized type of model, where the three most important types of damage (fiber failure, matrix cracking and delamination) are all represented through interfaces. To correctly model matrix cracks and their interaction with delamination, discrete elements seem necessary [START_REF] Wisnom | Modelling discrete failures in composites with interface elements[END_REF].

The best way to represent strong discontinuities, i.e. matrix cracking, should be the discrete crack methods. Numerical studies on open-hole laminate composites, conducted by Hallett et al. [START_REF] Wisnom | Scaling Effects in Notched Composites[END_REF] and Jiang et al. [START_REF] Jiang | A concise interface constitutive law for analysis of delamination and splitting in composite materials and its application to scaled notched tensile specimens[END_REF] using a simplified mesh coupled with a strong simplification of the hole geometry, confirm these assumptions. Swindeman et al. [START_REF] Swindeman | Strength Prediction in Open Hole Composite Laminates by Using Discrete Damage Modeling[END_REF] established a regularized X-FEM (Extended Finite Element Model) with a predefined minimum crack spacing. However, as the maximum number of representable cracks per ply is in strong correlation with the number of degrees of freedom, its use is limited. The model did not include fiber failure, but was able to capture matrix cracks and two load drops due to excessive delamination, captured by the cohesive zone method. Chen et al. [START_REF] Chen | Modelling the tensile failure of composites with the floating node method[END_REF] applied the floating node method and enriched 3D cohesive element to improve the connectivity between the plies and interlayer. Bao and Liu [START_REF] Bao | Progressive failure analysis on scaled open-hole tensile composite laminates[END_REF] used surface-based cohesive contact for major split cracks tangential to the hole as well as delamination. The distributed matrix cracks appeared to be captured by continuum damage mechanics. Higuchi et al. [START_REF] Higuchi | Numerical simulation of progressive damage and failure in composite laminates using XFEM/CZM coupled approach[END_REF] presented a coupled X-FEM/CZM (Cohesive Zone Model) approach where fiber failure was not modeled. A predefined crack spacing was used to model all potential matrix cracks, and delamination by interface elements. Joosten et al. [START_REF] Joosten | A cohesive network approach for modelling fibre and matrix damage in composite laminates[END_REF] developed a model in which cohesive elements are used for all failure modes, i.e. fiber, matrix and interface failure. Lu et al. [START_REF] Lu | Adaptive Discrete-Smeared Crack (A-DiSC) Model for Multi-Scale Progressive Damage in Composites[END_REF] used an adaptive discrete-smeared crack model, where the critical discrete matrix cracks were captured by the floating node method with a crack alignement algorithm.

To capture the interaction between splitting and delamination, delamination was modeled with separable cohesive elements. The model was able to predict major and minor distributed cracks and delamination migration.

Numerous models have been created to simulate open-hole laminate composite behavior but few have been used to model size effects. The Discrete Ply Model (DPM) used in this paper employs cohesive interfaces to model matrix failure and delamination, and 3D volume elements to predict fiber failure. It was initiated by Bouvet et al. for modeling low velocity impacts in composite panels [START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[END_REF], and afterwards enhanced to capture permanent indentation [START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF] and to simulate compression after impact [START_REF] Rivallant | Failure analysis of CFRP laminates subjected to compression after impact: FE simulation using discrete interface elements[END_REF]. The use of interfaces connected by a specific mesh makes the coupling between the intra-and inter-laminar damage natural. The DPM has also been used to represent pull-through cases [START_REF] Adam | Discrete ply model of circular pullthrough test of fasteners in laminates[END_REF] where the effects of splitting on the load redistribution were correctly predicted. In addition, after a modification of the hole contour mesh, the DPM has also been able to predict notched strength, failure scenarios and failure patterns on different stacking sequences [START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF] and a first study of the scaling effect in notched composites has been performed [START_REF] Serra | Scaling effect in notched composites: The Discrete Ply Model approach[END_REF].

To the knowledge of the authors, the only numerical work on a wide range of laminate thicknesses, hole diameters and stacking sequences on an open hole tension configuration has been performed by Hallet et al. [START_REF] Hallett | An experimental and numerical investigation into the damage mechanisms in notched composites[END_REF]. There are three main differences between the approach developed by Hallet and that of the work presented here (DPM). First, in Hallet et al.'s approach, the fiber failure is simulated using a stress criterion based on Weibull statistics.

Calibration is therefore needed to extract the two parameters of the Weibull law. In the DPM approach, the fiber failure is only determined by "physical parameters": the failure strain and the critical energy release rate. Second, the locations of the matrix cracks are explicitly determined from experimental observations in [START_REF] Hallett | An experimental and numerical investigation into the damage mechanisms in notched composites[END_REF], whereas the DPM meshing strategy situates cohesive interfaces everywhere in order to have a "blind approach" and to develop a more generic model. Third, the numerical predictions [START_REF] Hallett | An experimental and numerical investigation into the damage mechanisms in notched composites[END_REF] are not able to distinguish pull-out failure from brittle failure, whereas the DPM results predict them satisfactorily.

In this paper, the DPM is used to simulate 18 configurations tested experimentally by Green et al. [START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF]. The experimental strength results are compared to numerical ones. A thorough investigation on the origin of scaling effects on the strength of notched composites is performed. The Increasing m increases the effective ply thickness, i.e. it increases the number of plies of the same orientation that are blocked together. From here on, this is referred to as ply-level scaling, and noted PLS. Increasing n increases the laminate thickness by increasing the number of sublaminates but keeps the ply thickness constant. From here on, this is refered to as sublaminate-level scaling, and noted SLS. Three different scaling routines were used: one dimensional (1D), two-dimensional (2D) and three-dimensional scaling (3D), as presented on the testing program matrix of Table 1.

t(mm)

Sublaminate-level Scaling

Hole diameter (mm)

Ply-level Scaling

Hole diameter (mm) One-dimensional scaling consists in keeping the hole diameter constant and increasing only the thickness of the laminate. In the case of two-dimensional scaling, the thickness remains constant and the in-plane dimensions (hole diameter, and hence width and length) are increased. Threedimensional scaling consists in increasing the hole diameter while increasing the thickness of the laminate. As described by Green et al. [START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF], failure was taken as the first significant load drop on the load-displacement curve (greater than 5%). Table 2 and Table 3 t (mm)

Sublaminate-level Scaling

Hole diameter (mm)

Ply-level Scaling

Hole diameter (mm) The variation of displacement, named C between point A and B is also noted and allows the average strain, , to be calculated as:

(2) = = with ( the displacement of the point A (B) and AB the initial distance between A and B. The ) initial distance between the points A and B was chosen equal to 5D in order to evaluate the effect of the hole.

The numerical displacement speed was selected so as to ensure that the corresponding additional energy was negligible, while keeping a reasonable calculation time by using multiple CPUs in parallel. A total calculation time of 5 ms was used for the ''Small" specimens and of 10 ms for the larger sizes. To ensure that changes in the calculation time and consequent increases in the inertial forces do not alter the solution significantly, a calculation was performed with a numerical displacement speed 4 times less than the initial one and showed similar results. Computation time varied between 30 min (36 CPUs) for the ''Small" specimens with about 100,000 elements and 60 h (72 CPUs) for the largest specimens with about 5,000,000 elements. The DPM uses a small number of parameters that are all provided by experimental tests [START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[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][START_REF] Adam | Discrete ply model of circular pullthrough test of fasteners in laminates[END_REF][START_REF] Achard | Discrete ply modelling of open hole tensile tests[END_REF][START_REF] Serra | Scaling effect in notched composites: The Discrete Ply Model approach[END_REF]; no coupling parameter is required.

Modeling of matrix cracking

As presented in Figure 4, the matrix cracking is taken into account using vertical zero-thickness interface elements between 2 consecutive volume elements. The failure criterion is calculated in the volume elements neighboring the interface elements, to avoid stress concentration. The stiffness of the interface (initially 10 6 MPa/mm) is set to zero if the Hashin criterion (1) is reached in either of the neighboring volume elements:

(3)

[ 9 : + 9 < : ] 2 + > 2 ? + > 2 @ (> < ? ) 2 B 1
where is the transverse stress, and are the shear stresses in the (lt) and (tz) planes, is > ? > @ < the transverse failure stress and is the shear failure stress of the ply. > < ?

When interface elements are used to simulate matrix cracking, as only through-the-ply cracks are taken into account in this method, diffuse damage is not modeled. The supporting hypothesis is that in a thin ply, a matrix crack propagates fast and therefore diffuse damage should not be represented.

Modeling of delamination

Delamination is modeled using cohesive interfaces like matrix cracking (Figure 4). Unlike for the matrix cracking, in the case of delamination, energy is dissipated, since, if we were to use an energy based criterion for matrix cracking, the energy dissipated would be dependent on mesh density. However, energy spent in matrix cracking is still taken into account: when the critical energy release rate for delamination is measured, it encompasses delamination energy and most of the energy dissipated through matrix cracking as the classical tests used to evaluate the fracture toughness (Double Cantilever Beam or End Notched Flexure tests) include several types of damage.

The energy dissipated through matrix cracking is therefore not taken into account in the matrix crack behavior law but rather in the delamination one.

For the present discrete modelling, even if there is no parameter coupling delamination and matrix cracking, the discontinuity still allows this interaction to take place. Delamination normally occurs between differently oriented plies. It was therefore simulated in interface elements, and joining nodes of lower and upper volume ply elements (Figure 5). Mode I was in the thickness direction normal to the delamination plane, while mode II and mode III were in the in-plane direction. A mode I delamination is initiated when the distance between nodes reaches a critical distance , then degradation of the stiffness in mode I is established, and ends when the

Modeling of fiber failure

Fiber failure was modeled continuously as described in Figure 9: failure is initiated when the fiber strain reaches the failure strain and is propagated using the critical energy release rate . The

(FGD < C )
affected volume elements is gradually reduced until fiber strain reaches :

E 1 (4) H I ( H E 1 0 ? J K ?) J KI = FGD < C
where is the longitudinal stress (strain), V(S) is the volume (section) of the element, is the

? ( ? ) E 1
strain of total degradation of the fiber stiffness (Figure 9), and is the critical energy release rate in D < C opening mode in the fiber direction. First order volume elements (C3D8) were chosen, as two nodes along the thickness direction of the ply are sufficient to describe linear strain variation within the element.

In order to be able to dissipate the critical energy released due to fiber failure fracture per unit area of crack, a communication law between the 8 integration points was used within each element to determine the average energy release restitution rate (over 8 integration points) and the different (with i the integration point number). The black line in Figure 9 represents the / # 1 stress/strain curve evolution of the first integration point that reaches . The 7 grey lines / # 0 display the stress/strain curve evolution of the other 7 integration points. Independently of the The images presented in Figure 18 were taken at structural failure for each of the three configurations (represented by red triangles in Figure 17).

The stress at which delamination occurred increased from 300 MPa for a 4 mm thick laminate with a 3.175 mm hole diameter, to 425 MPa for a 4 mm thick laminate with a 25.4 mm hole diameter.

As specimens fail by delamination, there is a discontinuity on stress-strain curves, and it can be seen in Figure 17 that load drops depend on the size of the hole. The smaller the hole is, the earlier the load drops take place. Load drops correspond to massive delaminations before 0° fiber failure but the stress value in the first load drop is always taken as failure even if it is not the maximum value of stress on the curve, as the specimen has delaminated extensively.

PLS [45m/90m/-45m/0m]S

[454/904/-454/04]S

mm in PLS configuration

Figure 18 shows the difference between delaminated elements according to the size of the hole.

For a large hole, delamination is small, there is just initiation of splitting on the 0° ply, and matrix cracking just occurs around the hole boundary and free edges. For a smaller hole, delaminations are larger and accompanied by splitting, forming the extensive delamination. The extensive delamination on a smaller hole propagates quickly from the hole to the specimen edge and quickly reaches failure

Sublaminate level scaling

In specimens with thickness increased by the sublaminate-level scaling method, fiber failure occurs first, with extensive pull-out on smaller specimens, as for the 2 mm thick laminate with For pull-out type failure, the laminate failure is sudden and dominated by fiber failure, which triggers delamination while propagating. Figure 20 represents a typical stress-strain curve for a 2 mm thick sublaminate-level scaled specimen with a 3.175 mm diameter hole. It is seen that the stressstrain curve increases linearly and ends suddenly when the fiber failure occurs. In the pull-out failure mechanism, fiber failure within the 0° plies is accompanied by delamination between some of the off-axis plies. Depending on the thickness of the laminate, this delamination occurs on plies close to the surface or across the entire thickness of the specimen (Figure 20-a). Matrix cracking is observed in several plies, essentially on 90° plies, crossing from the hole edge to free edges (Figure 20-b).

Matrix cracking on off-axis plies is just located around the hole, in the initiation phase.

Brittle failure has a mechanism similar to that of pull-out type failure; it is fiber-dominated.

Figure 21 represents a typical stress-strain curve for a 4 mm thick sublaminate-level scaled specimen with a 12.7 mm diameter hole. The stress-strain curve is seen to increase linearly and ends when the fiber failure occurs. The failure begins from matrix cracking with delamination in all plies and then fiber failure in the 0° ply appears as shown in Figure 21. For fiber failure and delamination between -45°/0° ply, failure occurs on the edge of the hole, and then propagates adjacent to the hole, creating a clean fracture surface across the width of the specimen, at the mid-plane of the laminate. The fracture plane is the same in every ply of the laminate. Matrix cracking begins around the hole and then extends only on the free edges, due to the 90° plies (Figure 22). The failure scenario is similar to that of pull-out failure but with less damage (matrix cracks and delamination) extension before total failure. The brittle type of failure is correctly accounted for by the model since numerical failure patterns correlate well with exprimental ones (Figure 3

-Top right).

Increasing the laminate thickness by using extra sublaminates, from 4 mm to 8 mm, does not change the failure stress value significantly, as seen the 1D scaling results (Figure 16). This means that strength is globally independent of thickness. Increasing the thickness of the specimens by increasing the number of sublaminates results in distribution of 0° plies throughout the thickness. The curves obtained experimentally and numerically are plotted in Figure 12, Figure 16, and Figure 19, respectively, for 1D, 2D and 3D scaling. Good experiment-model correlation can be observed. However, the numerical results exhibit slightly higher failure stresses than the associated experimental results. The 1D scaling curve represents the variation of the strength as a function of the thickness. It can be seen that, when the total thickness of the laminate increases, the failure stress is much higher for SLS configurations than for PLS ones. The effect of ply thickness is much stronger in PLS than in SLS because the failure mechanism of laminates in PLS is delamination, while SLS laminates fail by pull-out or by brittle-type failure. Increasing the ply thickness in PLS laminates leads to fast delamination and this phenomenon affects the strength of the structure negatively. The 2D scaling curve is less straightforward. The reader will recognize a classical sizeeffect curve for the SLS configurations whereas the PLS ones exhibit a counter-intuitive increase of the failure stress when the hole diameter increases. This is due to the grouping of plies, which favors delamination. Indeed, delamination saves the specimen by delaying fiber failure. In 3D scaling, both configurations show a decrease in failure strength when the hole diameter increases. Thus, for PLS, the effect of increasing total thickness is dominant over the effect of increasing hole diameter. There is no significant difference between PLS and SLS, so, in the end, the classical size effect is seen: increasing the diameter of the hole leads to a decrease in the strength.

The final fitting curves for the SLS and PLS configurations are presented in Figure 23. The black curve represents the baseline, which is both PLS and SLS. The blue and black curves represent the SLS specimens, in which the main failure mechanisms are pull-out failure or brittle type failure. The effect of the hole diameter on the strength is quite logical: increasing hole diameter leads to a decrease in strength. Whatever the thickness of the laminate or the diameter of the hole, the strength decreases with the thickness and the size of the hole. On the other hand, the PLS variation of strength is less straightfroward. This is due to the fact that failure is predominantly influenced by delamination. It is also seen that PLS specimens have a different behavior below and above a thickness of 2 mm. The PLS effect follows the same logic as the SLS for thicknesses from 1 mm to 2 mm. Then, when the thickness increases from 2 mm onwards, the PLS effect is positive: increasing hole diameter leads to an increase of strength. This is due to the grouping of plies, which favors delamination, then delamination saves the specimen and delays the fiber failure. This thickness effect is beneficial at the beginning for large diameters but becomes catastrophic, with a very large drop of strength, when too many plies are grouped. This can be explained by the beneficial effect of grouping some plies together to promote delamination. However, the use of this effect should not be exaggerated: if the laminate is too thick, early delamination occurs and this is detrimental.

  Composite materials have become indispensable in the aerospace and astronautics sectors.

  details of the experimental procedure. Then, the DPM modeling strategy is detailed and the results are compared with experimental data. The DPM has been used as a prospective tool to establish charts for engineers to help select the best stacking sequence to use for a given laminate thickness/ hole diameter combination. a numerical investigation of the experimental data from research carried out by Green et al.[START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF]. They performed open hole tensile tests on IM7/8552 carbon fiber/epoxy quasi-isotropic laminates manufactured from prepegs supplied by Hexcel. The nominal ply thickness was 0.125 mm. Laminates were autoclave cured according to Hexcel's recommended cure cycle[START_REF]Epoxy Matrix Product Datasheet[END_REF], and specimens were cut from these using a water-cooled diamond saw. The holes were drilled with tungsten carbide drill bits, with the drill speed scaled according to the hole diameter so that a constant cutting speed was used. Specimens with a central circular hole were tested in quasi-static tension. Constant length to hole diameter (L/D) and width to hole diameter (W/D) ratios were used, as shown in Figure1. Six specimens for each condition were loaded to failure. A nominal strain rate of 0.78%/ min was used, which equates to using constant displacement rates of 0.5, 1, 2 and 4 mm/min for hole diameters of 3.175, 6.35, 12.7 and 25.4 mm, respectively. Many of the tests were recorded with a video camera using a high magnification lens, so that the ply-level damage during the test could be determined.

Figure 1 .Figure 2 .

 12 Figure1. Specimen geometry[START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF]

  summarize the values of failure stresses obtained by numerical simulation by DPM (presented in the next section) and experimental results from Green et al. [1], covering both laminates with all plies blocked together and multiple sublaminates of single plies. Only the failure stresses are given, without the experimental stress/strain curves. Examination of the specimens post failure showed three distinct types of failure mechanisms [1]. These have been termed pull out (fiber-dominated failure with extensive sub-critical damage -Figure 3/top left), brittle (fiber-dominated failure with little subcritical damage -Figure 3/top right) and delamination (matrix-dominated -Figure 3/bottom) type failures. The configurations for which they occurred are shown in

Figure 3 .

 3 Figure 3. Top -Photographs of pull-out failure (left) and brittle failure (right)Bottom -Schematic diagram of delamination failure[START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF]

0 CFigure 8 Figure 8 .

 088 Figure 8. DPM delamination behavior law [33] Left -Pure Mode I, Right -Mixed mode

  lengthen until consistent energy has been dissipated. The stiffness of the

  constant energy is released per unit area. The characteristic length used in the model, which is the volume element length, ensures the mesh independency[START_REF] Bazant | Scaling of Structural Strength[END_REF],[START_REF] Bazant | Size Effect and Fracture Characteristics of Composite Laminates[END_REF].

Figure 9 .Figure 13 .Figure 14 -

 91314 Figure9. Discrete Ply Modeling of fiber failure[START_REF] Hongkarnjanakul | Modélisation numérique pour la tolérance aux dommages[END_REF]

Figure 18 .

 18 Figure 18. Delaminated elements, splitting on 0° ply and matrix cracked elements for t = 4 mm, (a) d = 3.175 mm, (b) d = 6.35 mm and (c) d = 25.4 mm in PLS configuration

3 .

 3 175 mm hole diameter shown in Figure 20, and more brittle failure occurs on larger specimens like the 4 mm thick laminate with 25.4 mm hole diameter, as shown in Figure 22 and Figure 21.

  Therefore, the 0° plies have the capacity to arrest propagation of the sub-critical damage, and confine it to the two outermost sublaminates.

  framework has been presented whereby user-defined volume and cohesive element formulations have been developed to simulate the scaling effects acting on open hole tension tests. This has been achieved by taking the three main failure modes in laminated composites into account: matrix cracking, delamination and fiber failure. A specific meshing methodology has been used to allow for natural coupling between matrix cracks and delamination openings. The Discrete Ply Model needs a small number of parameters, all of which come from material tests. No numerical tuning parameter has been used. The numerical framework has been validated against a large set of experimental data. Three types of scaling (1D, 2D and 3D) have been studied on two types of stacking sequences (Sublaminate-level and Ply-level scaling). The 3 types of failure (brittle, pull-out and delamination) observed experimentally have been correctly simulated. Consistent average failure stress values and evolutions have been observed according to various scaling methodologies. Since the numerical results correlate well with the experimental observations, it can be considered that the DPM has been validated for evaluation of the failure of open hole tensile tests. This model has therefore been used to open up the field of investigations and simulate configurations not previously tested experimentally. Numerical simulations have been performed for all laminate thicknesses, from 1 mm to 8 mm, and for all hole diameters, from 3.175 mm to 25.4 mm. The idea is to provide charts to help engineers choose the best stacking sequence configuration bewteen PLS and SLS for a given laminate thickness/hole diameter combination.Globally, the thicker the laminate is, the more suitable the SLS confiuguration seems to be. However, in some specific combinations (small thickness/small diameter hole), PLS can be a smarter choice.It should be noted that the research presented in this article draws this conclusion for a specific type of material, stacking sequence and loading. Other combinations could lead to different behaviors and conclusions. It would be interesting to use the Discrete Ply Model to simulate a large number of cases, with different stacking sequences and compare the influence of stacking sequence with the hole and thickness effects. Moreover, the laminate has been studied in pure tension and it would be of interest to see whether the conclusions hold under multiaxial loadings.

Table 1 .

 1 Testing program matrix[START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF]

		3.175	6.35	12.7	25.4	3.175	6.35	12.7	25.4
	1
	2	1D	3D
	4		2D
	8
		Configuration tested

Table 4 .

 4

	Hole Diameter (mm)

Table 2

 2

	F o r
	P e
	e r
	R e
	v
	i e
	w

. Numerical vs experimental average failure stress (in MPa) of sublaminate-level scaled specimens (exp. results from

[START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF]

)

Table 3 .

 3 Numerical

vs experimental average failure stress (in MPa) of ply-level scaled specimens (exp. results from

[START_REF] Green | An experimental investigation into the tensile strength scaling of notched composites[END_REF]

)

Table 4 . Failure mechanisms in testing program [1]

 4

	F o r
	P e
	e r
	R e
	v
	i e
	w

http://mc.manuscriptcentral.com/jcm Journal of Composite Materials

Acknowledgment

The authors would like to express their gratitude to CALMIP (CALcul en Midi-Pyrénées, (https://www.calmip.univ-toulouse.fr) for access to the HPC resources under allocation p1026.

http://mc.manuscriptcentral.com/jcm Journal of Composite Materials