Jan Rojek
email: jan.rojek@mines-paristech.fr

Sébastien Joannès
email: sebastien.joannes@mines-paristech.fr

Jean-Christophe Teissèdre
email: jean-christophe.teissedre@mines-paristech.fr

Lucien Laiarinandrasana
email: lucien.laiarinandrasana@mines-paristech.fr

Anthony R Bunsell
email: anthony.bunsell@mines-paristech.fr

TWENTY-SECOND INTERNATIONAL CONFERENCE ON COMPOSITE MATERIALS (ICCM22) EFFECT OF THROUGH-THICKNESS COMPRESSIVE STRESS AND POROSITY ON THE TENSILE STRENGTH OF CARBON-FIBRE REINFORCED COMPOSITES

Keywords: biaxial, porosity, CFRP

Low and high-porosity cross-ply CFRP specimens are manufactured by modifying the autoclave curing cycle. A custom compression setup is designed and integrated with an electromechanical tensile rig to allow biaxial testing of the material. Tensile strength of the specimens is tested at different levels of applied through-thickness compression. A reduction of tensile strength due to through-thickness compression is observed. This effect is magnified by the presence of porosity. The ply level stress state is analysed through a three-dimensional finite element contact simulation.

INTRODUCTION

Understanding the mechanical behaviour of fibre-reinforced composites under multiaxial loads is important for their applications in thick-walled and geometrically complex parts. Existing experimental results suggest that through-thickness compression can significantly reduce the tensile strength in the fibre direction [START_REF] Gan | Effect of high through-thickness compressive stress on fibre direction tensile strength of carbon/epoxy composite laminates[END_REF][START_REF] Goto | Strength evaluation of unidirectional carbon fiber-reinforced plastic laminates based on tension-compression biaxial stress tests[END_REF].

Voids are perhaps the most widely studied type of defects in composites [START_REF] Mehdikhani | Voids in fiber-reinforced polymer composites: A review on their formation, characteristics, and effects on mechanical performance[END_REF]. Nonetheless, little is known about their influence on composite strength under multiaxial loading conditions. This type of loading could be potentially relevant in case of thick-walled filament-wound pressure vessels, which are known to have a significantly more irregular microstructure than that obtained through other methods, e.g. prepregging. In particular, mesoscale voids elongated along fibre direction are commonly observed [START_REF] Rojek | Modelling the effect of porosity on the mechanical properties of unidirectional composites. The case of thick-walled pressure vessels[END_REF].

Figure 1: CT view of a cross-section of a highly porous CFRP laminate.

The study described here aims at evaluating the effect of multiaxial stresses on the tensile strength of highly porous carbon fibre-reinforced composites. The response of a carbon-fibre reinforced polymer to a combination of longitudinal tension and through-thickness compression is studied. In addition to a reference low-porosity material, highly porous laminate specimens were manufactured and tested to assess the influence of voids on the material strength under these conditions.

MATERIAL

Carbon/epoxy laminate plates of dimensions 30 cm × 30 cm and a stacking sequence [0/90 2 /0] 2S were manufactured by manual lay-up using an airspace-grade preimpregnated unidirectional tape HexPly 6376C-HTS(12K)-5-35%, with a longitudinal tensile strength of approximately 2500 MPa. After being placed in the stack, each ply was manually rolled to remove excess air. Next, the plates were autoclave-cured. Reference plates (from now on referred to as non-porous) were made following the manufacturer-recommended cycle. Temperature was ramped up to 175 °C over 65 min. Then, it was held constant for 150 min, before starting the cool-down phase. Throughout the process, a vacuum and a 7 bar overpressure was maintained. Highly porous plates were manufactured by modifying the cycle explained above. The temperature was the same as for the reference case. Vacuum was maintained throughout the cycle. However, no overpressure was applied. Low curing pressure is known to lead to formation of voids in fibre-reinforced composites and a similar approach has been used by other authors to obtain a highly porous material [START_REF] Olivier | Effects of cure cycle pressure and voids on some mechanical properties of carbon/epoxy laminates[END_REF][START_REF] Kakakasery | Cure cycle effect on impact resistance under elevated temperatures in carbon prepreg laminates investigated using acoustic emission[END_REF].

Figure 1 shows a detail of a porous specimen observed using micro-computed tomography. The voids are strongly elongated along the fibre direction, and their cross-sectional dimensions are typically significantly larger than fibre diameter. The voids in the manufactured plates were not uniformly distributed. The highest porosity was usually present near the mid-thickness line. Typically, the porosity fraction observed through optical microscopy was in the order of 5%, although it varied significantly between specimens.

The laminate plates were cut down into rectangular specimens of dimensions 250 mm x 10 mm. The average thickness of the specimens was equal to 2.15 mm for the non-porous and 2.20 mm for the porous ones. Glass/epoxy end tabs of 60 mm length were used to reduce stress concentrations at the tensile machine grips.

METHODS

Tensile testing

The specimens described in the previous section were tested in tension on an electromechanical rig, at a constant cross-head speed of 1mm/min. Selected specimens had their edge polished to allow insitu observation of damage development at the micro and mesoscale. These specimens were loaded in tension with several intermediate hold phases. During those, an approximately 30mm long section of the edge was mapped using a motorized digital microscope, allowing observation of damage development. For a single porous specimen, the microscope was focused on an approximately 0.75 mm × 0.75 mm region of the specimen edge, where the damage was then recorded in situ without interrupting the test, allowing a more detailed local observation.

Biaxial testing

A custom hydraulic compression system was developed to allow biaxial testing of the specimens. The requirements at the design stage were that the system be compatible with a standard tensile rig and portable. A motorcycle brake caliper with a load range up to 100 bar was used as a base of the system. A custom-designed indenter support is presented in Figure 2. The load is transferred to the specimens using cylindrical indenters with a 10 mm radius. The alignment of the indenters is controlled through the use of bearing-supported guiding rods.

The compression system described above was mounted on a tensile machine, with four counterweights on pulley blocks used to avoid passing additional loads on the specimens, as shown in Figure 3. To carry out biaxial tests, the compressive load was first applied. The tested pressure levels were 25, 50 and 100 bar. Once the desired pressure was reached, it was maintained using an automatic control unit. The specimen was then loaded in tension until failure at a cross-head speed of 1 mm/min.

Additional measurements

Several supplementary measurements were taken. A video acquisition setup was used during the biaxial tests to photograph the edge of the specimen in the compressed section every 2 seconds. The indenter displacement was automatically evaluated using in-house image analysis software, based on cross-correlation analysis. Furthermore, a speckle pattern was applied on selected specimens, and a DIC analysis was carried out using VIC-2D software. Contact area was measured at different pressure levels, using the method described by Gan et al. [START_REF] Gan | Effect of high through-thickness compressive stress on fibre direction tensile strength of carbon/epoxy composite laminates[END_REF]. The contact area was measured at 25, 50, 75 and 100 bar for the non-porous specimens. For the porous ones, it was measured at 20, 40, 60 and 80 bar.

RESULTS

Tensile testingmechanical properties

Mean tensile strength, strain and modulus are summarized in Table 2. The failure strain was very similar for the two studied materials. However, the highly porous specimens had a lower strength and modulus, by 7% and 6% respectively.

Low porosity High porosity

Strength Voids were seen to act as preferential crack origin sites in the transverse plies. The first observed cracks always appeared at the voids, at approximately 30% of the ultimate load. Only later did cracks appear elsewhere. The initial cracks at voids propagated over longer periods of time, while the cracks not passing through the voids appeared instantaneously. In the non-porous specimens, cracks started to appear later, at approximately 60% of the ultimate load and their density at failure was significantly lower. Figure 4 shows the development of transverse cracking and fibre breaks in a porous laminate subjected to tensile loading (in the vertical direction). The first image shows the pristine material, with several voids visible. The second image, taken shortly prior to failure shows several transverse cracks that formed in the cross-ply as well as numerous fibre breaks in the longitudinal ply.

Biaxal testing

The porous specimens all broke before reaching the maximum compressive load, on average at 84.8 bar. Figure 5 shows the average values of tensile strength at the different levels of applied through-thickness load. For the non-porous specimens, no decrease in tensile strength was observed at 25 bar of compression. For higher pressures, the average tensile strength decreased abruptly, being only 144 MPa at 100 bar. For the porous specimens, the decrease was even stronger. While in tensile tests they did not differ much from the non-porous ones, already at 25 bar they showed an approximately 20% lower tensile strength. This effect was even more pronounced at 50 bar, where the porous specimens were approximately 60% weaker. The porous specimen results were also characterized by a higher scatter.

FINITE ELEMENT ANALYSIS

Finite element simulations were carried out to compute the stress state in the low-porosity specimens subjected to biaxial loading. The simulations were carried out in an in-house finite element code Z-Set [START_REF] Besson | Large scale object-oriented finite element code design[END_REF]. The composite was modelled using a transversely isotropic linear elastic material behaviour. Only a quarter of the specimen was meshed, taking advantage of the material symmetries. The indenter-to-specimen contact was modelled using a frictionless formulation. Displacement boundary conditions were applied at the top of the indenter (compression) and at the far end of the specimen (tension). A structured mesh of 8-node brick elements (C3D8) was used, with the mesh being most refined under the indenters.

Figure 6 shows the simulated stresses along a through-thickness line, joining the tips of the indenters at mid-width of the specimen. The applied pressure for the case was 100 bar, and the tensile load of 144 MPa was chosen based on the average value at which the non-porous specimens failed at this pressure level. The through-thickness compression leads to the appearance of axial tensile stress in the longitudinal plies. This stress is further increased after the tensile load is applied. However, even at failure it is still significantly lower than the uniaxial tensile strength.

CONCLUSION

The difference in tensile strength between the non-porous and porous specimens is limited and could be in part attributed to the lower effective volume fraction of porous specimens. This is further confirmed by an approximately proportional decrease in the modulus value. The transverse cracking process is much more intensive in the high-porosity specimens, with voids leading to local stress concentrations, accelerating crack formation.

The biaxial test results suggest that a uniaxial longitudinal strength criterion is not sufficient in the presence of important transverse compressive loads. For a through-thickness compression of 100 bar, the longitudinal plies can fail at a longitudinal tensile load of approximately 1000 MPa, more than two times lower than in the uniaxial case. While a mesoscale failure criterion based on a combination of the three principal stresses could be proposed, it is not clear what is the underlying microscale mechanism leading to the reduced tensile strength. The presented results demonstrate that mesoscale porosity can further strengthen the negative effect of through-thickness compression on the longitudinal tensile strength. This could be caused by local transverse stress concentrations caused by the presence of voids. In this case, the material response would strongly depend on the local void morphology (shape, size, orientation). This would also provide an explanation for the higher experimental scatter observed in porous specimens.

The approach presented allowed assessing the influence of a multiaxial stress state on the strength of a highly porous composite material. A custom compression setup was designed and manufactured, allowing biaxial testing on a standard tensile rig. Three-dimensional finite element analysis was carried out to simulate the stresses inside the specimens. The obtained results suggest that a compressive through-thickness stress negatively affects the longitudinal tensile strength of carbonfibre reinforced composites. High porosity, although relatively benign under uniaxial longitudinal tension, was found to exacerbate the effect of through-thickness compression.

Figure 2 :

 2 Figure 2: Custom compression setup components. 1) caliper body, 2) cylindrical indenter, 3) load cell, 4) guiding rod.

Figure 3 :

 3 Figure 3: The compression system mounted on a tensile machine. The inclination of the caliper allows direct camera observation of the compressed specimen section.

Figure 4 :

 4 Figure 4: Damage development in a highly porous specimen under tension (applied vertically) observed in situ using a digital microscope. Pristine state (left) and shortly prior to failure (right).

Figure 5 :

 5 Figure 5: Average values of measured tensile strength as a function of applied compression.

Figure 6 :

 6 Figure 6: FEM simulation: principal stresses across thickness of a low-porosity specimen at 100 bar compression and 144 MPa applied tension. σ 11 is oriented along the longitudinal direction, σ 22 through the thickness of the specimen, and σ 33 along the width.

Table 1 :

 1

		1044 MPa	972 MPa
	Strain	1.49%	1.47%
	Modulus 70318 MPa	66163 MPa

Tensile test results for low and high-porosity CFRP.

Jan Rojek, Sébastien Joannès, Jean-Christophe Teissèdre, Lucien Laiarinandrasana and Anthony R. Bunsell

ACKNOWLEDGEMENTS

The research leading to these results has been carried out within the framework of the FiBreMoD project and has received funding from the European Union's Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement No 722626. The authors express their gratitude for the assistance of the professors and staff of the Composite Materials Group at KU Leuven for providing the materials and access to their facilities, where the specimens used in this study were manufactured.


Patrice Loisel

Guillerme Duvillié

Denis Barbeau

Brigitte Charnomordic

A user-friendly software to evaluate forestry scenarii including natural risk

Keywords: forest management, simulation, storm risk, natural risk, optimization, Decision Support System, stand level, Faustmann

published or not. The documents may come

• Aims We propose a user-friendly software to simulate even-aged and monospecific forest at the stand level, in order to evaluate and optimize forest management. The software gives the possibility to run management scenarii with or without considering the impact of natural risk. The control variables are the dates and rates of thinning and the cutting age.

• Methods The risk model is based on a Poisson processus. The Faustmann approach, including tree damage risk, is used to evaluate future benefits, economic or ecosystem services. It relies on the calculation of expected values, for which a dedicated mathematical development has been done. The optimized criteria used to evaluate the various scenarii are the Faustmann value and the Averaged yield value.

• Results We illustrate the approach and the software on two case studies: economic optimization of a beech stand and carbon sequestration optimization of a pine stand.

• Conclusion Software interface makes it easy for users to write their own (growth-tree damage-economic) models without advanced programming skills.

Introduction

Forest managers have always been interested in the evaluation of silviculture practices.

Nowadays, they must consider criteria that cover the increase in natural risks [START_REF] Hanewinkel | Climate change may cause severe loss in the economic value of european forest land[END_REF] due to climate change. Furthermore, in the face of growing societal demand, forest managers must also consider biodiversity, carbon sequestration and more generally ecosystemic services. For all these criteria, risk must be taken into account and silviculture should evolve accordingly.

Focusing on forest management at the stand level, a lot of growth models are now available, which were historically dedicated to one or two given species [START_REF] Monserud | A basal area increment model for individual trees growing in even-and uneven-aged forest stands in austria[END_REF][START_REF] Coates | Use of a spatially explicit individual-tree model (sortie/bc) to explore the implications of patchiness in structurally complex forests[END_REF][START_REF] Lacerte | Calibration of the forest vegetation simulator (fvs) model for the main forest species of ontario, canada[END_REF][START_REF] Pretzsch | The individual-tree-based stand simulator silva[END_REF]. Most of them are implemented in software: we do not give details here, but an extensive review can be found in [START_REF] Dufour-Kowalski | Capsis: an open software framework and community for forest growth modelling[END_REF]. Recent efforts have been made to integrate growth models in software platforms in order to fulfill the growing user demand. However existing software suffers from two main limitations : i) it requires a lot of efforts for users to customize the models used for the simulations ii) it does not integrate natural risks (except ForestGALES for wind risk [START_REF] Gardiner | Management of forests to reduce the risk of abiotic damage -a review with particular reference to the effects of strong winds[END_REF], but ForestGALES is very specific and must be coupled with forest growth software).

The motivation behind the present work is to answer the limitations of forest simulation software cited above. The challenge is to easily evaluate scenarii while eventually taking the natural risks into account. We propose a software, called EvaSylv which constitutes a significant advance regarding both these limitations. This software will facilitate the emerging of new silviculture guidelines.

EvaSylv consists of a chain of interconnected models, including a tree growth model, a risk model, a tree damage model following risk, and several economic criteria. The chain of models works at stand level, and allows to compute analytical expressions for criteria (Faustmann Value, Averaged Yield Value) linked to a given technical itinerary. This analytical approach has the advantage to yield precise results in a short computing time, compared to a Markov-based approach for computing expected values, which requires a great number of time-consuming simulations.

In EvaSylv, users can easily write their own models through a user-friendly interface.

The functions and models are considered as a supplementary kind of data, and can be modified via the interface. Contrary to other forest simulation software, users do not have to get into the complex software code. They are written in basic Python and do not require advanced programming skills.

The paper is organized as follows: the framework is introduced in Section 2, including a recall of the mathematical models used for simulation and optimization. Software design is presented in Section 3. Case studies relative to simulation and optimization of beech and pine stands are given in Section 4. Section 5 gives some conclusions and perspectives.

Framework

The objective of the present work is to provide an operational easy-to-use simulation and optimization framework for testing forest management scenarii at stand level.

EvaSylv has been designed to allow users to easily write their own models through a user-friendly interface: growth functions, damage models (linked to natural risk), economic valorizations. The functions and models are considered as a supplementary kind of data, and can be modified via the interface, without having to get into the complex software code. They are written in basic Python and do not require advanced programming skills.

The evaluation requires the choice of a technical itinerary describing the authorized actions. For each itinerary, the Faustmann approach is available to evaluate future benefits, economic or not. It necessitates the calculation of their expected values, for which a dedicated mathematical development was written [START_REF] Loisel | Faustmann rotation and population dynamics in the presence of a risk of destructive events[END_REF][START_REF] Loisel | Impact of storm risk on faustmann rotation[END_REF]. software. The elements shaded in gray correspond to the in-depth program structure and the elements within orange-colored containers correspond to the interface. Particular attention has been paid to the choice of the latter in order to free the user from writing complex mathematical calculations whilst ensuring a broad range of possibilities for designing the scenarii.

Software core

The software core consists of the definition of technical itineraries (with their corresponding scenarii), the set of mathematical models and evaluated criteria, used in the simulation or optimization module.

Different technical itineraries are available, each of them corresponding to a management strategy, taking or not taking the natural risks into account. Within the technical itinerary, the scenario parameters provide the settings for the initial forest state and for the thinnings: tree-density, dates and rates, as well as the choice of the criterion to optimize, if optimization is required. The software can be used in three main modes:

• Unconstrained simulation: apply a technical itinerary,

• Constrained simulation : apply an admissible (which satisfies constraints for optimization) technical itinerary,

• Optimization : find the optimized settings according to a given criterion.

The simulation is run over the rotation time T . It is based on four interconnected models: tree growth model, natural (such as storm) risk model, tree damage model and gain/cost model. The tree growth model output feeds both the tree damage model and the valorization one. The evaluation criteria, based on the Faustmann approach [START_REF] Faustmann | Berechnung des wertes welchen waldboden sowie noch nicht haubare holzbestände für die waldwirtschaft besitzen[END_REF], are calculated differently whether or not natural risk is taken into account in the technical itinerary. In the no risk case, the criteria are deterministic [START_REF] Faustmann | Berechnung des wertes welchen waldboden sowie noch nicht haubare holzbestände für die waldwirtschaft besitzen[END_REF], in case of risk [START_REF] Reed | The effects of the risk of fire on the optimal rotation of a forest[END_REF], they are based on expected values or CVaR.

Technical itineraries

The technical itineraries constitute a simple framework, based on production rules which its own consequences in terms of income.

In the present version, two technical itineraries are available, the first one for the no risk case and the second one for the case of risk. In the no risk case, stand management follows a basic non ambiguous itinerary: a clearcutting and a regeneration of the stand are done at the end of rotation period (or cutting age) T . In presence of natural risk, different technical itineraries may be used by forest managers. The software implements an adaptation to storm risk [START_REF] Loisel | Impact of storm risk on faustmann rotation[END_REF] of the itinerary described in [START_REF] Reed | The effects of the risk of fire on the optimal rotation of a forest[END_REF], originally designed for forest fire risk and widely used as a reference in presence of natural risk. It handles storm risk by defining a rule base, based on the fact that only when trees reach a certain height, storms cause direct damage but also weaken the trees, increasing future damage. The rule base is as follows:

• If a storm occurs before Threshold time (user-defined)t L : no or neglectable impact on trees, then forestry management is unchanged.

• If a storm occurs at time τ between t L and T : the ratio of damaged trees is high, then a clearcutting and a regeneration of the stand are performed at time τ .

• if no storm occurs before time T : a clearcutting and a regeneration of the stand are done at time T .

In the work presented in this paper, the considered natural risks are storms, though our approach remains valid for other natural risks, provided that suited itineraries are elaborated. To each technical itinerary corresponds associated criteria. The criteria are of different kinds, deterministic or expected actualized value of future incomes (Faustmann Value) or expected income (Yield Value), For the expected cases, the expressions of the expectations are hard-coded in the software. In order to evaluate the different criteria, it is necessary to evaluate incomes and costs. To this issue, various models must be given.

A sequence of models

The various model parameters and functions (see Table 1) are of four different kinds, detailed below:

• tree growth model (at stand level). It describes the dynamics of the tree basal area s: increment growth rate and the function describing the evolution of the dominant H and average Hmoy tree heights. It is a continuous model which is then discretized.

• natural risk model. The natural events occur independently of one another, and randomly in time.

• tree damage model. The proportion of damaged trees depends on the stand state when the event occurs (in case of storm, see [START_REF] Schmidt | An inventory-based approach for modeling single-tree storm damage experiences with the winter storm of 1999 in southwestern germany[END_REF]). The expected proportion of damaged trees may depend on the time, the tree basal area and the tree height.

• stand valorization model. Two submodels may be considered:

economic model. The economic incomes depend on the timber price functions for thinnings and on the final cutting age.

sequestred carbon model. The sequestred carbon model is related to the instantaneous income and the decreasing following thinning or final cut.

Each model is known with a different accuracy. Growth models are accurate and usually well validated by forestry researchers, while the risk and tree damage models are much more difficult to assess and generally have a low degree of fiability. The economic criteria are usually easy to compute, contrary to the sequestred carbon criteria. It is important to underline that, in a model sequence, the overall accuracy is governed by the less precise model. This justifies the choice made in this first version of EvaSylv to consider tree average-based growth models, instead of modeling basal area structured tree growth which describes more precisely the tree growth.

Silviculture

Two phases are to be considered in silviculture. The first phase consists in respacing. It is user defined by a set of two dates and the corresponding respacing rates (

The criteria used for evaluation

We distinguish the no risk case and the storm risk case.

(i) In absence of storm risk: For a cutting age T , a discount rate δ, the Faustmann

Value J taking into account thinning incomes of a stand is the discounted value of cutting incomes minus cost of regeneration C 1 (n 0 ) and actualized annual costs c a :

J = - c a δ -C 1 (n 0 ) + W e (T ) -C 1 (n 0 ) e δT -1 (1)
where W e (T ) is the total economic income on [0, T ] composed of the sum of thinning incomes R k at time u k actualized at time T and the final income V (T ):

W e (T ) = N k=1 R k .h k e δ(T -u k ) + V (T ) (2)
For a simplified carbon sequestration model, J C = Wc(T ) e δT -1 where:

W c (T ) = δ γ + δ N k=1 Cb(u k ).h k e δ(T -u k ) + Cb(T ) + δ T 0 Cb(u)e δ(T -u) du (3)
(ii) In presence of storm risk:

Let E V (τ ) the expected final income and E Cn (τ ) the expected clearing costs in case of a storm at time τ . E R k is the expected potential thinning income at time u k . The expected Faustmann Value with storm risk becomes:

J R = E[F] = - c a δ -C 1 (n 0 ) + E W (T ) -C 1 (n 0 ) -C 1 (n 0 )a(T ) b(T ) (4)
where

(δ + λ)a(T ) = λ(e (δ+λ)(T -t L ) -1), b(T ) = e (δ+λ)T -λt L -a(T ) -1 and E W e (T ) is
an expected income, with

E W e (T ) = N k=1 β k δ,λ E R k .h k + λ T t L [E V (τ ) -E Cn (τ )]e (δ+λ)(T -τ ) dτ + V (T ) (5)
where

β k δ,λ = β δ,λ (u k ) and β δ,λ (u) = e (δ+λ)(T -u)+λ(u-t L )-.
The expressions of the expected values of criteria are hard-coded in the software.

So, the optimization procedure for expected values does not require simulation of the technical itineraries.

For a simplified Carbon sequestration model, we can express a Faustmann Value F C , the corresponding expected Faustmann value with storm risk

J RC = E[F C ] = E W c (T ) b(T )
where:

E W c (T ) = δ γ + δ N k=1 β k δ,λ Cb(u k ).h k + λ T t L E Cb (τ )e (δ+λ)(T -τ ) dτ + Cb(T ) + δ T 0 β δ,λ (u)Cb(u)du (6)
In the case of an economic evaluation, the software gives both the Faustmann criteria and the Long Run Average Yield, with or without risk. With risk, CV aR(p) is also available.

Optimization

Forest management depends on several control variables (Table 2), some of them being scalars, the others being vectors. The scalar ones include the cutting age T and the final n f tree density at time T . The user has the possibility to optimize the different criteria with the cutting age T imposed or not, the final tree density n(T ) set to n f or only constrained: n(T ) ≥ n f . The respacing dates and rates are control variables. The vector of control variables is organized in two sets. The first set is composed of the thinning dates (u k ) k=1..N such that 0 < u 1 < u 2 < .. < u N < T and the second set is composed of the thinning rates (h k ) k=1..N . A minimum time lag ∆ between two thinnings (starting no sooner than t 0 ) is imposed in the constrained and optimization modes (see Section 2.1).

The results are obtained using algorithms to optimize the different considered criteria.

User selects the choosen algorithm. The provided algorithms do not manage constraints on control variables. In order to satisfy constraints, some artificial variables are introduced (see Appendix A). These artificial variables are required in both the constrained simulation and optimization modes, but are not used in the unconstrained simulation mode.

Software design

EvaSylv is designed to let the user set up the runtime behaviour by modifying scalar values and defining custom functions to be used by the mathematical models. Parameters are organized in a hierarchical tree with nodes. A key point is that custom functions are treated exactly as parameters, giving a great flexibility to the models.

In its present version, the program is a stand alone program, distributed in open source license. The software architecture is schematically composed of three layers At the deepest level, the core layer contains all components of the simulation: model structures and communication between them, technical itineraries and analytical expressions for the criteria. At the most external level, the interface layer handles the (graphical or command line) interface that allows to enter a model and its parameters and to build a scenario. In between, the Build layer fills in the model structures and scenario parameters for transmission to the core layer. The software is entirely written in Python. Python version 3.0 has been used, with the main following Python libraries:

• numpy(www.numpy.org) provides tools and data structures for high precision scientific computation,

• scipy(www.scipy.org), based on numpy provides, among other scientific programming tools, a toolbox for optimization,

• sympy(www.sympy.org) provides symbolic computation tools.

• matplotlib(www.matplotlib.org) is used for plotting functionalities.

The software can be run in script mode with a configuration file or using a graphical user interface (see screenshot in Fig. 2). Parameterizing a simulation/optimization is done by double-clicking an item belonging to this tree and editing it. Customizable functions are compiled before the execution of the simulation/optimization functions.

Users write their own functions using the software interface, which checks syntax validity. These user data are automatically interpreted by the software to build the full them to describe other species. The most important functions are summarized in Table 3.

The function arguments are limited to the ones given in Table 3, with an extra parameter vector, as illustrated in 4.1.1 (named m or p in respective functions).

Models

Results and discussion

EvaSylv offers the possibility to implement many models, depending on the species and on the taking into account of storm risk, or not. The software functionalities are illustrated

with two examples. The first example aims to optimize the economic valorization of a beech stand. The second one is related to the carbon evaluation of a pine stand. In the first case, silviculture is evaluated for the four types of criteria: Faustmann and Average value with risk and without risk. In the second one, the Faustmann criteria are given.

To enter the sub-models, the user can either use the interface or a configuration file, starting from a default file.

Economic evaluation of silviculture for a beech stand

We consider a beech stand. We are interested in the financial incomes obtained from the timber sales.

Growth model

This sub-model is the one that requires the greatest number of functions. The model core is related to the tree basal area. The tree basal area is governed by a dynamic model, following a differential equation. The user must provide the second member of that equation. Let us take as example a growth function [START_REF] Le Moguedec | Fagacees: a tree-centered growth and yield model for sessile oak (Quercus petraea L.) and common beech (Fagus sylvatica L.)[END_REF] with one parameter m 1 and the following mathematical expression:

(1 -e -m1n(t)

√ 4πs(t) )V (t, H 0 )/n(t)
where V (t, H 0 ) is the potential total basal area increase.

The software interface (see Fig. 2) displays a help text with the growth function purpose and its possible dependencies with other model components, in that case the tree density n, the basal area s, the potential height H 0 . The expression of the growth function is given below, written in Python code, very close to the mathematical formulation. The potential total basal area increase is for instance equal to: Here is an example of a H function.

V (t, H 0 ) = p 1 +
def fct_H(self, t, p, H0):

return H0 * ( 1 -exp(-p[1] * t))
The derivative is automatically generated (using the symbolic computation tools of sympy). Hder is given by: def fct_Hder(self, t, p, H0):

return p[1] * H0 * exp(-p[1] * t)
The user may check it and modify it, if necessary.

Optimization

Results given in Table 4 are obtained with the classical Nelder Mead algorithm. The Faustmann value is calculated with a discount rate δ = 0.025 per year. Whatever the chosen criterion, Faustmann value with actualization or Average value, the presence of storm risk reduces the optimal cutting age and therefore the optimal thinning dates occur earlier. This is in agreement with the commonly admitted management strategy.

By changing the timber price function P s (see Table 3), the user can evaluate the criterion behaviour with respect to the market price.

The introduction of risk in the scenarii has a much greater impact onto the Average value. This is due to the fact that the loss is increasing with respect to stand age.

The evolution of n, the tree density, is plotted in logarithmic scale at the right side in Fig. 2.

Evaluation of silviculture for Carbon sequestration of a pine stand

We now consider a pine stand and the impact of silviculture on aerial carbon sequestra- levels. We chose a simple model, where the carbon function is assumed to be proportional to the tree volume, in order to illustrate the interest of the simulation. With or without risk, the criterion to optimize has an integral form (see Eq.3 and Eq.6), so the optimization must be performed with a minimum tree density target.

The software is run with a constraint on the final tree density n(T ) , that must be ≥ n f = 200 stems/ha. Due to the fact that the chosen model for carbon sequestration is an increasing function of time, the optimization will always try to reach the minimum tree density value. The optimal cutting age is found to be 116.4 years, in both cases, with or without risk. The Faustmann criterion value, calculated with a discount rate of 0.03 per year, is found to be equal to 293 tons of CO 2 for the unrisky case, and 273 tons of CO 2 for the risky case. We arbitrarily assumed that a cubic meter of wood stores a ton of CO 2 , according to generally admitted values, that range between half a ton and a ton.

Conclusion

The EvaSylv software presented in this paper hopes to fill a gap in the available softwares for forest management. Even though a lot of forest models are available on tree growth at the stand level, and are included in modeling platforms, they are hard-coded and not user-editable. EvaSylv gives the users an easy way to simulate their own models at stand level, without needing to have a deep knowledge of computer science languages and of the program structure. This is the main strength of the software. EvaSylv functionalities include simulation and optimization tools. Moreover, it allows to take into account natural risks, which is not the case with other forestry simulation platforms.

The software relies upon a sequence of integrated models: growth, risk, tree damage and valorization.

In the present version, limited to storm risk, the software provides two kinds of technical itineraries: a deterministic one (without risk) and a basic one in presence of risk. The modeling in presence of risk requires the calculation of expected values for criteria, which change depending on the technical itinerary. Implementing more complex itineraries requires an interaction with software designers, in order to establish the corresponding evaluation criteria.

The software is likely to participate in the improvement of silviculture guidelines, by taking into account available field data and updated storm risk information.

The approach is flexible enough to open several perspectives. First the elaboration of more complex technical itineraries together with forest managers, which implies mathematical developments to obtain the corresponding analytical criteria. Second the extension to more realistic growth models (structured into classes of basal area). Third the building of ecosystemic service forest-related models.

u k =u k-1 + ∆ + x 2 k N j=1 x 2 j + 1 |x N +1 |, k = 1..N T =t 0 + ∆.N + |x N +1 |, h k =1 -( n f n 0 (1 -h 0 )(1 -h 1 )
)

y 2 k N j=1 y 2 j , k = 1..N
Although the Nelder Mead algorithm has a better behavior than the steepest descent method, it is likely to get stuck in local minima (maxima in our work). The implemented algorithm does not ensure to obtain the solution corresponding to the global maximum.

To avoid this problem, it is recommended to run the optimization with several initial conditions.

The automatic transformation back and forth between artificial and real control variables is implemented in the software workflow, as represented in Fig. 1.

Fig. 1

 1 Fig.1 represents the workflow of the simulation and optimization processes within the

Figure 1 :

 1 Figure 1: Workflow of the proposed modeling approach.

Figure 2 :

 2 Figure 2: EvaSylv graphical user interface.

  def fct_Growth(self, t, n, s, m, H0): g=((1-exp(-m[1]*n*sqrt(4*pi*s)))*self.fct_V(t,H0))/n return gThe growth function definition is prefilled with a default expression. The first line must not be edited. The function takes as argument the m parameter, which is a vector of any length, whose elements will be accessed by using m[1], m[2] . . . in the function body, which can be edited to fulfill the user needs. The components of the m vector are specified by the user in a predefined box located above the corresponding function. The m vector being of variable length, there is a lot of flexibility, as the user can add any parameter of interest. The reference to a pre-defined function must be preceeded by the self. prefix.

Table 1 :

 1 Variables and parameters used for simulation.

	Scope	Variable/Parameter Description
		n	tree density (stems/ha)
	Growth	s	average tree basal area (m 2 )
		H 0	potential height (linked to fertility) (m)
	Risk	λ	risk rate (year -1 )
	Tree damage	t L H L	limit time for tree damage (year ) limit height for tree damage (m)
		L s	damage rate (year -1
		δ	discount rate (year -1 )
	Economic	C n	clearing costs after a destructive event (e)
	parameters	C 1	regeneration cost (e)
		c w	silvicultural cost (e)
		c a	annual cost (e)
	Carbon	γ Cb	carbon discount rate (year -1 ) sequestred carbon (ton)

Table 2

 2

	Scope	Variable/Parameter Description
	Respacing	h 0 , h 1	fixed respacing rates (year -1 )
		td 0 , td 1	fixed respacing dates (year )
		N	maximum number of thinning dates
	Control	h k	kth thinning rate (year -1 )
	variables	u k	kth thinning date (year )
		T	cutting age (free or imposed) (year )
		n f	final tree density (constrained or imposed) (stem/ha)
	Constraints	∆ t 0	minimum elapsed time between thinnings (year ) reference time for thinning date (year )

). The second phase is relative to thinnings. It is also user defined first by the maximum number of thinning dates N , then by the thinning dates (u k ) k=1..N and the corresponding rates of thinnings (h k ) k=1..N . Users must also specify the cutting age T .

Table 2 :

 2 Variables and parameters used in silviculture and optimization.

Table 3 :

 3 Functions used in models.

  p 2 Hder(t, H 0 ), where Hder(t, H 0 ) is the derivative of H(t, H 0 ), the dominant height with respect to time t.

	return p[1]+p[2]*self.fct_Hder(t,H0)
	def fct_V(self, t, p, H0):

Table 4 :

 4 Optimal Value (e/ha) with respect to thinning rates h k , dates u k and cutting age T (years).

	tion. The carbon model is currently an open research topic in itself. Many models exist
	in the literature (Price and Willis, 2011); (Susaeta et al., 2014) with various complexity

Acknowledgments

This research received financial support from the Economic, Human and Social Science Researchers Network of Ecofor, (grant number ECOFOR 2012-24). We thank our funders as well as those who provided assistance with writing this paper.

A Optimization using artificial variables

The constraints on thinnings dates (u k -u k-1 ≥ ∆ years) and, if necessary, the final tree density constraint (n(T ) = n f ) are managed using artificial variables x k and y k . For exemple, in the case without mortality: tration under fluctuating economic environment, forest management and technological changes: An application to forest stands in the southern united states. Journal of Forest Economics, 20(1):47 -64.