I Olmedo

F Bourrier
email: franck.bourrier@irstea.fr

D Bertrand

F Berger

A Limam

Dynamic analysis of wooden rockfall protection structures subjected to impact loading using a Discrete Element model

Keywords: Rockfall, Wooden structures, Protection structures, Fresh wood, Impact, Discrete Element Method

Wooden structures made of felled trees are used in forested slopes as protection structures against rockfall. A model based on the Discrete Element Method is developed to analyse their response to normal impacts of blocks. The interaction between the blocks and the stems and the stem response are explicitly integrated. After the model calibration using impact tests, a wide range of impact scenarii representative of real configurations is explored. The influence of three parameters (stem diameter, impact velocity and block diameter) is analysed.

Three impact types were identified for different ratios RD between the block and stem diameters. For small RD , the impact is limited to a brief contact. For intermediate RD, successive contacts of similar duration are observed. After the last contact, the block is sent back with a large velocity. For large RD, two contacts are observed. The second contact can be assimilated to a quasi-static loading of the stem. In addition, structure damage increases for increasing RD ranging from almost no damage for the first impact type to rupture of the stem for the third type. The largest RD associated with the first impact type is the target configuration that favors block energy decrease and limits structure damage.

Introduction

In the field of civil engineering, numerous models have been developed to assess the efficacy of a wide range of protection structures against rockfall hazard: metallic fences [START_REF] Bertrand | Full-scale dynamic analysis of an innovative rockfall fence under impact using the discrete element method: from the local scale to the structure scale[END_REF][START_REF] Cazzani | Dynamic finite element analysis of interceptive devices for falling rocks[END_REF][START_REF] Gottardi | Full-scale modelling of falling rock protection barriers[END_REF][START_REF] Volkwein | Flexible rockfall barriers subjected to extreme loads[END_REF] or embankments [START_REF] Lambert | Design of rockfall protection embankments: a review[END_REF][START_REF] Plassiard | Rockfall impact parameters on embankments: a discreteelement method analysis[END_REF][START_REF] Ronco | Design of reinforced ground embankments used for rockfall protection[END_REF], for example. These numerical approaches are of interest as they allow accounting for the physics involved in the modeling of block impacts onto complex structures. Numerical simulations enable to analyse the structures and to propose recommendations for their design. However, experimental data remain necessary to calibrate and validate the numerical models. Given the technical difficulties and economic costs to carry out in-field real-scale impact experiments, reduced scale laboratory tests are an interesting alternative to get accurate experimental data.

In alpine forested slopes, specific wooden structures made of felled trees are increasingly installed as a complement to the protection provided by the forest stand itself [START_REF] Dorren | Mechanisms, effects and management implications of rockfall in forests[END_REF]. Contrary to the impact of blocks on standing trees, [START_REF] Dorren | Real-size experiments and 3-D simulation of rockfall on forested and non-forested slopes[END_REF][START_REF] Lundstrom | The root-soil system of Norway spruce subjected to turning moment: resistance as a function of rotation[END_REF], the impact of blocks on these wooden protection structures has not been investigated.

Recent experimental laboratory studies enabled to characterize the dynamic response of small stems (up to 10 cm of diameter) to impacts [START_REF] Bertrand | Experimental and numerical dynamic analysis of a live tree stem impacted by a Charpy pendulum[END_REF][START_REF] Olmedo | Experimental analysis of the response of fresh wood stems subjected to localized impact loading[END_REF] with the final aim of understanding the response of wooden protection structures to a block impact. A first Discrete Element Method (DEM) based model [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF] of the stems response to impact loading was also proposed and allowed describing the dynamic response of the stems tested in the laboratory. In the simulations, the loading was imposed by the experimentally measured time evolution of the force at mid-span. Thus, the proposed model did not explicitly account for the block and its contact with the impacted system.

These studies focused on the dynamic response of fresh wood stems only. In a first step, the time evolution of the structure response, which is related with the degradation of the mechanical properties of the wood material along the stem composing the structure and in the vicinity of the structures boundaries, was not investigated. Consequently, these studies can only be used to analyze the response of a numerical model of protection devices made of felled trees considering their mechanical properties after setting up the structure.

The objective of this paper is to analyse the response of real scale wooden protection structures to the normal impact of a block after setting up the structure for design purposes. To meet this objective, the DEM model proposed and calibrated in [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF] has been improved by integrating a contact model describing the local interaction between the block and the fresh wood stem. After a calibration phase, the improved model was used to simulate real scale impacts on felled trees protection structures (tree diameters ranging from 35 to 80 centimeters) for a wide range of impact conditions. Numerical predictions were used to analyse the physical processes involved for different impact configurations and to identify the most favorable structure configurations in terms of block momentum reduction as well as structure damage and rupture. These observations were finally discussed to propose general recommendations for the design of these structures.

Numerical modelling

Discrete Element model

The DEM was used to model the mechanical response of a fresh stem to a dynamic loading [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF]. This method allows modeling the dynamic motion of particles that interact by contact or remote interactions [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. At each time step, the particles' locations, orientations, geometries, and physical properties are used to calculate the sum of forces on each particle. Newton's law are then solved for all the particles, allowing the update of their positions and orientations. The DEM has been widely applied for the study of granular materials [START_REF] O'sullivan | Particle-based discrete element modeling: Geomechanics perspective[END_REF]. It has also been used to model deformable structures (e.g. [START_REF] Bertrand | Full-scale dynamic analysis of an innovative rockfall fence under impact using the discrete element method: from the local scale to the structure scale[END_REF][START_REF] Effeindzourou | Modelling of deformable structures in the general framework of the discrete element method[END_REF].

The numerical developments were implemented in the YADE-dem code [START_REF] Smilauer | The Yade project[END_REF]). In the model, the fresh wood stems are considered homogeneous cylindrical bodies composed of nodes (N ) regularly located along the stem longitudinal axis x (Figure 1). The nodes are linked by cylindrical connections (C) (Bourrier, Kneib, Chareyre, & Fourcaud, 2013) that can interact with external bodies, such as spheres. The model of the interactions between the nodes (N ), called stem model in the followings, governs the mechanical response of the deformable cylinder (i.e. a fresh wood stem). The model of the interactions of the cylindrical connections (C) with external bodies, called contact model in the followings, controls the interaction between the cylindrical element (i.e. a fresh wood stem) and an external body (i.e. a block).

Fresh wood stem model

The stem model was initially developed and calibrated in [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF]. It allows calculating the forces and moments between contiguous nodes depending on their relative positions and orientations, which determines the macroscopic response of the cylinder (i.e. stem). In the context of the study, the stems are mainly subjected to bending loads. Thus, the relation between node positions (resp. orientations) and interaction forces (resp. moments) associated with tensile, shear, or twist efforts are modeled using an elastic formulation that depends on the longitudinal modulus of elasticity (MOE) and the Poisson's ratio of the stem material.

A multi-linear relationship is used to model the response of the stem to bending loadings. The moments between contiguous nodes associated with bending loadings are the moments M y and M z (called M in the following) along the y and z axes of the local coordinate system associated with the first node of the two contiguous ones. These moments are calculated incrementally depending on the nodes relative orientations θ y and θ z (called θ) along the y and z axes (Figure 1b) and applied to both of the contiguous nodes.

A bi-linear relation describes the M (θ) loading phase (i.e. for increasing θ)(Figure 2). An elastic threshold moment (M Elast ) allows defining the transition from the first linear relation, defined by the stiffness coefficient k Elast , to the second one, associated with the coefficient k N L . The unloading phase (i.e. for decreasing θ) is also governed by a linear relation described by the coefficient k U nld . Bourrier et al. (2013) and [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF] proposed relationships that relate the stiffness coefficients k Elast , k N L , and k U nld to the stem cross section geometry, distance between the contiguous nodes, and material properties:

k Elast = M OE I l (1) k N L = E N L I l (2) k U nld = E U nld I l (3)
where M OE, E N L and E U nld are the longitudinal modulus of elasticity, the mod-ulus governing the second linear phase, and the modulus for the unloading phase, respectively. l and I are the distance between contiguous nodes and the area moment of inertia (I = πD 4 /64), respectively. D is the stem diameter at the node location.

When the bending moment M Rupt is reached, which corresponds to stem's failure, the interaction is considered broken and the simulation is stopped.

The elastic threshold (σ Elast ) and the ultimate stress before failure (σ Rupt ) are calculated from the geometrical properties, M Elast , and M Rupt such as:

σ Elast = M Elast I D 2 (4) σ Rupt = M Rupt I D 2 (5)
σ Rupt assumption, even if not valid when material non-linearities develop, is widely accepted for the mechanical characterization of wood [START_REF] Bergman | Wood handbook -Wood as engineering material[END_REF]. Indeed, it corresponds to the definition of the modulus of rupture (M OR).

A viscous damping of the bending moment M was finally implemented following an incremental formalism:

dM Damp = βk Elast d θ (6)
where the increment of the damping contribution to bending (dM Damp ) is calculated from the increment of the relative angular velocity of the nodes (d θ) and a damping coefficient (β).

Finally, the calibration parameters of the stem model are: M OE, E N L , E U nld , σ Elast , σ Rupt , and β.

Block -stem contact model

The block-stem contact model added to the initial stem model consists in modelling the interactions between the cylindrical elements and external bodies. In the DEM framework, the contact model between two particles usually relates the contact forces to the interpenetration and to the relative velocities at the contact point. The contact plane, tangent to the contact surface at the contact point, enables to define the normal and tangential components of the contact force. In the context of the study, only the normal component is considered because impacts are assumed orthogonal to the cylinder longitudinal axis.

Previous results from [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF] showed that including a simple linear contact model was sufficient to significantly improve the relevance of the impact model. However, in this study, we choose to use a contact model integrating energy dissipation at the contact which allowed us to analyse energy dissipation processes during impact.

The normal contact force (F ) depends on the indentation (α) of the impactor into the cylinder (Figure 3). From the numerical point of view, it is described by the overlap between the block and the cylindrical element. In the absence of additional information concerning the highly variable block shape, a formulation corresponding to the Hertz-Crook contact theory for spheres developing an elasto-plastic response was chosen [START_REF] Abrate | Impact on composite structures[END_REF][START_REF] Crook | A Study of Some Impacts between Metal Bodies by a Piezo-Electric Method[END_REF]. During the loading phase, the contact force can be expressed as:

F Load = k c .α 3/2 (7)
where k c is the contact stiffness, which can be calculated from the physical and geometrical properties of the colliding bodies (bodies 1 and 2):

k c = 4 3 E R √ R (8)
where,

1 R = 1 R 1 + 1 R 2 (9) 1 E R = 1 -ν 2 1 E R1 + 1 -ν 2 2 E R2 (10)
R j , E Rj and ν j are the radius, radial modulus of elasticity, and Poisson's ratio of the body j.

For the unloading phase, the formulation proposed by Crook [START_REF] Abrate | Impact on composite structures[END_REF][START_REF] Crook | A Study of Some Impacts between Metal Bodies by a Piezo-Electric Method[END_REF], assuming negligibles permanent indentations, was implemented:

F U nld = F m ( α α m ) 5/2 (11)
where F m is the maximum force reached during the loading phase and α m is the corresponding indentation.

Finally, for this study, the calibration parameters of the contact model are E R , and ν of each colliding body. et al. (2016) proposed a calibration of the stem model where the contact model was not taken into account, i.e. the loading was modelled by a time varying force measured experimentally. In this paper, a new calibration is proposed that accounts for the contact model to describe the interaction between the block and the stem.

Model calibration

Olmedo

Impact experiments

The chosen species for the tests was European beech (Fagus sylvatica) for being representative of the French Alps forests exposed to rockfall (i.e. European beech is present in the 25% of the french protection forests) (IGN, 2014). The wood pieces were extracted from the same forest patch at the Col de Vence forest (Rhone Alpes, France). Tested trees had similar growing conditions and were cut at the same period of the year to reduce mechanical properties variations. Stems were tested maximum one day after their cut to limit changes in wood moisture content.

The laboratory experiments on fresh wood stems consisted in impacting simply supported stems at mid-span with a Mouton-Charpy pendulum. Stems of diameters between 3 and 9 centimeters and of 1.5 m length between supports were impacted. Force cells (Siscam MC-100) enabled to measure the evolution in time of the impact force and of the reaction forces at the supports. A high speed camera (Photron Fastcam, Nikkor 35-80 lens) was installed to measure the impactor and stem's mid-span displacements. Each stem was impacted twice, at two different impact energy levels. Additional details on the experimental set-up can be found in [START_REF] Olmedo | Experimental analysis of the response of fresh wood stems subjected to localized impact loading[END_REF].

Calibration procedure

The model calibration was based on experimental results from [START_REF] Olmedo | Experimental analysis of the response of fresh wood stems subjected to localized impact loading[END_REF]. The stems were modeled by deformable cylinders and the time evolution of the measured impact force was imposed at mid-span. The calibration was done by adjusting the calibration parameters to achieve a good fit between experimental and numerical impact work calculations (W Ext ). The external work is defined such as:

W Ext = 0 t F dy b ( 12
)
where F is the force applied onto the stem and y b is the impactor displacement at the instant t [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF]. Thirty-seven impact tests were calibrated, which allowed to estimate the mechanical properties of the tested wood.

The calibration focused on the adjustment of the parameters concerning the stem (M OE, E N L , E U nld , σ Elast , σ Rupt , β).

For the pendulum, classical values for steel were used as parameters for the contact model : E R1 = 210 GP a and ν = 0.3. For the stem, the parameter E R2 of the contact model was set at 160 M P a, as proposed by [START_REF] Quetel | Approche multi-chelle du comportement mcanique d'un arbre soumis l'impact d'un bloc rocheux[END_REF], and the Poisson's ratio ν was set at 0.3 [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF].

The experimental data from [START_REF] Olmedo | Experimental analysis of the response of fresh wood stems subjected to localized impact loading[END_REF] allowed adjusting the parameters of the stem model. The simulations consisted in impacting at mid-span a simply supported stem, represented by a cylinder, with a spherical discrete element, representing the pendulum. The measured impact velocity of the pendulum was imposed to the sphere at the beginning of the simulation. The parameters values obtained from the previous calibration of the stem model [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF] were initially used. The calibration consisted in slight adjustments of some of these parameters to obtain the best agreement with experimental results. The numerical simulations were adjusted to fit the experimental results in terms of energy, W Ext .

Impact scenarii

Description

After the calibration phase, the model was used to simulate blocks impacts on felled trees made structures corresponding to real scenarii in terms of field impact conditions. The analysis focused on the impact response of a simplified structure made of one tree. The structure was modeled as a deformable cylinder with simply supported boundary conditions (Figure 4), which is representative of real set-ups. Impacting blocks were considered as rigid spherical bodies with a translational velocity perpendicular to the longitudinal axis of the tree, and the block angular velocity was considered nil. Structures were impacted at mid-span and impacts were centered on tree's longitudinal axis.

For the design of structures, engineers and practitioners consider the characteristics of the expected rockfall events (in terms of block mass and velocities) to propose an ad-hoc rockfall protection device capable to resist to the incoming impact energy. In particular, stem diameter needs to be defined for the setup of the structure.

Simulations were performed for several impact scenarii, characterized by the block diameter (D b ) and impact velocity (v in ). Moreover, different stem diameters (D t ) were considered. Thus, D t , D b and v in were considered as input variables. The stem length L was set at 10m for all simulations.

Concerning the impacting block, D b and v in ranges are representative of rockfall events Volkwein et al. (2011)(Tab.1). Tree diameters (D t ) ranges are representative of the current felled trees made protection devices installed in the European Alps (Tab.1). When Db Dt > 4, it is assumed that the stem cannot absorb a significant amount of block energy [START_REF] Dorren | Mechanisms, effects and management implications of rockfall in forests[END_REF]. The sets of parameters presenting Db Dt > 4 were thus not considered. 1386 scenarii, corresponding to different values of D b , D t , and v in , were finally considered for the simulations.

Mechanical properties of fresh wood for in-field scale simulations

The fresh wood properties depend on many factors and the characteristics of each forest stand cannot be directly related to the experimental data presented in this paper. Despite this large variability, the fresh wood properties were assumed constant for all simulations. Given the very reduced amount of data in the literature concerning the mechanical properties of fresh wood at the tree scale, the choice was made to set up these mechanical properties depending on the results obtained during the model calibration.

The parameters E R and ν describing the contact model were set at the same values as in the calibration phase. The block's and stem's densities (ρ b and ρ t ) were set at 2800 kg.m -3 and 1200 kg.m -3 , respectively.

For the stem mechanical properties, the previous model calibration [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF] provided a characterization of fresh wood mechanical properties for small diameter trees (D t < 10 cm). Given the complexity of the tree's growing mechanisms, the mechanical properties of wood evolve along the tree [START_REF] Keith | Annual layers effect resistance of wood to impact[END_REF] and radially (Lundstrom et al., 2007). In addition, the irregularities associated with growing conditions, branch insertions, or scars lead to significant variations on the mechanical properties of wood [START_REF] Lundstrom | Analysing the mechanical performance and growth adaptation of Norway spruce using a non-linear finite-element model and experimental data[END_REF]. Consequently, the previous model calibration for small diameter trees (D t < 10 cm) cannot be directly extrapolated to larger diameter trees (D t > 35 cm).

The values obtained have been compared to previous studies found in the literature concerning large diameter stems ranging between 20 cm and 50 cm. [START_REF] Lavers | The strength properties of timber[END_REF] proposed values of MOE and modulus of rupture (M OR = σ Rupt = 32MRupt πD 3 ) [START_REF] Bodig | Discrete modeling of granular soils reinforcement by plant roots[END_REF]) for fresh European beech stems with a 50 % moisture content. They estimated the MOE at 9790 MPa and the MOR at 64.8 MPa from quasi-static bending tests. The differences between the previous model calibration values [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF] and [START_REF] Lavers | The strength properties of timber[END_REF] estimations are not significant in view of the variability of wood properties [START_REF] Bergman | Wood handbook -Wood as engineering material[END_REF]. The MOE measured by [START_REF] Lavers | The strength properties of timber[END_REF] is 9.6% larger than the calibration value. The MOR measured by [START_REF] Lavers | The strength properties of timber[END_REF] is 10.3% lower. Provided that the values from the calibration described in section 2.2.2 do not differ substantially from these obtained in [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF], this comparative analysis will lead to use the values obtained from the calibration of the model that explicitly account for block-stem contact for the simulations for in-field structures made of large diameter trees.

Energy transfers

The main goal of the design of a rockfall protection device is to maximize the block's energy loss due to the impact. In addition, the study of the energy transfers and dissipations on the protection structure is mandatory to develop optimal and adapted protection devices and, at the same time, guarantee their durability.

During the impact of a block on a protection device, an energy exchange between the block and the structure occurs. The external work related to the forces applied onto the block at a given time t (W Ext ) is equal to the block's kinetic energy variation from the impact beginning

(∆E b k = E k b t -E k b in )
, where E k b in and E k b t are the block's kinetic energy before impact and at time t, respectively.

W Ext = -∆E b k (13)
where W Ext is calculated from equation 12. As regards to the protective structure, the energy balance can be expressed as follows for any time during impact:

W Ext = W Int + W Cont + ∆E t k + W Damp ( 14
)
where W Int , W Cont , ∆E t k , and W Damp are the stem's strain energy, the contact strain energy, the stem's kinetic energy variation from the impact beginning, and the energy dissipated by damping, respectively.

W Int is calculated as follows:

dW Int = n-1 a=1 N a .du a + n-1 a=1 T a .dv a + n-1 a=1 M a .dθ z,a (15)
where a is the node interaction identifier and n is the number of nodes. N , T and M are the normal and tangential components of the interaction force, and bending moment between nodes participating to the interaction a. u a and v a are the normal and tangential components of the relative displacements between nodes taking part to the interaction.

The stem kinetic energy variation (∆E t k ) is calculated from the nodes masses and velocities (translational and angular) [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF].

The contact strain energy W Cont is assessed incrementally. The contact energy increment dW Cont is calculated as a function of the impact force F and of the indentation increment dα:

dW Cont = F.dα (16)
The energy dissipated by damping W Damp is also assessed incrementally using the increment dW Damp calculated from θ and the damping contribution to the bending moment (M Damp ).

dW Damp = M Damp .dθ (17)

Results and discussion

Model calibration

The results obtained from the model calibration are compared to those obtained from simulations which do not consider explicitly the block-stem interaction [START_REF] Olmedo | Discrete element model of fresh wood's dynamic response to impact[END_REF], and to the experimental data. Figure 5 depicts the time evolutions of the external work (Eq.12) for a typical situation. The integration of the contact model allows a better description of the experimental observations, especially at the early beginning of the contact and around the maximum mid-span displacement. Similar values of the parameters M OE, σ Elast , and σ Rupt are obtained for model calibrations that explicitly integrate or not block-stem contact (Tab.2).

The explicit modeling of the block-stem contact leads to reduce the damping parameter β (Tab.2). Indeed, the contact model tends to limit vibrations with high frequencies. Thus, if the stem/block contact is not explicitly taken into account, high frequency vibrations have to be damped with larger values of β parameter.

Differences between calibration values of E N L and E U nld are also observed. E N L and E U nld are the leading parameters for the modelling of correct values of the maximum stem deflection and contact force. And, near the maximum stem deflection, different structural responses occur if the block-stem contact is explictly modelled or not.

Impact scenarii

Impact classification

3.2.1.1. Impact duration (t c ). t c is defined from the instant when the block-stem contact is established to the end of the last contact. The evolution of t c as a function of the ratio R D (R D = Db Dt -Figure 6) shows that low diameter ratios (R D < 1.5) are associated with short contact durations (t c < 30 ms). In these cases, t c presents a small variability (from 10 ms to 30 ms). For R D values between 1.2 and 2, t c presents a larger variability, ranging from 70 ms to 150 ms, with no clear dependency to R D . Finally, for R D > 2 and t c > 150 ms, a linear increase of t c with R D is observed. This analysis tends to show a potential classification of the simulations in three groups (A, B, and C) presented in Figure 6.

The time evolutions of the impact forces between the block and the stem are depicted in Figure 7 for groups A, B, and C. For the group A, the block-stem interaction presents only one contact. For the group B, successive contacts (2 or 3) of similar durations can be identified. t c values are significantly higher for simulations presenting three contact phases, even if the latest is not significant in terms of force magnitude. Finally, group C simulations always present two contacts, the second being substantially longer than the first one. These clear differences in terms of time evolutions of the impact forces between the groups lead to classify each impact simulation.

3.2.1.2. Momentum variation (∆p). The block's momentum variation ∆p (∆p = m b (v in -v out ), with v out the block's velocity after impact) is calculated to evaluate the variation in the block kinematics. The relation between the block's momentum variation ∆p and the diameter ratio R D is analysed. The normalized value ∆p = ∆p pin is defined, where p in is the block momentum just before the impact (p in = m b v in ).

The ∆p values for all simulations range between two particular cases. The first one is an elastic impact implying that v out is opposite to the incoming velocity (v out = -v in ). In the second one, the block is stopped (v out = 0). For the first case, ∆p = 2 whereas, for the second one, ∆p = 1. Thus, the lower ∆p, the larger the loss of block's momentum.

Figure 8 underlines that the groups A, B and C present different trends in the evolution of ∆p with R D . Group A results present a linear decrease in ∆p with increasing R D values. Group B results globally present higher ∆p values and larger variability than group A data for a given R D value. Group C results, observed for large R D , exhibit similar values as for group B with a slightly smaller variability.

The R D values delimiting the groups A, B, and C depend on the impact velocity (Figure 9): for impacts at v in = 5 m.s -1 (resp. v in = 30 m.s -1 ), the transition from group A to B is observed at

R AB D ≃ 1.5 (resp. R AB D ≃ 1.2). The transition from group B to C appears at R BC D ≃ 2. (resp. R BC D ≃ 1.55) for impacts at v in = 5 m.s -1 (resp. v in = 30 m.s -1 ).
3.2.1.3. Structure damage. The maximum bending moment M M ax achieved during each impact simulation is compared with M Elast and M Rupt for each simulation, allowing the assessment of stem damage for each case. If M M ax < M Elast , no material non-linearities developed and no structure damage occurred. When M Elast < M M ax < M Rupt , material non-linearities appeared but the structure did not break. Finally, if M M ax > M Rupt , the rupture of the structure occurred. These three scenarii are called, zone 1, 2, and 3.

Stem rupture (zone 3) is observed for groups B and C, specially for high impact velocities. On the contrary, for most of group A simulations, the struture is not damaged at all (zone 1), only few high speed (i.e. v in > 20 m.s -1 ) impacts present material non-linearities (zone 2) (Figure 10).

Relative displacements between the block and the stem

Three impact configurations representative of each response group are analysed to highlight the differences in terms of relative displacements between the block and the stem. The input parameters for each example case are summarized in Table 3. They correspond to impacts at 5 m.s -1 .

The time evolutions of the block center of mass (y b ) and the stem's mid-span (y t ) displacements and velocities are presented in Figure 11. Block and stem are in contact when α = y b -y t 0 (c.f. red surfaces in Figure 11. Specific attention is paid to the block and stem velocities just after the first contact and at the beginning of the second contact.

As observed previously, group A impacts present only one and short contact (Figure 11a,d). At the end of this short contact, the block velocity is negative. This value increases toward positive values for increasing R D . Thus, the block is sent back more and more slowly for increasing R D which favors the occurrence of a second contact (Figure 11b,e). Group B impacts are observed if the block velocity after the first contact is positive. In addition, the second contacts occur for negative values of the stem velocity, i.e. during the stem way back from its initial oscillation.

For Group C impacts, i.e. for larger R D values, the block velocity after the first contact is positive and keeps increasing for increasing R D values (Figure 11f). Thus, the second contact occurs earlier in the stem oscillation process. It takes place for positive stem mid-span velocities. In these configurations, the block and the stem velocities directions are similar after the first contact. In consequence, the stem's midspan and the block displacements follow analogous evolutions after the first impact (Figure 11c) and till the end of the impact.

Energy exchanges

The energy balance during impact is analysed for the three cases representatives of each group.

For group A, the block's energy is nil at the end of the initial contact loading stage (i.e. t ≃ 6 ms) and most of the energy is transferred into energy at the contact (W Cont ). The contact energy is transferred back from the contact to the block and the stem during the unloading phase of the initial contact (Figure 12a). As the interaction between the block and the stem is limited to the initial contact, the block's energy loss is mainly due to momentum transfer during the first contact and, possibly, to energy dissipation at the contact (up to 20% of the impact energy in some cases).

For group B, energy transfers during the first contact are very similar to those observed for the group A (Figure 12b). The momentum transfer during the first impact induces that the block energy is almost nil at the end of the first contact (i.e. t ≃20 ms). On the contrary, the stem energy is significant between the first and second contact. It is periodically converted from kinetic to strain energy following the stem oscillation and substantial dissipation of energy occurs due to damping between the two contacts. During the second contact (i.e. 70 ms t 90 ms), a significant amount of energy is transferred back from the stem to the block. When the second contact takes place, the stem kinetic energy is high and decreases to almost zero during the second contact. Energy dissipation may occur at the contact during the first contact, but not significantly during the second one.

For the group C, as for the groups A and B, the first contact leads to a significant momentum transfer. Similarly to group B, given the high stem kinetic energy after the first contact, the energy dissipation by damping between the two contacts is significant (between 10% and 20% of the impact energy). But, contrary to group B, when the second contact occurs, the stem strain energy is high and its kinetic energy is low (Figure 12c). Thus, the energy exchange during the second contact is driven by a quasi-static interaction between the stem and the block.

The block energy loss is due, for Group A, to a significant momentum transfer and to energy dissipation at the contact during the first contact. The occurence of a second contact leads to a second energy transfer from the stem to the block for groups B and C. For these groups, a significant energy amount is also dissipated by damping of the stem motion between the two contacts. For group B, the second contact induces a momentum exchange from the stem to the block because the stem kinetic energy is large and the block one is small at the beginning of the second contact. On the contrary for group C, the stem strain energy is large at the beginning of the second contact and, thus, the second contact can be assimilated to a quasi-static interaction.

Towards design recommendations

The validity of the model is limited to the impact conditions of the study in terms of block and stem properties as well as loading type, i.e. normal impact of a block far from the boundary of the structure, and structure size, in particular. As the model was not calibrated directly from real scale impact experiments, it seems more reasonable to use the results for a qualitative analysis rather than for quantative assessments. Despite these limitations, relevant information is provided for the design of wooden protection structures.

The design of rockfall protection structures is generally composed of two steps related to the structural and functional design. The functional design is dedicated to assess the capacity of the structure to stop the blocks and the structural design aims at quantifying the capacity of the structure to sustain the loading by the block [START_REF] Lambert | Design of rockfall protection embankments: a review[END_REF].

The functional capacity of the structure relies on the analysis of the block momentum reduction while the structural design is based on the study of structure damage. Concerning block momentum reduction, the choice was made to consider the stopping of the rock (∆p = 1) as a target to assess the functional efficacy of the structure. On the contrary, a substantial velocity (∆p = 2) is considered as the worst case in terms of functional efficacy among the configurations observed. For the structure damage, the optimal design corresponds to no damage of the structure.

A design approach can be proposed with the objective of favoring one or several of the three impact groups identified. Both the functional and structural efficacies of the structure can be estimated for the different groups. For group A, the block momentum reduction depends on a single and brief interaction. The momentum reduction increases for increasing diameter ratio to reach a maximum value among all groups observed for R D around 1.3. In addition, almost no structure damage is observed. On the contary, the momentum reduction is substantially larger for both group B and C while significant structure damage may occur. Thus, Group A impacts correspond to the target configuration both in terms of functional and structural designs.

In practice, as the block volume, i.e. diameter, is given for a study site, designing the structure to reach the target configuration consists in setting the tree diameter of the structure to obtain the target range of R D . The analysis of the simulations allowed identifying optimal R D values. However, the results also show that the transition from an optimal configuration (group A) to an unfavorable one can occur for small variations of R D . The choice of the tree diameter should thus be made with caution.

Finally, the analysis of the physical processes shows that the optimal configuration corresponds to a single and brief interaction between the block and the structure. Technical solutions can be prescribed to favor the occurrence of a single impact, i.e. to avoid the occurrence of a second momentum exchange leading to a group B impact. For this purpose, changes in the boundary conditions of the structure, such as fixing the structure to loose cables anchored on stumps instead of placing it directly on the stump, can be envisaged for example.

Conclusions and perspectives

The DEM model developed integrates the description of the physical processes both at the stem and contact scales. After a calibration phase, the combined model, accounting for the contact between the block and the stem, has been used for real scale simulations.

Three types of impact have been identified depending on the ratio R D between the block radius and the stem diameter. A unique and brief contact occurs for small R D values. For intermediate R D values, at least two successive contacts of similar duration occur. The block continues in the same direction after the first contact and is sent back with a large velocity after the second one. For large R D values, two contacts are also observed. During the second contact, longer than the first one, a quasi-static loading of the stem occurs. Structure damage increases for increasing R D values : the first impact type entails almost no damage whereas substantial damage occur for the third type.

The optimal configuration in terms of structure design is a unique and brief contact between the block and the structure with a substantial decrease in the block momentum and almost no structure damage. This configuration corresponds to the largest R D values associated with the first impact type.

The model proposed is a computationally efficient tool usable for the study of the efficacy of wooden protection structures for a wide range of impact conditions and structure configurations. Further model developments are foreseen to simulate noncentered and oblique impacts and to integrate the interaction of the structure with the ground that may represent a major source of energy dissipation.

In addition, the relevance of the calibration of the model, done using experimental results on smalll diameters stems (D t < 10 cm), for large tree diameters has been assessed based on comparisons to previous studies found in the literature only. This aspect of the development of the model could be improved by detailed studies of the mechanical properties of large diameter trees. These studies could be based on experimental studies found in the literature and dedicated to the assesment of the mechanical properties of trees at the real-scale [START_REF] Ruel | Relating mechanical strength at the stem level to values obtained from defect-free wood samples[END_REF][START_REF] Schonenberger | Effect of timber removal from windthrow slopes on the risk of snow avalanches and rockfall[END_REF] or to field experiments carried out at the tree scale (Lundstrom et al., 2007).

Finally, a major improvement of the model would consist in the integration of an appropriate mechanical description of decayed wood into the model to estimate the time evolution of the protective capacity of these structures. This improvement would concern both the model of the stem mechanical properties and of the structure boundaries. Indeed, both of these components of the wooden structures are affected by wood degradation [START_REF] Benfratello | Testing mechanical characteristics of chestnut stakes used in bed sills for stream restoration[END_REF][START_REF] Romano | Residual life and degradation assessment of wood elements used in soil bioengineering structures for slope protection[END_REF]. For that purpose, detailled studies have to be carried out to characterize the time evolution of the mechanical properties of the stems and of the structure boundaries due to degradation processes (i.e. fungi and entomofaune aggression in particular).

Figure 1 .

 1 Figure 1.: Yade-DEM elements are composed of nodes and connections (a). The forces and moments applied to contiguous nodes are calculated from their relative positions and orientations (b).

Figure 2 .

 2 Figure 2.: Relation for the calculation of bending moments between contiguous nodes.

Figure 3 .

 3 Figure 3.: Block-stem contact : the indentation α is calculated from the overlap between the bodies.

Figure 4 .

 4 Figure 4.: DEM model accounting for the block-stem interaction. Face (a) and side (b) views.

Figure 5 .

 5 Figure5.: Time evolutions of the external works for the experimental results (red) and numerical predictions with imposed loading force on the stem (blue) and explicit block-stem interaction (black) for a typical impact test.

Figure 6 .

 6 Figure 6.: Contact durations (t c ) vs. diameter ratio (R D = D b /D t ). Identification of the groups A (red), B (blue) , and C (green).

Figure 7 .

 7 Figure 7.: Time evolutions of the contact force for three typical simulations belonging to the groups A (red), B (blue) and C (green).

Figure 8 .

 8 Figure 8.: Normalized momentum variation (∆p) vs. R D . The groups are depicted in red (A), blue (B) and green (C). The dashed lines show the cases of a block stopped (v out = 0) and of no block energy loss (v out = -v in ).

Figure 9 .

 9 Figure 9.: Evolution of ∆p with R D emphasizing impact velocities of 5 m.s -1 (red) and 30 m.s -1 (blue).

Figure 10 .

 10 Figure 10.: Evolution of ∆p with R D emphasizing the elastic (z1), non-linearities occurrence (z2) and rupture (z3) zones.

Figure 11 .

 11 Figure 11.: Top: time evolutions of the mid-span (dashed line) and block (solid) displacements for group A (a), B (b) and C (c) simulations ; Bottom: time evolutions of the stem's mid-spam (dashed line) and block's (solid) velocities normalized by the impact velocity v in for group A (d), B (e) and C (f) simulations. Red surfaces indicate the block-stem contacts (α > 0).

Figure 12 .

 12 Figure 12.: Time evolutions of block and stem energies normalized by Ek b in : kinetic (-∆E b k and -∆E t k ), bending strain (W B ), damping dissipation (W Damp ) and contact (W Cont ) energies. Shear and normal strain energies are negligible and, thus, not represented on the graphs.

Acknowledgements

This work was supported by the Natural Hazards and Vulnerability of Structures (VOR) research network. This network is funded by the French Ministry of Research and joins different laboratories in the Rhônes-Alpes region.


Laurent Robert

Gilles Dusserre

Anisothermal Thermosetting Resin Cure Monitored by Optical Fiber Refractometer

The cure monitoring of composite parts is an important key to optimize the industrial curing cycles. In this study, it is proposed to monitor the cure of a thermosetting resin by a Fresnel refractometer. The measurements are first compared to dielectric analyses and predictions of a thermo-kinetic model that uses parameters from the literature. Under isothermal conditions, the results show that the Fresnel power can be directly exploited to obtain qualitatively the degree of conversion. This is also demonstrated in anisothermal conditions by correcting the Fresnel signal from its temperature dependence. If the initial degree of conversion is known, its quantitative assessment is achieved using the Fresnel and temperature measurements.

INTRODUCTION

Thermosetting-based composite materials are widely used as structural materials because of their good ratio between weight and strength. However, some well-known issues concern damage [1], residual stresses [2,[START_REF] Mulle | [END_REF], or mechanical properties of these materials [4]. These issues meet some obstacles partly related to the variability of the properties induced by the manufacturing processes. Indeed, the quality of the composite part depends on the control of the process cycle and its parameters (temperature, pressure/vacuum, curing time). In order to improve the curing of thermosetting resins and therefore the quality and reliability of these composite parts, high performance in situ sensing systems are required for the monitoring of processes [5][6][7][8][9][10][11][12][13][14][15][16].

Many different techniques have been developed to characterize the cure process of thermosetting resins. One can cite (see, for example, Ref. 15] and references within) thermal measurements (such as differential scanning calorimetry (DSC) and thermal mechanical analysis (TMA)), dielectric analysis (DEA), nuclear magnetic resonance, dynamic mechanical analysis, ultrasonic techniques and optical measurements (by infrared or Raman spectroscopy, grating-based optical fiber sensors (OFS) and Fresnel refractometers). Among the methods cited above, mainly DEA and OFS methods are able to monitor the cure during the manufacturing process.

Further studies have demonstrated the advantages of OFS over conventional sensors: they are small and flexible, light, immune to electromagnetic interference, and can withstand corrosive environments and high temperatures. Most of the published works on monitoring of epoxy cross-linking for composite processes by OFS use a specific sensing element [5] or the Fresnel reflection. The Fresnel reflection at the end of an optical fiber has been reported in Refs. 8, 10, 11. 15, and 17, while other Fresnel devices have been also developed [12,14]. Another kind of refractive index OFS is based on in-fiber grating structures, e.g., long period gratings [16] or tilted fiber Bragg gratings [16,18]. However, the sensitivity range is generally limited to refractive index below the refractive index of the optical fiber cladding. Optical-based sensors are in most cases calibrated by additional indirect measurements by DSC or TMA [6,10,11,15]. The drawback of this procedure, except the fact that both measurements are not done simultaneously, comes from the experimental conditions of DSC or TMA (size of the specimen, thermal environment) that are quite different from those encountered in more industrial situations. In this context, an interesting approach consists in using in the same experiment both heat flux sensors and OFS [17] for the determination of the degree of conversion.

Among the commercially available sensors in the context of industrial composite processes, the dielectric sensors allow in situ monitoring of the degree of cure. These sensors utilize the changes induced in the dielectric behavior of the resin to provide real-time and in situ information on the advancement of the process, generally also derived from some calibrations with the help of a quantitative analysis, for example, by DSC [9,13,19]. While advantages lie in its robustness, drawbacks come from the need of external calibration from DSC, and the size of the rigid sensors that provide only surface information.

Except methods based on heat flux measurements (DSC), all the devices described above are also sensitive to the temperature variations. Therefore, in cases of anisothermal curing, the calibration of the sensor needs to be temperature dependent. A supplementary appropriate evolution law is thus required for the measured quantity. For instance, in Ref. 20, the non-isothermal thermoset cure is monitored by an OFS, and the refractometer signal is expressed as the sum of two linear functions of temperature and chemical conversion.

Another way to determine the curing cycle of composites involves the modeling of the polymerization kinetics. For thermosetting resins, it is usually performed thanks to phenomenological models, requiring parameter identification from experimental data. In the case of epoxy/amine systems, Horie et al. [21] assumed an identical reactivity for primary and secondary amine functions leading to a model with a global reaction order of 3. A generalization of this model was proposed by Kamal and Sourour [22] including partial orders of catalytic and autocatalytic reactions to be identified, making the model very suitable for simulation of the curing of epoxy/amine systems. The main limitation stands for high degree of conversion when the reaction is mainly controlled by diffusion. This phenomenon can be taken into account by adding a diffusion factor, as proposed, for example, by Lam [23].

In this study we report on the cure monitoring of an epoxy resin suitable for liquid composite processes, by using simultaneously a multi-instrumentation based on an OFS (a Fresnel refractometer) and DEA. Two mold materials and sample thicknesses are used to enhance or reduce the exothermic phenomenon, providing isothermal and anisothermal curing data. The isothermal curing is exploited to correlate the optical signal, the dielectric signal, and the results of a kinetic modeling. It allows to base our analyses of the curing only on the signal from the weak-intrusive OFS. In the anisothermal case, an evolution law of the Fresnel power with temperature is proposed, based on the experimental results. Quantitative values of the degree of conversion are deduced from the experimental Fresnel signal by correcting it using the temperature changes. This approach is finally validated with a numerical model of the cure.

EXPERIMENT

Dielectric Analysis

Dielectric responses of thermoset resins are caused by the orientation of dipoles existing in functional groups and mainly (at low frequency or low temperature) by migration of ions present in epoxy resins due to their inherent impurities [24]. In a typical dielectric test, the resin sample is placed in contact with two electrodes of the dielectric sensor, and a sinusoidal voltage is applied to one electrode. The resulting sinusoidal current measured at the second electrode is attenuated in amplitude and shifted in phase in relation to the mobility of ions. One can show [24] that the ionic conductivity is proportional to the concentration of ions and to their mobility, which is a function of the (dynamic) resin viscosity. Therefore, the change in ion conductivity, related to the inverse of the ionic viscosity, can be used as an indicator of the curing state of resins and therefore composites. In this work, a 10-channel dielectric sensor system, DEA 230/10, from Netzsch Instruments was used. The sensors chosen for the experiments are single use sensors comprised of two interdigitated electrodes on a polyimide substrate. These sensors of size 10 3 20 mm 2 are simply placed in the mold in the vicinity of the resin. Dielectric measurements were conducted at various frequencies from 0.1 Hz to 10 kHz.

Fresnel Refractometer

An alternative approach makes use of the changes in the refractive index of a thermoset resin that occur during the curing due to resin shrinkage, by measuring the Fresnel reflection at the end face of an optical fiber [6,10,11,[15][16][17]. For a typical epoxy resin, refractive index shows an increase of 0.05 for a fully cross-linked resin [17]. The principle of Fresnel reflection is briefly recalled: a light coming out from an optical source is launched into an optical fiber through a coupler. The tip of the optical fiber is cleaved very precisely and immersed into the resin. At the interface between the fiber core and the resin, there is an abrupt change in the refractive index. Since the fiber end face is cleaved perpendicular to its axis, the reflection coefficient is given by Fresnel laws. A small part of light that arrives is backreflected at the interface and passes again through the coupler and is detected by a photoreceiver. Its intensity depends directly on the resin refractive index. This technique has been demonstrated to be relevant to monitor the evolution of the cure processing [7,8,10,15,17]. Drawbacks come from the facts that this technique needs to be externally calibrated to give quantitative results and is not totally reproducible as the reflection coefficient depends on the quality of the optical fiber end cleavage. Note finally that temperature produces similar changes in the signal and needs to be compensated in cases of an anisothermal curing that can be a difficult task. In this study, the Fresnel reflection setup consists of a very stable laser source (Nettest V R Tunics Plus CL) with an output power of 10 mW operated at 1550 nm, connected to an optical circulator, and a high resolution Rifocs 575L power meter. An optical fiber pigtail of SMF28 was used, and attention was paid for the fiber optic end cut. The end face of the fiber is positioned in the mold in contact with the resin and in the vicinity of the dielectric sensor in cases of cross-correlation of both measurements. As there is no need for knowing the absolute value of the resin refractive index, the Fresnel reflected power is simply recorded versus time and is enough to monitor the relative evolution of the cure.

Materials

The resin used in this work is the HexFlow V R RTM6 from Hexcel V R Composites, a premixed epoxy specifically dedicated to liquid composite processes, and commonly used in the manufacturing of composite reinforced by carbon fibers mainly in the aeronautic industry. In order to enhance or not the exothermal manifestation of the resin during its polymerization, two kinds of molds were used. A soft metallic (SM) mold, which was a 0.05 mm thick aluminum foil pan easily deformable, enabled good thermal transfer, containing 2.4 mm thick resin sample. A 1 mm thick silicon (SI) mold that exhibited low thermal conductivity was also used. In this case, the thickness of the resin sample was 12.8 mm. These molds were simply placed on a steel plate in a thermally controlled oven (V€ otsch VT4002). Type K micro-thermocouples were used with a Pico V R TC-08 controller device to assess the temperature in the resin and the oven during all curing cycles.

CURE MODELING

This section details the equations used to calculate the degree of conversion a of the resin in order to compare our results to those of the literature. The RTM6 resin curing kinetics has been widely studied [17,[25][START_REF] Lecointe | Caract erisation, simulation des processus de transferts lors d'injection de r esine pour le proc ed e RTM[END_REF][START_REF] Karkanas | [END_REF][28], and reliable models are available in the literature. Let us consider a thin plane infinite volume whose thickness is a. This volume is occupied by a thermosetting resin of density q, specific heat c P , thermal conductivity k, and conversion enthalpy DH T . The two faces exchange heat with the environment at temperature T env with a different heat exchange coefficient on the upper and the lower boundaries, respectively, h u and h l . The one-dimensional heat transfer equation given by Eq. 1 governs the resin temperature T subjected to a source term proportional to the conversion rate.

qc P oT ot 5r krT ð Þ1qDH T oa ot : (1)
The heat flux q C exchanged with the environment through the boundaries C gives

q Cu;l 5h u;l T2T env ð Þ on C u;l : (2)
If the resin sample is thermally thin, T can be assumed to be uniform in the thickness. The total heat energy balance in the volume delimited by an in-plane surface of area A simplifies to Eq. 3, by neglecting the heat conduction term:

qc P oT ot aA5qDH T oa ot aA12h à T2T env ð Þ A h à 5 h u 1h l 2 : (3)
T can easily be derived from Eq. 3 as a function of the conversion rate and environment temperature, by the following relationship:

oT ot 5 DH T c P oa ot 1 2h à a:q:c P T2T env ð Þ : (4)
In Eq. 4, DH T , c P , and q are the material parameters whose values are reported by Lecointe [START_REF] Lecointe | Caract erisation, simulation des processus de transferts lors d'injection de r esine pour le proc ed e RTM[END_REF]. DH T is set to 414 JÁg 21 , other authors report similar values of 432 JÁg 21 [17]. The specific heat capacity c P is assumed to be linearly dependent on the degree of conversion [29]. Thus, the step change around the vitrification/devitrification point is neglected. It has been verified by a parametric study that the results are not significantly affected by this assumption:

c P a; T ð Þ5c P a51; T ð Þ a1c P a50; T ð Þ12a ð Þ: (5)
The specific heat temperature dependence of the uncured and fully cured resins, respectively, c P a50; T ð Þ and c P a51; T ð Þ , is fitted with a second-order polynomial as given in Ref. 29. The linear relationship between the specific volume and a established by Aduriz et al. [17] allows to calculate the resin density as a function of a:

q a ð Þ5q a51 ð Þa1q a50 ð Þ 12a ð Þ (6)
using the density of the uncured and fully cured resin, respectively, q a50 ð Þ5 1117 kgÁm 23 and q a51 ð Þ5 1141 kgÁm 23 , reported in Ref. 29. In this approach, the thermal expansion is neglected. A modified Kamal-Sourour model [22] is used to calculate the conversion rate at each time step. A diffusion term has been added to the initial law in order to evaluate more accurately the slowed down conversion kinetic of the vitrified material [23]. The values of the coefficients of Eq. 7 used in the present work are taken from Refs. 25, 27, and 28 and given in Table 1. R is the ideal gas constant.

oa ot 5 k 1 k d k 1 1k d 1 k 2 k d k 2 1k d a m 12a ð Þ n k i 5A i e 2E i =RT ; i51; 2 k d 5A d e 2E d =RT e 2b=f f 5f 0 1d T2T g À Á : 8 > > > > < > > > > : (7)
In order to take into account a diffusion mechanism in the modified Kamal-Sourour law, the glass transition temperature T g of the resin must be evaluated as a function of a. This is achieved by using the so-called DiBenedetto equation [START_REF] Pascault | Thermosetting Polymer[END_REF], Eq. 8, where the glass transition temperature T g of the uncured and fully cured resin, respectively, T g,0 5 215 C and T g,1 5 220 C, as well as the value of the adjustable coefficient k 5 0.2784, is reported by Aduriz et al. [17].

T g 5T g;0 1 T g;1 2T g;0 À Á ka 12 12k ð Þa : (8)
A finite difference integration scheme is applied to Eq. 4, coupled to Eqs. 5-8, in order to evaluate the resin temperature evolution with time when the thickness of the resin sample a, the initial temperature T 0 , and the initial degree of conversion a 0 are known, as well as the temperature of the environment versus time, T env (t), measured by a thermocouple located in the oven.

The exothermal source term related to resin cfonversion enthalpy leads to a resin temperature arising at a higher value than the environment. The value of this exothermal peak is directly related to the heat transfer coefficient between the resin and its environment, h*. It is important to note that h* is actually an effective heat transfer coefficient that includes not only convection transfer effects, but also the contribution of the mold and the conductive transfer through the thickness of the resin. This parameter is the only unknown in the present model and depends on the experiment conditions, in particular on the thickness of the mold and the sample. It is proposed to identify the value of the coefficient h* by minimizing, using a Newton optimization algorithm, the gap between the measured and calculated maximum temperatures reached during the cure.

The value of h* has been identified for each experiment and allows to calculate the Biot number. It turns out that the assumption of uniform temperature field is only valid for the SM experiment (Biot number around 0.3). For the SI experiment, the large thickness of the sample involves a nonnegligible thermal gradient.

To summarize the material parameters being provided by the literature, the only data required to fully adjust the presented model is the resin maximal temperature. A validation of the applicability of the model to the experimental cases is obtained by comparing the calculated and measured resin temperatures versus time during the whole experiment. Once the validity of the simulation is established in "Isothermal cure monitoring assessment" section, the response of the model will be used as a reference in "Comparison to the model" section.

RESULTS AND DISCUSSION

In "Isothermal cure monitoring assessment" section, the responses of DEA and OFS are compared, allowing to analyze the curing comparatively to the results of the aforementioned model. In "Anisothermal cure monitoring assessment" section, it is experimentally demonstrated that a can be evaluated even in the case of anisothermal curing by correcting the Fresnel signal from its thermal variations. Finally, it is shown in "Comparison to the model" section that after an updating of the initial conditions, the model fits well the experimental results in anisothermal curing conditions.

Isothermal Cure Monitoring Assessment

Figure 1 shows the evolution of the ionic viscosity from the DEA, the Fresnel reflected power, and the thermocouple temperature measured in the vicinity of the other sensors during a thermal cycle. The cure is done at a relatively low temperature level (140 C) in the SM mold with a thin resin sample in order to limit the temperature overshoot due to the exothermal reaction. The thermal cycle consists in an increase of temperature from 80 C to 140 C at a heating rate of 6 CÁmin 21 to reduce the temperature gradients (heating phase), followed by a cure plateau where the resin temperature remains almost constant. The third period is a simple cooldown by natural convection to ambient temperature (cooling phase).

We show clearly that the correlation between the dielectric and the optical signals is very good, for at least the two first phases. During the heating phase, the viscosity decreases until a minimum corresponding to the onset of cross-linking (point A in Fig. 1). The optical signal shows also a decrease since the refractive index is only influenced by the thermal changes. During the cure plateau, the cross-linking occurs (increase of the degree of conversion), involving an exothermal reaction which is dissipated by the metallic mold. These phenomena lead to an increase of the molar weight of the resin macromolecules and correspond physically to a densification of the resin network and shrinkage. It is manifested by an increase of the ionic viscosity in Fig. 1a. The reflected power follows the same evolution with time, and as an isothermal state can be considered, the increase in the reflected power is directly linked to the increase in the refractive index of the resin. It is noteworthy that the refractive index of the resin is proportional to the density and to the molecular weight of polymer repeat units following the Lorentz-Lorenz law [11]. It was also shown in Ref. 17 that there is a linear correlation between the resin refractive index and the degree of conversion a. Furthermore, the relationship between the refractive index and the Fresnel reflected power can be considered, in a first approximation, as linear [11], because the amplitude of the refractive index variation is of 10 22 order that makes possible the linearization. This has been verified, as refractive index values ranging between 1.54 and 1.58 are reported in Ref. 17 for the RTM6 during its curing. It is thus concluded that the Fresnel power is a linear function of a at constant temperature, as also the ionic viscosity. After an inflection point of the DEA and Fresnel curves (point B in Fig. 1a), the attainment of the plateau of both signals correlates with the end of the cross-linking (point C in Fig. 1a). For the third phase (cooling), the dielectric and the optical signals differ because the ionic conductivity is slightly dependent on the temperature, and the sensor can be saturated by cured resin. On the contrary, the reflected power signal shows a monotonic function of the temperature as the polymerization reaction is finished.

In order to assess the correlation between the electric and optical signals, both are plotted altogether as depicted in Fig. 2. This figure shows that the better correlation between the two signals is get for the second phase (isothermal) and an almost linear correlation is found, demonstrating experimentally that the Fresnel refractometer can be used instead of the dielectric sensor, with all the advantages of OFS. In the first phase (heating), the correlation is also good whereas not so linear. Nevertheless, the minimum value of the ionic viscosity corresponds to the minimum value of the reflected power. Note that small differences between both signals can be attributed to the fact that the sensors are not exactly situated at the same location, and their spatial resolution is different. Finally, the dielectric sensor is unable to detect correctly the thermal variation of the cured resin during the cooling, and no correlation is seen in Fig. 2

(ellipse).

During the whole curing cycle, an analysis of the Fresnel reflected power can be done, as illustrated in Fig. 3. In the fist part of the curve, the resin moves from viscous to fluid, and the plot in Fig. 3 shows a decrease of the refractometer signal with the temperature, corresponding to a decrease of the refractive index with the temperature. During the cure phase, an increase of the optical signal from 16.6 mW to 20 mW is measured, whereas the temperature remains almost constant. In the cooling phase, as suggested by Giordano et al. [START_REF] Giordano | [END_REF] the glass transition temperature (T g ) can be evaluated by the inflection point in Fig. 3 during the cooling, around 110 C, which is related to the change in the thermo-optic effect by crossing the glass transition zone.

As explained previously, the Fresnel reflection power is a linear function of a but cannot provide directly quantified values without calibration. A quantitative comparison of this experimental data to the degree of conversion is however possible by using the thermo-kinetic model presented in "Cure modeling" section. The model is adjusted to the actual experimental conditions by identifying the effective heat transfer coefficient only using the maximum temperature reached by the resin (139.8 C provided by the thermocouple). After convergence, a reasonable value of h* 5 65 WÁm 22 ÁK 21 is found. To analyze the results, we first focus on the temperature evolution versus time. Figure 4 shows that the results of the model fit almost perfectly to the whole experimental curve, ensuring the relevance of the model. A maximum deviation of 1.4 C is found in the temperature increase, but the most significant bias is located just after the weak exothermal manifestation (bias lower than 0.3 C in cure plateau), as can be seen in the inserted plot in Fig. 4.

The model also provides the evolution of the conversion degree a which is plotted versus the Fresnel reflected power in Fig. 5. As expected, in the almost isothermal condition involved in this case, a first rough analysis shows that the Fresnel reflected power is a linear function of a. In detail, before the onset of polymerization, the decrease of the Fresnel reflected power is directly linked to the thermal sensitivity of the signal. Moreover, small deviation to the linearity at the end of the cure can be associated to an overestimation of the calculated temperature after the exothermic temperature peak, from about 120 to 180 min. It leads to an overestimation of a given by the model. The model also provides an estimation of the evolution of T g during the curing process. The final value at the end of the cure reaches 127 C which is necessarily also overestimated. This is however consistent with the glass transition manifestation deduced from the Fresnel power curve in Fig. 3, which exhibits a slope change around 110 C.

Anisothermal Cure Monitoring Assessment

Another curing experience has been conducted in the insulating SI mold with a substantial thick resin sample, allowing large exothermal manifestations and an anisothermal cure. Note that before this experiment, the RTM6 sample has experienced an inadvertent thermal exposure suspected to have slightly modified the initial degree of cure. The thermal cycle consists in an increase of temperature from 80 C to 150 C at a heating rate of 6 CÁmin 21 , followed by a cure plateau where the temperature setting stands at 150 C, and then by a cooldown by natural convection to room temperature. Figure 6 shows the evolution of the Fresnel reflected power as a function of the thermocouple temperature during the whole cycle. During the heating step, the Fresnel power reflected in the uncured resin (a 5 0) decreases linearly with the temperature up to 140 C when the exothermal reaction begins to occur. During the exothermal polymerization (the resin temperature is above 150 C during about 50 min), both temperature and a evolve. This results in a soft increase of the Fresnel power. Contrarily to the curing in the SM mold, as expected it was observed a huge overshoot of the temperature which corresponds to the exothermal manifestation of the chemical bonding: the temperature reaches more than 240 C. It is explained because the diffusion of the temperature was not at all facilitated in the SI mold, and only low heat flux was able to be exchanged between the sample and the mold. The large thermal exposure undergone by the resin allows to assume that the polymerization was fully completed in the linear part of the temperature decrease (below 210 C). This is consistent with the fact that almost no evolution of the Fresnel power was observed when the resin temperature was maintained at the temperature setting during about 1 h (see ellipse in Fig. 6).

It was shown in "Isothermal cure monitoring assessment" section that the Fresnel power P is an affine function of a. The Fresnel power also evolves linearly with the temperature so that P can be written as:

P a; T ð Þ5P a51; T ð Þ a1P a50; T ð Þ12a ð Þ: (9)
Both terms P a50; T ð Þ and P a51; T ð Þ are deduced from the experimental data presented in Fig. 6 by fitting an affine function represented, respectively, by the dashed and dot-dashed lines in this figure. By inverting Eq. 9, the temporal evolution of a can be easily deduced from the experimental curve P a; T ð Þ in Fig. 6. The results are plotted in Fig. 7 together with the measured temperature. It is clearly shown that the plot of a is consistent with the temperature evolution. The measured degree of conversion oscillates around zero up to a temperature about 125 C. At this time (15 min), the polymerization starts and the conversion rate accelerates. After the exothermal peak, the temperature decreases while a still increases more slowly until it reaches its maximal value, 1. The signal then oscillates around this last value during the cooling phase due to noise in the Fresnel signal (see also Fig. 6) without any physical meaning.

Comparison to the Model

It was shown in the previous section that the Fresnel reflected power provides an evaluation of a even under anisothermal conditions. In order to assess the relevance of a measured by this way, a comparison to a model is also proposed. However, the model described in "Cure modeling" section is not relevant because of the thermal gradients. For that purpose, a two-dimensional (2D) finite element numerical simulation has been achieved in order to take into account the temperature gradient through the thickness of the sample. The model includes the 1 mm thick SI mold (conductivity 0.23 WÁm 21 ÁK 21 , density 1200 kgÁm 23 , and specific heat 1590 JÁkg 21 ÁK 21 ) and the 12.8 mm thick RTM6 sample (conductivity ranging from 0.1 to 0.23 WÁm 21 ÁK 21 and expressed as a function of a, according to Ref. [START_REF] Lecointe | Caract erisation, simulation des processus de transferts lors d'injection de r esine pour le proc ed e RTM[END_REF]). The model solves the 2D heat transfer equation Eq. 1 coupled to Eqs. 5-8. Boundary conditions are temperature and flux continuity at the mold/sample interface and convective heat transfer between the external boundaries and the environment at T env controlled by the heat transfer coefficient h. Figure 8a summarizes the finite element modeling and the meshes. Figure 8b andc presents, respectively, the temperature field and the degree of conversion through the thickness at the exothermal peak for which the gradient is maximal. The temperature deviation exceeds more than 30 C, and the degree of conversion ranges from 71% to 94%, requiring to consider accurately the actual sensor location. Figure 8b also locates the thermocouple and the Fresnel sensor in the thickness of the resin sample. The values extracted at points A and B in Fig. 8b (respectively, C andD) are used to compare FEA results to the experimental temperature (respectively, Fresnel signal) in order to take into account the sensor location uncertainties.

As previously conducted in "Isothermal cure monitoring assessment" section, the heat exchange coefficient h was identified by adjusting the maximum experimental temperature of the resin (i.e., at the exothermal peak) to the simulated one at points A and B. The optimized value found for h is 45 WÁm 22 ÁK 21 , a relevant order of magnitude for natural convection in air. In Fig. 9a are plotted the FEA results (degree of conversion) at both points C and D and in Fig. 9b the FEA results (temperature) at both points A and B. Despite large gradients observed in the SI curing experiment, small discrepancies are finally noticed in the vicinity of both sensors. Figure 9a also compares the experimental and calculated conversion kinetics with the assumption that the initial degree of conversion is zero. It is shown that both curves agree qualitatively. However, it turns out that the model is delayed of about 14 min. This is obviously related to the delay observed also between the experimental and calculated temperatures (see Fig. 9b). In our opinion, this discrepancy is necessarily related to a nonzero initial degree of conversion of the RTM6 sample at the beginning of this experiment, as mentioned previously. For that reason, the initial value of the degree of conversion a 0 was also identified to adjust the experimental and calculated time-temperature data. The last results are plotted in Fig. 10. As shown in Fig. 10b, a very proper agreement is obtained with a reasonable value of the initial degree of conversion a 0 5 6.25%, explained by a long exposure to the initial temperature before to conduct the SI experience. Naturally, the Fresnel data need to be processed also using this value of a 0 . This is achieved using the following equation:

P a; T ð Þ5P a51; T ð Þ a 0 2a a 0 21 1P a5a 0 ; T ð Þ 12a 12a 0 ; (10)
where P a5a 0 ; T ð Þis fitted from Fig. 6 using the linear part of the curve before the onset of polymerization. It is shown in Fig. 10a that a good quantitative agreement is obtained between the measured and calculated a using these initial conditions. Both temperature and degree of conversion evolution are very properly described by the model.

It was thus demonstrated that an experimental evaluation of the degree of conversion of a thermosetting resin could be obtained using a Fresnel refractometer even under anisothermal conditions. An absolute value of the degree of conversion is available if the measurement is calibrated. An in situ calibration was performed using the evolution with temperature of the measured Fresnel power for two different constant degrees of conversion. The initial degree of conversion is constant up to the onset of polymerization. A first calibration set of data is then available at the beginning of the experiment if the initial degree of conversion is known. The second set of data available for the calibration of the measurement is obtained at the end of the experiment, when the resin is fully polymerized. The cooldown step then provides the variation of the Fresnel power as a function of temperature for a constant a. If the initial degree of conversion is known, this method provides a posteriori an absolute value of a once the polymerization has ended. It could also be used in real time to obtain a relative value of a and then monitor the process. To enhance the method and obtain an absolute value of a, it is mandatory to calibrate the setup a priori by monitoring a cure of the resin during a cycle with high temperature, as performed here in the SI mold. Nevertheless, it is not sure that the cleaning of the fiber end after calibration does not affect the cleaved surface and thus the refraction of the light. The second way to provide a priori an absolute value is to ensure a perfect cleavage of the fiber end allowing to determine the refraction index of the resin from the Fresnel power. If a reproducible evolution of the refraction index with temperature and degree of conversion could be available, an absolute value of the degree of conversion would be provided directly from the Fresnel power measurement.

SUMMARY AND CONCLUSIONS

In this paper, the cure of the RTM6 thermosetting resin was investigated using a thermocouple and an optical fiber Fresnel refractometer. The quasi-perfect correlation between the optical signal and the signal from a commercial DEA intended to the composite manufacturing industry was first experimentally demonstrated. It allowed to base our analyses only on the Fresnel power and temperature evolutions. Two kinds of mold materials and sample thicknesses were used to enhance or reduce the exothermic phenomenon. A cure simulation was also performed in order to compare it to our experimental results. Two assumptions have been made: for thin samples, temperature gradient is neglected and a semi-analytical model is relevant, whereas gradients have been taken into account for thick samples thanks to a finite element analysis. For the isothermal cure, it was shown that the Fresnel signal evolves linearly with the calculated degree of conversion, and a perfect agreement was achieved for the whole temperature curves. In the case of anisothermal cure caused by a confined exothermic reaction, the degree of conversion was experimentally found by correcting the Fresnel power from the temperature variations, provided that the initial and final degrees of conversion are known. A comparison to the FEA model showed a mismatch with the experimental results, which has been cancelled out by adjusting the initial degree of conversion. An absolute value of the degree of conversion is thus provided a posteriori once the polymerization has ended. This approach can be directly applied in real time to obtain a relative value of the degree of conversion of a composite part instrumented by Fresnel and temperature sensors, during the manufacturing process.

FIG. 2 .

 2 FIG. 2. Fresnel reflected power as a function of the ionic viscosity during the whole curing process (SM). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FIG. 4 .

 4 FIG. 4. Measured and calculated resin temperatures as functions of time during the curing process in the SM mold. A detailed view in the 134-140 C temperature range is inserted. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FIG. 7 .

 7 FIG. 7. Experimental degree of conversion and temperature evolution during the curing process (SI). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FIG. 8 .

 8 FIG. 8. Finite element analysis. (a) FEA model, meshes, and boundary conditions. Results at the exothermal peak: (b) temperature and (c) degree of conversion through the thickness of the sample. The locations of the sensors are also presented in this figure. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

TABLE 1 .

 1 Values of the coefficients of the Kamal-Sourour model (Eq. 7).

	Coef. (unit)	A 1 (min 21 )	E 1 (JÁmol 21 )	A 2 (min 21 )	E 2 (JÁmol 21 )	A d (min 21 )	E d (JÁmol 21 )
		4.50 3 10 6	7.47 3 10 4	1.30 3 10 6	5.84 3 10 4	1.42 3 10 28	2.16 3 10 5
	Coef. (unit)	d (-)	f 0 (-)	b (-)	m (-)	n (-)
		4.80 3 10 4	0.025	0.227	1.216	1.284

FIG. 1. Ionic viscosity and Fresnel reflected power (a) and temperature (b) evolution recorded for the whole temperature cycle applied to the resin sample (SM mold). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

ACKNOWLEDGMENTS

The authors are grateful for the financial support by the French Ministry of the Economy, Finance and Industry in the framework of the Carnot M.I.N.E.S. Institute and acknowledge C. Vergnes for his great contribution to the experimental work.