Dimitri Justeau-Allaire
email: dimitri.justeau-allaire@cirad.fr

Philippe Birnbaum
email: philippe.birnbaum@cirad.fr

Xavier Lorca
email: xavier.lorca@mines-albi.fr

Unifying Reserve Design Strategies with Graph Theory and Constraint Programming

The delineation of areas of high ecological or biodiversity value is a priority of any conservation program. However, the selection of optimal areas to be preserved necessarily results from a compromise between the complexity of ecological processes and managers' constraints. Current reserve design models usually focus on few criteria, which often leads to an oversimplification of the underlying conservation issues. This paper shows that Constraint Programming (CP) can be the basis of a more unified, flexible and extensible framework. First, the reserve design problem is formalized. Secondly, the problem is modeled from two different angles by using two graph-based models. Then CP is used to aggregate those models through a unique Constraint Satisfaction Problem. Our model is finally evaluated on a real use case addressing the problem of rainforest fragmentation in New Caledonia, a biodiversity hotspot. Results are promising and highlight challenging perspectives to overtake in future work.

Introduction

Human activities are exerting pressure on natural habitats, which generally results in a loss of surface and an increase of fragmentation. As a consequence, many species depending on those habitats are threatened, sometimes with extinction. In this context, it is essential to devote an important part of conservation efforts in the protection of natural habitats through the establishment of nature reserves [START_REF] Beier | Toward Best Practices for Developing Regional Connectivity Maps[END_REF][START_REF] Baguette | Individual dispersal, landscape connectivity and ecological networks[END_REF][START_REF] Haddad | Habitat fragmentation and its lasting impact on Earths ecosystems[END_REF][START_REF] Prendergast | Rare species, the coincidence of diversity hotspots and conservation strategies[END_REF]. Designing a reserve system is a difficult process involving a tradeoff between the conservation targets and the socioeconomic constraints. This problem is known as the reserve design problem. The associated questions are at the crossroad between conservation biology, geography, mathematics, computer science, decision theory and environmental philosophy [START_REF] Sarkar | Environmental philosophy: From theory to practice[END_REF]. In this paper, we focus on the mathematical modeling and the computational solving of the reserve design problem. From this point of view, it is a decision and/or optimization problem. In almost all cases, the combinatorial complexity justifies the need of a systematic approach based on mathematical modeling and computational tools.

In the literature, two major aspects of the reserve design problem usually stand out: the feature covering and the spatial configuration. The first is often referred as the reserve (or site) selection problem [START_REF] Pressey | Beyond opportunism: Key principles for systematic reserve selection[END_REF][START_REF] Revelle | Counterpart Models in Facility Location Science and Reserve Selection Science[END_REF][START_REF] Billionnet | Solving the probabilistic reserve selection problem[END_REF][START_REF] Watts | Marxan with Zones: Software for optimal conservation based land-and sea-use zoning[END_REF]. In extension, we refer to the reserve design problem when spatial attributes are considered [START_REF] Diamond | The island dilemma: Lessons of modern biogeographic studies for the design of natural reserves[END_REF][START_REF] Williams | Spatial attributes and reserve design models: A review[END_REF][START_REF] Billionnet | Designing Connected and Compact Nature Reserves[END_REF][START_REF] Dilkina | Trade-offs and efficiencies in optimal budget-constrained multispecies corridor networks[END_REF][START_REF] Jafari | Achieving full connectivity of sites in the multiperiod reserve network design problem[END_REF]. Current models usually focus on a few aspects of the problem because: (1) they provide an ad-hoc solution to a specific instance of the problem, or (2) they are limited by the modeling paradigm. However, there is a need for a more unified and flexible framework [START_REF] Rodrigues | Flexibility, efficiency, and accountability: Adapting reserve selection algorithms to more complex conservation problems[END_REF] which, in our opinion and based on our experience in New Caledonia, could help to reduce the gap between computer scientists, conservation scientists, and practitioners.

In this paper, we show how the combination of graph-based models with CP can be the basis of such a framework. After a detailed description of the reserve design problem (Section 2), we presents two graph-based models (Section 3). One model is dedicated to the constraint representation of the features covering issues (Section 3.1) and the other one is dedicated to the constraint representation of the spatial issues (Section 3.2). We then unify the models throughout a single CP model based on the Choco constraint solver [START_REF] Prud'homme | Choco Documentation[END_REF] (Section 4). Finally, a realistic operational use case on the problem of rainforest fragmentation in New Caledonia is depicted and first results are discussed (Section 5).

Description of the Problem

The reserve design is a decision and/or optimization problem in the discretized geographical space. Given a set of geographical features (e.g., Fig. 2), we are looking for a reserve system satisfying several criteria, in accordance to conservation targets. In this section, we describe and formalize the problem precisely. We start by defining the characteristics of the problem and then define a set of criteria that can be required for a reserve system.

Characteristics of the Problem -Input Data

The Discretized Geographical Space. The geographical space is tessellated into n granular parcels, which are the decision variables of the problem. Several tessellation methods are possible [START_REF] Sahr | Geodesic Discrete Global Grid Systems[END_REF][START_REF] Birch | Rectangular and hexagonal grids used for observation, experiment and simulation in ecology[END_REF]. The most commonly used is the regular square grid (illustrated in Fig. 1). We choose to this method in this paper.

We denote the number of rows by r, the number of columns by c and the set of parcels by P. We identify a single parcel with the letter i, and index the parcels with integers from 0 to n -1: P = { i | i ∈ 0, n }. While this indexing is not the most convenient for a grid, it has the advantage to be independent of the tessellation method and thus offers extensibility for future work. Finally, we use the 8-connected (cf. Fig. 1) neighborhood to define the adjacency between the parcels, in opposition to the 4-connected neighborhood.

The Environmental Features.

The geographical space is characterized by a set of m environmental features. A feature can be anything that can be spatially represented (e.g. the presence of a species, a certain type of habitat, human constructions). We denote by F the set of environmental features and use the letter j to identify the features:

F = { j | j ∈ 0, m }.

The Values of the Features.

To each feature j is associated a set V j , representing the available data about j among the parcels:

V j = { v ji ∈ R + | i ∈ P }.
Each v ji ∈ V j corresponds to a value describing the feature j in the parcel i. Three types of data are possible: the presence-absence data, the abundance data and the probability of presence data. An example for each data type is given in Fig. 2, and below is a short description for each of them:

-Presence-absence: if j is present in the parcel i, v ji = 1, else v ji = 0. For each (j, i) ∈ F × P we then have v ji ∈ {0, 1}. The presence-absence data is often used to describe the occurrence distribution of a species or a particular characteristic of the landscape (e.g. forest, savanna, fields, roads, cities). -Abundance: in this case, v ji represents a quantitative value about the feature j in the parcel i (e.g. density of trees per parcel, average annual rainfall). For each (j, i) ∈ F × P we then have v ji ∈ [0, +∞[. -Probability of presence: it can be possible to evaluate the probability of presence of a feature j for every parcel in P. The most common situation is the use of Species Distribution Models (SDMs), that are able to combine observations of a species with environmental data to predict its spatial distribution [START_REF] Guisan | Predictive habitat distribution models in ecology[END_REF][START_REF] Elith | Species Distribution Models: Ecological Explanation and Prediction Across Space and Time[END_REF]. With probability of presence data, for each (j, i) ∈ F × P we have

v ji ∈ [0, 1].
The Domains/Anti-Domains of the Features. As illustrated in Fig. 3, to each feature, is associated a set D j and its complement D j . D j represents the domain of j, that is, the parcels where j is present or where the probability of presence of j is not null:

D j = { i ∈ P | v ji > 0 }.
Conversely, D j represents the anti-domain of j, that is, the parcels where j is not present or where the probability of presence of j is null: The domains are represented with solid lines and the anti-domains with dashed lines.

D j = { i ∈ P | v ji = 0 }.

The Reserve System -Solution of the Problem

In the first place, we define the terms "parcel" (sometimes called "site" in the literature), "reserve" and "reserve system". As defined in the previous subsection, a parcel is a granular selection unit of the discretized geographical space. On top of that, a reserve is a set of spatially continuous selected parcels (note that a single selected isolated parcel is a reserve). Finally a reserve system is a set of spatially disjoint reserves (note that a reserve system can be composed of a single reserve). We illustrated the previous definitions in Fig. 4. Given that, a solution to our problem is a reserve system whose attributes are satisfying a set of criteria, themselves depending on the conservation question. We denote such a reserve system by S, its number of reserves by n r and its k th reserve by X k : S = { X k ⊆ P | k ∈ 0, n r } (where P is the set of parcels).

Required Criteria for a Reserve System

According to the underlying conservation questions, several criteria can be required for a reserve system. We distinguish between the feature covering criteria and the spatial criteria.

Feature Covering Criteria. By providing one of the first formalization of the reserve selection problems, ReVelle et al. [START_REF] Revelle | Counterpart Models in Facility Location Science and Reserve Selection Science[END_REF] introduced three fundamental feature covering criteria.

-Covered Features. Among the features that are covered (with certainty) by the reserve system S, we want a set of mandatory features F to be represented (e.g. rare or endangered species). α-Covered Features. Assuming that the v ji 's are pairwise independent, we want a set of features F to be covered by S with a probability of at least α. This criterion is helpful when probability of presence data is available. -k-Redundant Features. A feature j is k-redundant in the reserve system S if and only if it is covered (with certainty) by at least k distinct parcels. We want to enforce this property for a set of features F (e.g. for increasing the chances of persistence of vulnerable species).

Spatial Criteria.

A list of six geometric principles had been defined by Diamond [START_REF] Diamond | The island dilemma: Lessons of modern biogeographic studies for the design of natural reserves[END_REF] and Williams et al. [START_REF] Williams | Spatial attributes and reserve design models: A review[END_REF] summarized them into six spatial attributes to take into account when designing a reserve system: the number of reserves, the reserve areas (by extension we define the reserve system area), the reserve proximity, the reserve connectivity, the reserve shape and core areas and buffer zones. Here we consider three of those spatial attributes (expressed as criteria) and keep the remaining ones for future work.

-Number of reserves. Determining if the best suited is a "single large or several small reserves" (SLOSS), or a "few large or many small reserves" (FLOMS) is a well known debate in ecology ( [START_REF] Diamond | The island dilemma: Lessons of modern biogeographic studies for the design of natural reserves[END_REF][START_REF] Etienne | On optimal size and number of reserves for metapopulation persistence[END_REF]). The conclusion is that the answer strongly depends on the context, and that flexibility is needed. We therefore want to set a minimum value N min and/or a maximum value N max for the number of reserves. -Reserve Areas. Following the previous criterion, it is also essential to provide control on the reserve areas by setting a minimum area A min and a maximum area A max . -Reserve System Area. It should also be possible to express this criterion on the whole reserve system area, by setting a minimum total area A Tmin and a maximum total area A Tmax .

The Graph-Based Models

In this section, we present two graph-based models. The first one is a resource allocation model that will enable us to express the feature covering criteria in the form of constraints. In the same way, the second one is a spatial model that will enable us to control the spatial criteria. In both models, each parcel is represented by a vertex. This common characteristic is essential since it is the one that makes the aggregation of the models possible, through a set of appropriate channeling constraints.

The Resource Allocation Graph

We consider parcels as resources that can be allocated to the conservation of features, then considered as tasks, and thus define the directed graph G r = (V r , A r ), also called the resource allocation graph. The vertices of G r are partitioned into three disjoint sets F r , P r and {s, t}. F r represents the feature (or task) vertices, P r represents the parcel (or resource) vertices, s is the source vertex and t the sink vertex.

V r = F r ∪ P r ∪ {s, t};

F r = { f j | j ∈ F }; P r = { p i | i ∈ P }. (1)
Furthermore, using A r (X, Y ) as the notation for the set of all X-Y arcs, we define the arcs of G r in the following way:

A r = A r (s, F r ) ∪ A r (F r , P r ) ∪ A r (P r , t).
(

A r (s, F r ) and A r (P r , t) are defined such that there is an arc from s to each feature vertex and an arc from each parcel vertex to t:

A r (s, F r ) = { (s, f j ) | f j ∈ F r }; A r (P r , t) = { (p i , t) | p i ∈ P r }. (3)
Moreover, A r (F r , P r ) represent the possible allocations between F r and P r . More precisely, there is an arc from a feature vertex f j to a parcel vertex p i if and only if the feature j is represented in the parcel i, that is i ∈ D j . We then have:

A r (F r , P r ) = j∈F { (f j , p i ) | i ∈ D j }. (4)
On the arcs of G r , we define a lower bound (or demand) function l :

A r → R + ∞
and an upper bound (or capacity) function u :

A r → R + ∞ such that if f is a flow in G r : ∀a ∈ A r , l(a) ≤ f (a) ≤ u(a). ( 5
)
Finally, we define H r : P(P r ) → P(V r ) × P(A r ) that associates to a set X ⊆ P r the subgraph of G r induced by {s, t} ∪ F r ∪ X, that is, the resource allocation graph obtained when only considering a subset of parcels. We denote by V r [H r (X)] the vertices of H r (X) and by A r [H r (X)] the arcs of H r (X). An example is provided in Fig. 5.

H r (X) = G r [ {s, t} ∪ F r ∪ X ]. ( 6
) X = {p 2 , p 8 } p8 F r s f0 t f2 p2 f1 l(s,f 2 ) u(s,f 2 ) l(f 0 ,p 8 ) u(f 0 ,p 8 ) l(f 2 ,p 8 ) u(f 2 ,p 8 ) l(f 1 ,p 2 ) u(f 1 ,p 2 ) l(p 2 ,t) u(p 2 ,t) l(f 2 ,p 2 ) u(f 2 ,p 2 ) l(s,f 1 ) u(s,f 1 ) l(s,f 0 ) u(s,f 0 ) l(p 8 ,t) u(p 8 ,t)
Fig. 5. Hr({p2, p8}) associated with the example in Fig. 2. Lower and upper bounds are represented on the arcs.

Expressing Feature Covering Criteria as Constraints. From this point, to each feature covering criterion (as defined in the previous section) we associate a constraint that can be applied on the resource allocation model. More precisely, if S is a reserve system and X s its associated set of parcel vertices, a criterion is satisfied by S if and only if its associated constraint if satisfied by H r (X s ). We express these constraints as flow constraints by defining the value of l and u on certain arcs. When the value of l is not explicitly defined, it is unconstrained and then set to 0. Similarly, u is set to +∞ when its value is not explicitly defined. Covered Features. In our resource allocation model, we can easily express this criterion as a flow constraint on H r (X s ).

Constraint 1: Covered Features.

Input parameter(s):

A set of features F ⊆ F .

The set of features F is covered by S if and only if Hr(Xs) admits a feasible flow f verifying (5) when:

     l(s, fj) = 1, ∀fj ∈ F r ;
l(fj, pi) = 1, ∀(fj, pi) ∈ Ar(Fr, Pr); u(fj, pi) = 1, ∀(fj, pi) ∈ Ar(Fr, Pr) such that vji ≥ 1; u(fj, pi) = 0, ∀(fj, pi) ∈ Ar(Fr, Pr) such that vji < 1.

α-Covered Features. To express this criterion, we assume that the probabilities of presence v ji are pairwise independent. We then rely on the probability of absence q ji = (1 -v ji ) and express the constraint as:

∀j ∈ F , i ∈ S q ji ≤ 1 -α. ( 8
)
We then express the α-presence constraint in the following way:

Constraint 2: α-Covered Features.

Input parameter(s):

A set of features F and a real α ∈ [0, 1].

The set of features F is covered by S with a probability of at least α if and only if Hr(Xs) admits a feasible flow f verifying (5) when:

l(s, fj) = -log(1 -α), ∀fj ∈ F r ; u(fj, pi) = -log(qji), ∀(fj, pi) ∈ Ar(Fr, Pr) such that vji < 1. (9)
k-Redundant Features. Since the k-redundancy is actually a generalization of the covering features criterion, we can also express it as a flow constraint on H r (X s ).

Constraint 3: k-Redundant Features.

Input parameter(s):

A set of features F and a positive integer k.

The k-redundancy of the set of features F in the reserve S is satisfied if and only if Hr(Xs) admits a feasible flow f verifying (5) when:

     l(s, fj) = k, ∀fj ∈ F r ; l(fj, pi) = 1, ∀(fj, pi) ∈ Ar(Fr, Pr); u(fj, pi) = 1, ∀(fj, pi) ∈ Ar(Fr, Pr) such that vji ≥ 1; u(fj, pi) = 0, ∀(fj, pi) ∈ Ar(Fr, Pr) such that vji < 1. ( 10
)

The Spatial Graph

We now define the undirected graph G s = (V s , E s ), the spatial graph, which is a representation of the discretized geographical space P (a r × c regular square grid in our case). Once again, to each parcel i of P, we associate a vertex p i , we then have:

V s = { p i | i ∈ P }. (11)
Moreover, the edges of G s are defined such that if p u and p v are two vertices, there is an edge between p u to p v if and only if the parcels u and v are spatially adjacent. The edges of G s can be partitioned into four disjoint sets: the horizontal edges (E H ), the vertical edges (E V ), the north-west to south-east diagonal edges (E N W S E ) and the north-east to south-west diagonal edges (E N E S W ).

E s = E H ∪ E V ∪ E N W S E ∪ E N E S W ; E H = { (p i , p i+1 ) | i ∈ P ∧ ¬ (i + 1) ≡ 0 (c) }; E V = { (p i , p i+c ) | i ∈ P ∧ i < c(r -1) }; E N W S E = { (p i , p i+c+1 ) | i ∈ P ∧ i < c(r -1) ∧ ¬ (i + 1) ≡ 0 (c) }; E N E S W = { (p i , p i+c-1 ) | i ∈ P ∧ i < c(r -1) ∧ ¬ i ≡ 0 (c) }. ( 12
)
See Fig. 6 for an illustration of the above equation. Also note that it takes into account the extremal positions of the grid. In fact, the parcels located in the first column are the one whose index is a multiple of c, that is i ≡ 0 (c). Moreover, the parcels located in the last column are the ones preceding those that are located in the first column, that is (i + 1) ≡ 0 (c). Finally, the parcels located in the last line are the ones satisfying i < c(r -1).

pi+c pi+c+1 pi pi+1

Fig. 6. Illustration of a portion of a spatial graph Gs associated with a r × 4 square grid, using the 8-connectivity neighborhood definition.

Expressing Spatial Criteria as Constraints.

Similarly to what had been defined for the resource allocation graph, to a solution S of the problem we associate X s ⊆ V s . Moreover, to each reserve X k ∈ S we associate X s (k) ∈ X s , the vertices associated to the parcels of X k . We now express each spatial criterion as a constraint that can be applied on G s [X s ].

Number of Reserves.

We easily express this criterion by bounding the number of connected components (NCC, [START_REF] Dooms | The CP(Graph) Computation Domain in Constraint Programming[END_REF][START_REF] Beldiceanu | Graph Invariants as Necessary Conditions for Global Constraints[END_REF][START_REF] Beldiceanu | Graph Properties Based Filtering[END_REF]) in G s [X s ]. Ensuring that the number of reserves in S is bounded by Nmin and Nmax is equivalent to bounding the NCC of Gs [Xs] with Nmin and Nmax.

Nmin ≤ NCC(Gs[Xs]) ≤ Nmax. ( 13
)
Reserve Areas. We express this criterion as a constraint on the number of vertices of the smallest connected component of G s [X s ] (MIN NCC, [START_REF] Beldiceanu | Global Constraint Catalog, 2nd Edition[END_REF][START_REF] Beldiceanu | Graph Invariants as Necessary Conditions for Global Constraints[END_REF][START_REF] Beldiceanu | Graph Properties Based Filtering[END_REF]) and on the number of vertices of the largest connected component of G s [X s ] (MAX NCC, [START_REF] Beldiceanu | Global Constraint Catalog, 2nd Edition[END_REF][START_REF] Beldiceanu | Graph Invariants as Necessary Conditions for Global Constraints[END_REF][START_REF] Beldiceanu | Graph Properties Based Filtering[END_REF]).

Use case: Rainforest Fragmentation in New Caledonia

New Caledonia is biodiversity hotspot located in the South Pacific, slightly north of the tropic of the Capricorn. The flora of this large archipelago is distinguished by an exceptionally high rate of endemism. Like most of the world's remaining natural forests, New Caledonian rainforests are endangered with surface loss and fragmentation. A case study had been conducted in the south of New Caledonia in order to highlight "how does forest fragmentation affect tree communities" [START_REF] Ibanez | How does forest fragmentation affect tree communities? A critical case study in the biodiversity hotspot of New Caledonia[END_REF]. We relied on this case study and its associated dataset (up to date) for our use case, and considered the following fictive but realistic operational scenario:

"We want to establish a reserve system in which a pool of endangered species must be present. In addition, most of the other species known in the area must have a high probability to occur, or a high habitat suitability. The reserve system must be mostly covering rainforest areas. Its area and its number of reserves must be limited because of budget limitation. Moreover, each reserve must be large enough to ensure the persistence of the species." Note. In this scenario, the objective is to protect both existing and potential rainforest areas. To do so, we relied on SDM layers that were generated with presence-only data and thus produce a score of habitat suitability rather than a standardized probability of presence. A high habitat suitability in a nonrainforest zone can then be interpreted as an adequate zone for recolonization.

Input Data, Constraints and Parameters

The original dataset consists of the mapping of a 60 km 2 landscape where 97 tree communities had been sampled in 88 digitized rainforest fragments (forest/nonforest). The dataset gathers 5431 identified trees belonging to 223 species. Moreover, an SDM raster layer was available for 173 of the species [START_REF] Pouteau | Accounting for the indirect area effect in stacked species distribution models to map species richness in a montane biodiversity hotspot[END_REF][START_REF] Schmitt | Ssdm: An r package to predict distribution of species richness and composition based on stacked species distribution models[END_REF]. Arbitrarily, we considered the 50 species without SDM as the endangered ones. We then prepared this dataset by tessellating the study area into a 46 × 75 regular square grid and by rasterizing the dataset according to this grid. Each parcel then has an area of about 1.7 ha. Note that we also defined a set of forbidden parcels corresponding to lakes and mining sites.

From this point, we defined a feature for each observed species in the area. When available, we relied on the SDM layer for the feature data (probability of presence data). We forced the values to 1 for the parcels where an observation is available. When no SDM was available, we only relied on the occurrence dataset (presence-absence). We represented the rainforest coverage as a presence-absence feature.

We then applied the Covered Features constraint for the set of endangered species, and the α-Covered Features constraint for the other species with α = 0.8. In order to ensure a minimum rainforest area of 340 ha in the reserve system, we applied the k-Redundant Features constraint for the rainforest coverage, with k = 200 parcels. Moreover, we enforced the forbidden parcels on the envelope of the graph variable g. We then set the minimum area of the reserves to A min = 176 parcels (about 300 ha) with the Reserve Areas constraint. In addition, we limited the reserve system area using the Reserve System Area constraint, with A tmax = 589 parcels (about 1000 ha). According to those restrictive parameters, we allowed the number of reserve to be at most two, using the Number of Reserves constraint, with N min = 1 and N max = 2.

Questioning and Results

In the first place, the number of reserves and the number of parcels are critical parameters of our use case: the less the better. This is why we started by implementing a search strategy that starts by branching on the lower bound of the nbCC variable and continues by selecting the lower bound of the parcels variables, sorted in descending order by a score corresponding to the number of features with a value greater than 0.6 (cf. 7). The solver quickly found a solution to the decision problem (DP), as summarized in Table 1. Given that, we ran a first optimization problem (OP1) where we tried to minimize the total area of the reserve system, that is the nbParcels variable. We limited the computation time to 4 hours and retrieved the best solution found, which reduced the total area by 8% in comparison to DP (cf. Table 1). In order to cover more rainforest parcels, we ran a second optimization problem (OP2) in which we forced the nbParcels variable to be within 15% of the best value found in OP1, and tried to maximize k (the number of forest parcels), thus defined as an integer variable. After a limited run of 4 hours, we could increase the area of rainforest by 12% (cf. Table 1). A mapping of our results is provided in Fig. 7.

Conclusion and Challenges

To the best of our knowledge this paper tackles, for the very first time, the reserve design problem from a constraint programming point of view. It is also the first time that a reserve design model integrates such a diversity of constraints, simultaneously involving decisions based on occurrences, SDMs and spatial attributes with an exact approach. Although performance enhancements are needed, the combination of graph-based modeling and constraint programming reveals as a powerful and promising framework for dealing with the reserve design problem.

Based on a challenging use case, our model highlighted a solution compatible with the conservation strategy, namely a trend to link isolated forest patches in order to enhance the functioning of tree communities. However, in this use case we restrained to a binary landscape only composed of forest/non-forest while it is often assumed that a reserve system must include an assemblage of several landscape types. In such a mosaic, an important challenge lays in weighting and balancing the reserve system characteristics and shape in order to maintain (or restore) the functional connectivity inside and between the reserves. In fact, the functional isolation of an habitat leads to a reduction of biological flows, which tends to amplify its spatial isolation. Moreover, since the underlying processes are dynamic, robust solutions must rely both on the current state and future scenarios. It also remains to model the impacts of a reserve system on the offreserve area, such as the creation of boundaries or enclosed areas.

These elements lead us to identify several lacks and challenges. First of all, the main lack is that Choco solver does not offer an implementation of the flow constraint [START_REF] Steiger | An Efficient Generic Network Flow Constraint[END_REF][START_REF] Downing | Integration of AI and OR Techniques in Contraint Programming for Combinatorial Optimzation Problems[END_REF]. We will focus on its implementation in future work. Next, a bottleneck in the constraint propagation is the interaction between the constraint on the number and the size of the connected components [START_REF] Beldiceanu | Global Constraint Catalog, 2nd Edition[END_REF]. We actually treat each one independently but we think that there is a possible enhancement of the filtering by dealing with their interaction. A first challenge for future work concerns our capacity to model constraints on the shape of the reserves by using graph properties, such as graph diameter, in order to design reserves that are compatible with the long-term persistence of species. A second challenge is more oriented to decision making aspects such as identifying key areas that have to be present in any solution. Finally, a last challenge is related to our capacity to take a dual point of view: is it possible to take into account managers' needs on the off-reserve area by adding constraints on the same graph representation?

Fig. 1 .Fig. 2 .

 12 Fig. 1. Square grid tessellation and 8-connected neighborhood illustrations.

Fig. 3 .

 3 Fig.3. The values, domains and anti-domains associated with the features in Fig.2. The domains are represented with solid lines and the anti-domains with dashed lines.

Fig. 4 .

 4 Fig. 4. Illustration of a reserve system, composed of four reserves, themselves made of several adjacent parcels.

Constraint 4 :

 4 Number of Reserves.Input parameter(s):Two positive integer Nmin and Nmax.

Fig. 7 .

 7 Fig. 7. Mapping of the use case best solutions. The parcel scores correspond to the heuristic score and the mandatory occurrences to rare species observed only once in the study area (they must then be covered by any solution).

Table 1 .

 1 Use case results: resolution times and solution characteristics. All experiments were run on an Intel Core i5-5200U CPU (2.20GHz×4), with 7.7GB of RAM.

		DP	OP1	OP2
	Resolution time	28s	3h24m	1h5m
	Number of solutions found	1	8	25
	Number of reserves	1	1	1
	Number of parcels	318	292	328
	Number of rainforest parcels	200	200	224

https://gist.github.com/dimitri-justeau/8098af35824bbf8d52ef21282291e621

Constraint 5: Reserve Areas.

Input parameter(s):

Two positive integer Amin and Amax.

Ensuring that the area of every reserve X k ∈ S is bounded by Amin and Amax is equivalent to constraining the lower bound of MIN NCC(Gs [Xs]) to Amin and the upper bound of MAX NCC(Gs[Xs]) to Amax. ∀k ∈ 0, nr , MIN NCC(Gs[Xs]) ≥ Amin; MAX NCC(Gs[Xs]) ≤ Amax. [START_REF] Jafari | Achieving full connectivity of sites in the multiperiod reserve network design problem[END_REF] Reserve System Area. In the current case of a regular tessellation method, we can control the whole reserve system's area by bounding the norm of X s .

Constraint 6: Reserve System Area.

Input parameter(s):

Two positive integer AT min and AT max .

Ensuring that the total area of the reserve system is bounded by AT min and AT max is equivalent to bounding |Xs|.

The CP Model

In this section we present our CP model for the reserve design problem. For its implementation, we rely on the solver Choco [START_REF] Prud'homme | Choco Documentation[END_REF] and its extension Choco-graph [START_REF] Fages | Choco Graph Documentation[END_REF], which provides graph variables and constraints.

The Decision Variables. We naturally model the parcels with a boolean variable array, named parcels. If the parcel i is selected in the reserve system, parcels[i] = 1, else parcels[i] = 0.

BoolVar[] parcels = model.booVarArray("parcels", n);

These decision variables are the cornerstone of our CP model because they allow us to aggregate the two models we introduced in the previous section.

The Feature Covering Constraints. Given the particular configuration of the resource allocation graph, we are able to express each feature covering constraint with several local flow conservation inequalities, one for each feature involved in the constraint. Note that we would certainly benefit from the filtering of a global flow constraint [START_REF] Bockmayr | Network Flow Problems in Constraint Programming[END_REF]. However, there is no such constraint implemented in Choco at the time we are writing this paper. We thus keep this idea for future work. Constraint 1, Covered Features [START_REF] Revelle | Counterpart Models in Facility Location Science and Reserve Selection Science[END_REF]: with local flow conservation inequalities, [START_REF] Revelle | Counterpart Models in Facility Location Science and Reserve Selection Science[END_REF] becomes:

Below is the implementation with Choco 4, using the scalar constraint.

for (int j : featuresToCover) { int[] coeffs = Arrays.stream(V[j]) .mapToInt(v -> (v >= 1) ? 1 : 0) .toArray(); model.scalar(parcels, coeffs, ">=", 1).post(); } Constraint 2, α-Covered Features ( 9): the coefficients in the scalar constraint must be integers. We then retain only two digits of precision for the probabilities of presence. If α ∈ [0, 0.99] thenlog(1-α) ∈ [0, 2], moreover, with this precision the order of the smallest variation between two values (α = 0 and α = 0.01) is 10 -3 , multiply our local flow inequality by 10 3 in order to stay in the integer domain. If v ji ≥ 1, we set the flow upper bound to -10 3 log(1 -0.999) = 3000 as a replacement for +∞ .Consequently, we reduce (9) to:

Below is the implementation with Choco 4.

for (int j : featuresToCover) {

int[] coeffs = Arrays.stream(V[j]) .mapToInt( v -> (v >= 1) ? 3000 : (int) (-1000 * Math.log10(1 -v))) .toArray(); int scaled = (int) (-1000 * Math.log10(1 -alpha)); model.scalar(parcels, coeffs, ">=", scaled).post(); } Constraint 3, k-Redundant Features [START_REF] Diamond | The island dilemma: Lessons of modern biogeographic studies for the design of natural reserves[END_REF]: similarly, we reduce [START_REF] Diamond | The island dilemma: Lessons of modern biogeographic studies for the design of natural reserves[END_REF] to:

And implement it the following way with Choco 4:

.mapToInt(v -> (v >= 1) ? 1 : 0) .toArray(); model.scalar(parcels, coeffs, ">=", k).post(); } The Spatial Constraints. We rely on Choco-graph to express the spatial constraints in our CP Model. First, we use a graph variable g to model the reserve system. Its kernel is the empty graph (GLB in the code), and its envelope is G s (GUB in the code).

UndirectedGraph GLB = new UndirectedGraph(model, n, BIPARTITESET, false); UndirectedGraph GUB = new UndirectedGraph(model, n, BIPARTITESET, false); for (int i = 0; i < n; i++) { GUB.addNode(i); for (int ii : getNeighbors(i)) { GUB.addEdge(i, ii); } } UndirectedGraphVar g = model.graphVar("g", GLB, GUB);

Then, we link the graph variable g with the boolean variables parcels using the nodesChanneling constraint.

model.nodesChanneling(g, parcels).post();

We also force the existence of an edge between two selected adjacent parcels through an edgeChanneling constraint with a reified and constraint between each pair (i1, i2) of adjacent parcels. Doing so, we ensure that every existing edges between two selected vertices are also present in our graph variable. IntVar nbCC = model.intVar("nbCC", Nmin, Nmax); model.nbConnectedComponents(g, nbCC).post(); Constraint 5, Reserve Areas [START_REF] Jafari | Achieving full connectivity of sites in the multiperiod reserve network design problem[END_REF]: at the time we are writing this paper, there is no constraint in Choco-graph for controlling the MIN NCC and MAX NCC graph properties. We thus implemented the sizeConnectedComponents 5 constraint, which allows us to bound MIN NCC and MAX NCC.

model.sizeConnectedComponents(g, Amin, Amax).post(); Constraint 6, Reserve System Area [START_REF] Rodrigues | Flexibility, efficiency, and accountability: Adapting reserve selection algorithms to more complex conservation problems[END_REF]: we can control the number of vertices of G s (that is, the number of parcels) through the nbNodes graph constraint, or through the sum constraint over the decision variables.

IntVar nbParcels = model.intVar(Atmin, Atmax); model.nbNodes(g, nbParcels).post(); // Option 1 model.sum(parcels, "=", nbParcels).post(); // Option 2


Laurent Robert

Gilles Dusserre

Assessment of Thermoset Cure-Induced Strains by Fiber Bragg Grating Sensor

Composites based on thermosetting resin and reinforcement fibers present generally unwanted residual internal stresses inherent to their elaboration process. In this article, different curing experiments of thermosets (isothermal and anisothermal) were monitored using optical fiber Bragg grating (FBG) sensors and thermocouples, to assess the cure-induced strains. A thermokinetic model evaluates the degree of conversion of the resin. At the onset of stress transfer to the optical fiber, the degree of conversion ranges between 0.63 and 0.68. During curing, the FBG deforms under chemical shrinkage with an amplitude widely lower than the expected chemical strain, showing that the FBG signal is not directly related to the actual chemical shrinkage. However, once the resin is cured, the FBG sensor provides directly the coefficient of thermal expansion of the resin, as function of temperature and degree of conversion, reached in the different experiments.

INTRODUCTION

Fiber-reinforced composites are generally used as structural materials because of their excellent specific strength. However, they can suffer for sensitivity to damage and impact (delamination) [1], and their failure properties are not always well predicted. This is caused partially by residual internal stresses [2] and flaws, originating from the elaboration process (porosity, curing defect). The process can be monitored by including adequate high performance in-situ sensing systems [3], allowing a better knowledge of the process cycle and the related parameters.

An important issue that has promoted the development of in-situ composite process instrumentation is the knowledge of cure residual stresses (CRS) [2][3][4]. Indeed, process parameters have an important impact on the properties of the composite structure, and particularly in the development of CRS. CRS can cause important warpage, they can affect the final material properties and have an influence on the sensitivity to damage. The fine knowledge of the development, characterization, and amplitude of these CRS is thus necessary to reduce their effects thanks to an optimization of the cure cycle by adjusting the process parameters.

CRS have three main origins: chemical, thermal, and mold-laminate interaction origins. The chemical origin comes from the network structure of the thermoset which is built during polymerization. The reaction of reticulation and the formation of these networks induce volume variations of the matrix usually named chemical shrinkage. In a fiber-reinforced composite, fibers are opposed to the matrix chemical contraction, they are thus submitted to compressive stresses whereas the matrix is submitted to tensile stresses. During reticulation, at temperature above the glass transition, the viscosity of the network allows the relaxation of an important part of the internal stresses. Nevertheless, the level of chemical contraction of high performance epoxy systems remains relatively low, for example, volumetric shrinkage of RTM6 is approximately 2.7% [5]. However, very few studies report quantitative values of stress transfer onset on laminate fibers as well as chemical shrinkage level transferred to the fibers.

The CRS from thermal origin comes from the mismatch between the coefficients of thermal expansion (CTE) of all the involved components in presence, as manufacturing structural composites implies curing at a relatively high temperature (typically around 180 C). At the scale of the components, differential CTE between fibers and matrix produces internal stresses in the composite material. At the scale of the laminate, CRS are induced because of the anisotropic character of plies in terms of expansion behavior, stresses appearing mainly during the cooling to room temperature.

The third cause of CRS is linked to interactions between mold and composite caused by their differential thermoelastic behaviors. Stresses due to differences of CTE between tools and reinforcement give rise to strain in the adjacent reinforcement ply and gradually to the other plies.

From two decades, it has been reported the use of several optical fiber sensors for the measurement of residual strains in composite, most often extrinsic Fabry-Perot interferometric sensor [6] and gratings-based sensors [2][3][4][7][8][9]. Optical fiber sensors offer good solutions to address measurements of internal strains because of their intrinsic advantages: they are small and flexible, light, immune to electromagnetic interference, and can withstand corrosive environments and high temperatures. One of the most interesting candidate is the fiber Bragg grating (FBG) sensor, which can provide localized measurements with high sensitivity and good linearity. Moreover, strain gradients can be measured using the optical low-coherence reflectometry technique [10].

An important point is to know how and when the resin can transfer a load to the reinforcement fibers, giving rise to the onset of CRS. For instance, it has been reported in [11,12], using fiber pull-out tests, that the onset of gelation is the key point of the stress transfer between fibers and matrix. Using embedded optical fiber sensors during composite processing is a great advantage, because even if the diameter of the optical fiber (125 mm) is larger than the reinforcement fibers, the matrix may behave in the same way in the vicinity of both reinforcement fibers and optical fibers. The latter is thus and in a similar way submitted to induced strains, and can give information about the beginning, the origin, and the intensity of the CRS.

In this article, both chemical and thermal cure-induced strains are investigated from the measurement of temperature (thermocouples), strain (FBG sensors), and degree of conversion (thermokinetic model). The occurrence of the stress transfer is highlighted during the cross-linking of a thermoset resin subjected to several thermal cycles, mainly isothermal reticulation, cooling at a fixed degree of conversion, and anisothermal reticulation. The next section describes the materials and methods used for this work. The results are shown in section RESULTS for isothermal and anisothermal curing experiments and discussed in the last section.

EXPERIMENTAL METHODS AND MODELIZATION

FBG Sensor

A FBG is a local permanent submicron period refractive index modulation in the core of an optical fiber, created by an appropriate sideways illumination with UVlight of a photosensitive optical fiber [START_REF] Othonos | Fiber Bragg Gratings: Fundamentals and Applications in Telecommunications and Sensing[END_REF]. When light with a broadband spectrum travels into the optical fiber and interacts with the grating, it acts to couple the fundamental forward propagating mode to the contrapropagating core mode. This coupling occurs at the specific socalled Bragg wavelength k B , with the FBG acting as a narrow-band reflection filter. The Bragg wavelength k B is defined by the grating period K (typically 0.5 mm) and the effective refractive index of the propagating mode n eff (1.4682 for a Corning SMF-28 TM fiber at 1550 nm), according to the well-known Bragg equation, Eq. 1.

k B 52n eff K (1)
A change in either one of these parameters changes the wavelength. A temperature increase causes a thermal expansion of the grating (and thus a change of its period) and also a change in the refractive index. A mechanical strain deforms the FBG and thus changes its period, and modifies also the refractive index by photoelastic effects.

In a general point of view, let us consider the case of a FBG embedded in a host material that submits it to a three-dimensional (3D) thermomechanical stress field. Under the assumptions of uniform strain field along the grating and negligible shear and bending, the strains of the fiber core are the same as the strains of the optical fiber. The relationship that expresses the wavelength shift versus the three principal strains e x , e y , e z (z is the direction of the fiber axis) and the temperature change DT is as follows [START_REF] Tao | [END_REF]15]:

Dk B;x k 0 B 5e à z 2 n 2 eff 2 p 11 e à x 1p 12 e à y 1e à z h i 1CTE Á DT1n Á DT Dk B;y k 0 B 5e à z 2 n 2 eff 2 p 11 e à y 1p 12 e à x 1e à z À Á h i 1CTE Á DT1n Á DT 8 > > > > < > > > > : (2)
where the strain components e x , e y , e z are decomposed into a sum of the thermal expansion and the mechanical strains e Ã

x , e à y , e à z induced by stress. In Eq. 2, k 0 B is the initial Bragg wavelength, CTE 5 1 K @K @T is the coefficient of thermal expansion of the optical fiber (0.55 3 10 26 K 21 for silica), n5 1 n eff @n eff @T represents the thermo-optic coefficient (7.3 3 10 26 K 21 for SMF-28 TM silica-core fiber), and p 11 and p 12 are the Pockel's strain-optic coefficients (p 11 5 0.113 and p 12 5 0.252 [START_REF] Tao | [END_REF]). Note that the index change is by far the dominant thermal effect. As shown by Eq. 2, due to the birefringence effect, the Bragg peak is split into two distinct ones. In the present study, we consider only an axisymmetric configuration and thus e r 5e x 5e y is the radial strain, and e z is the axial strain. In this case, both peaks are superimposed and the Bragg wavelength variation is given by Eq. 3:

Dk B k 0 B 5e à z 12 n 2 eff 2 p 12 1 p 11 1p 12 ð Þ e à r e à z ! & ' 1CTE Á DT1n Á DT (3)
The thermal contraction/dilatation of a host material in which the FBG is embedded, is here included in the mechanical strain e à as it results in an occurrence of thermomechanical stresses. The information that is expected from the FBG measurement is the axial mechanical strain e à z that reveals the state of stress of the fiber in its environment. It is therefore mandatory to know the temperature change DT and the biaxiality ratio defined as B5e à r =e à z . A first issue appears when the two solicitations, mechanical strain and temperature change, coexist simultaneously. Many publications have been devoted to the discrimination of the temperature and strain responses of FBG sensors and a lot of techniques have been reported to overcome this problem [START_REF] Othonos | Fiber Bragg Gratings: Fundamentals and Applications in Telecommunications and Sensing[END_REF]16,17]. The most classical extrinsic scheme consists of two FBGs, with one encapsulated to remain free from mechanical stresses [8]. Optical intrinsic schemes are also reported that ensure an optimum spatial resolution as the measurement is performed exactly at the same location, for example, by superimposing a FBG and a long-period grating [9,18]. A thermocouple can also be associated to a FBG to collect the thermal component [4]. This last method was applied in the present study. From Eq. 3, a thermal calibration by keeping the FBG free from mechanical stresses allows to evaluate the thermal sensitivity K T 5k 0 B CTE 1n ð Þ % 12 pm C 21 at 1550 nm [START_REF] Othonos | Fiber Bragg Gratings: Fundamentals and Applications in Telecommunications and Sensing[END_REF]18].

The second issue comes from the knowledge of the biaxiality ratio B5e à r =e à z , function of the thermomechanical loading, that can only be assessed by considering a particular configuration of the fiber embedded in the host material. In our case, we consider a host material (i.e., the epoxy resin) subjected to isotropic volume change due to chemical shrinkage and thermal contraction. An analytical model is proposed in the following to evaluate B. This model involves an infinite straight optical fiber (radius R, material parameters: Young's modulus E, Poisson's ratio m, coefficient of thermal expansion CTE) embedded in an infinite coaxial resin cylinder (radius R e , material parameters: Young's modulus E r , Poisson's ratio m r , CTE r ). The problem is axisymmetric and invariant by translation along the fiber axis z. The momentum balance of the system subjected to a temperature change is solved by assuming the continuity of the displacements and of the radial stress at the interface between the two materials, and a free external boundary of the resin. Results allow to calculate the radial strain e r and the axial strain e z function of the temperature change DT, of the resin cylinder R e and of the resin parameters E r , m r and CTE r , the parameters relative to the optical fiber being known and given in Table 1.

The first remark is that the biaxiality ratio B5e As expected, for an infinite amount of resin surrounding the optical fiber and an infinite small elastic modulus, B tends toward 2m. In this particular case, Eq. 3 can be simplified and gives the strain sensitivity of a naked FBG, K e % 1.23 pm (me) 21 at 1550 nm [START_REF] Othonos | Fiber Bragg Gratings: Fundamentals and Applications in Telecommunications and Sensing[END_REF]18]. In all other cases, B should be evaluated correctly knowing the actual values of the parameters of the surrounding host material. On the one hand, it remains very complex to process the FBG signal by taking into account the variation of B during the experiment using the results of a thermomechanical modeling. On the other hand, it is very easy to calculate the axial strain error committed when assuming a constant B ratio. The axial error err(B) plotted in Fig. 1b function of R e for various resin elastic modulus is defined by Eq. 4.

err B ð Þ5 e à z B ð Þ2e à z m ð Þ e à z B ð Þ Ã 100 (4)
The actual configuration of the experiments involves a 1-mm resin radius (thickness of the sample higher than 2 mm), and a resin Young's modulus E r varying between 0 and 3.5 GPa along with the polymerization. The value of B may thus vary between 20.13 and 20.17 for the whole set of experimental data under study (Fig. 1a). It was chosen to consider B 5 20.17 (fiber Poisson's ratio) that leads to a negligible error when the resin modulus is close to zero. The maximal error in the measured axial strain is lower than 2% (Fig. 1b) and is obtained when the resin is fully cured. Note that this error level is acceptable in many cases (typically for isotropic surrounding material with relatively low modulus compared to the one of the optical fiber), but in some cases considering B as the optical fiber Poisson's ratio can induce very erroneous data for instance when the host material is anisotropic as in the case of a composite material [19].

Materials and Experimental Methods

The resin under study is the HexFlow V R RTM6 from Hexcel V R Composites. RTM6 is a premixed epoxy resin specifically dedicated to liquid composite molding and commonly used industrially to manufacture aeronautic carbon-reinforced composite parts.

Three kinds of mold were used to analyze the effect of the mold interaction on the resin cure-induced strain as well as to enhance or not the exothermicity of the reaction. A rigid metallic (RM) mold, which consists basically of a 0.9-mm thick 120 3 60 mm 2 steel plate, was used as a standard material for mold. It was treated with a mold release agent (Loctite V R Frekote V R ). A soft metallic (SM) mold was made of a simple 0.05-mm thick circular ($80mm diameter) aluminum foil pan, thus easily deformable, and enabling large thermal transfer. The thickness of the resin samples was, respectively, of 2.85 and 2.4 mm in these two first molds. The last mold was a 1-mm thick silicon (SI) circular ($80-mm diameter) mold that exhibited high thermal insulation. This mold was used to enhance the exothermicity of the reaction and to ensure an anisothermal curing. 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. Type K microthermocouples (110-mm diameter) were used with a Pico V R TC-08 controlled device to assess the temperature in the resin and the oven throughout the curing cycles.

The FBGs used in this work are commercially available apodized 5-mm length Bragg gratings from Welltech instrument Company Limited (Hong Kong), written on bare (uncoated) SMF-28 V R Corning V R optical fibers. Before embedment, each FBG was stabilized by heat treatment at 180 C for 3 h. The FBG were simply placed in each mold at the symmetry center of the specimens (as the thermocouples), in the middle of the resin thickness. Bragg wavelength shifts were measured with a Micron Optics si425 optical sensing interrogator (wavelength resolution 2.5 pm, accuracy 1 pm, and repeatability within 0.5 pm).

Cure Modeling

To study the signals provided by the FBG sensors, a thermokinetic model is used to assess the degree of conversion and the glass transition temperature against time. This section describes briefly this model that has been implemented for each experiment. This approach was detailed and experimentally validated in a previous paper [20] by comparison to dielectric analysis and Fresnel refractometer measurements. This model solves the heattransfer equation including a heat source term resulting from the polymerization, and calculated by a thermokinetic model (Kamal-Sourour modified to take into account diffusion mechanisms below the glass transition temperature). The variation of glass transition temperature is described by a Di-Benedetto law. All the coefficients used in this model come from the literature (see [20][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][23][24][25] for details) as well as the different thermal properties of RTM6.

Two different resolution schemes are used. For thin resin samples (SM and RM experiments), it was demonstrated that an assumption of uniform temperature in the thickness is relevant. A semianalytical resolution based on a finite differences integration scheme is thus possible, which can easily be coupled to an optimization algorithm. For thick resin samples (SI experiment), significant thermal gradients require a 2D finite element resolution (the commercial software COMSOL Multiphysics V R was used).

Whatever the resolution method, the model must be adjusted by identifying the heat exchange coefficient. An optimization allows to find the value of the heat exchange coefficient between the resin and its environment that minimizes the gap between the computed and measured resin temperature (the temperature in the oven and in the resin sample are recorded during each experiment). Excellent results are obtained for isothermal and anisothermal curing, as shown in section Strain Induced during the Curing, where the calculated resin temperature perfectly fits the experimental one. In some cases, the identification requires also to evaluate the initial degree of conversion of the resin, when setting it to zero provides time shifted results. All these developments are discussed in [20].

The last additional information that is needed in this study is the chemical strain e chem of the resin. This data are provided by Aduriz et al. [24] who performed shrinkage experiments of RTM6 at 150 C, and report a linear relationship between the specific volume and the degree of conversion.

RESULTS

This section is devoted to the presentation of the main results. First, the whole curing cycles are presented and global residual strains are analyzed. Then, the cureinduced strains deduced from the FBGs are systematically and quantitatively investigated along with the degree of conversion of the resin, for the different experiments.

Curing Cycles

Several curing experiences were conducted in the different molds cited above (RM, SM, and SI). In each case, a FBG sensor together with a thermocouple was inserted into the mold to monitor the curing reaction. The thermocouple was also used in the processing of the FBG signal to discriminate strain effects from thermal effects according to Eq. 3. Figure 2 presents the temperature change with time for the three curing experiments, whereas Fig. 3 depicts the corresponding Bragg wavelength shifts.

The thermal cycle consists of an increase in temperature from 80 C at a heating rate of 6 C min 21 to limit the temperature gradients (heating phase), followed by a cure plateau where the temperature setting stands at 140 C for 200 min for the SM mold, and 150 C for 130 min for the SI and RM molds. The third period is a simple cool down by natural convection to ambient temperature (cooling phase). Note that for the SM experiment the convection was enforced.

During the cure plateau, the behavior depends on the type of mold. For the SI mold, it is observed in Fig. 2 a huge thermal overshoot which corresponds to the exothermal manifestation of chemical bonding: the temperature reaches more than 240 C as a result of a thick resin sample and an insulating mold. For the RM mold, a small exothermal effect is seen. For the SM mold, the exothermal effect is nearly not detectable because the temperature setting was lowered by 10 C.

A first rough analysis of the Bragg wavelength variations during the experiments (see Fig. 3) allows us to make the following comments.

First, a small wavelength increase occurs, corresponding to the temperature increase of the liquid state resin. Considering the SI experiment, after the temperature overshoot, a large decrease in wavelength is observed at the same time as the decrease in temperature after the exothermal peak. For both metallic molds experiments, no significant wavelength variation occurs during the cure plateau. At the end of the cure plateau, a large decrease in wavelength is observed, corresponding to the thermal contraction of the cured resin in a solid state. Note that in the RM mold the cooling phase shown in Fig. 3 is separated in two parts. This is explained by a debonding of the cured resin from the metallic plate, and also illustrated in Fig. 4. For this figure, the mechanical strain experienced by the FBG sensors immerged in the resin during the cure cycles was calculated from Eq. 3 using the measured Bragg wavelength variations together with their associated temperature variations measured in each mold.

The thermal sensitivities were estimated from temperature and wavelength shift during the heating phase of each experiment (see Fig. 2 and3), under the assumption that FBGs are free of mechanical strain in a liquid surrounding material. Quite usual values were found, respectively, 11.25, 11.4, and 11.45 pm K 21 for the FBG located into the RM, the SM, and the SI molds. In Fig. 4, it is verified that all curves are superimposed during the heating phase and an almost zero strain is measured, due to the fluid state of the resin whatever the experiment.

The large exothermal peak is clearly emphasized for the SI experiment. Before the temperature reaches its maximum (about 225 C), the strain begins to decrease as the temperature increases. This shows that stress transfer has occurred, and at that time the measured strains originate from a chemical shrinkage. After the exothermal peak, an inflexion point corresponding to the glass transition temperature T g SI is highlighted, and will be discussed further. The cure plateaus are shown by the circles in Fig. 4. The strains developed during this phase seem particularly small. During the cooling phase, the behavior depends on the considered mold. It is supposed here that after the cure plateau, a certain degree of conversion is obtained in reference to the curing temperature reached, and therefore that the material undergoes only contraction from thermal origin. This will be confirmed further.

For the RM mold, the cool down is divided into three steps: (i) the strain first decreases linearly with temperature, leading to what is usually considered in the literature as an apparent CTE of the twin-material of 30.9 me K 21 . In this case, the thermal contraction of the resin is hindered by the lower thermal contraction of the metallic mold. (ii) A sudden jump of 3000 me occurs because the resin debonds from the RM mold. (iii) From the jump to the end of the cooling, the slope of the line in Fig. 3 (apparent CTE of the cured resin) is found to be 58.1 me K 21 , as it is considered that the mold does not constrain the cured resin.

For the SM mold, a careful observation of Fig. 4 reveals a nonlinear variation of the mechanical strain after reticulation. At the beginning of the cool down, the temperature is above the glass transition temperature T g of the reticulated resin. Around 110-120 C, the temperature crosses down T g SM and the resin becomes glassy, resulting in a change in the apparent CTE. Below 100 C when the strain varies linearly with the temperature, the apparent CTE of the cured resin can be estimated by the slope of the straight line portion of the curve at 61.7 me K 21 .

For the SI mold allowing totally free deformations of the resin during the whole curing cycle, from the end of the cure plateau at 150 C to the cooling end, the total strain shifts also linearly with temperature. In this case, the apparent CTE of the cured resin is estimated at 56.6 me K 21 , a value of the same order as the other values found.

Strain Induced During the Curing

In this section, are detailed the measured strains from the FBG signal together with the degree of conversion a of the resin resulting from the aforementioned thermokinetic cure model. For the three experiments, the models were calibrated as previously explained [20]. The results fit well the measured temperature cycle, showing that the models are able to physically represent both isothermal as anisothermal curing. It allows to give an estimation of the degree of conversion and also of the glass transition temperature T g . The chemical strain e chem is deduced from the degree of conversion a using the change in the specific volume with a provided by the result of an experiment at 150 C reported in [24].

Figure 5 presents the change in temperature, mechanical strain, and degree of conversion with time for the heating phase and the cure plateau, in the SM mold. During the heating phase, the small change in strain is due to the difference in the positioning of both sensors, and the presence of temperature gradients. The strain remains almost zero while the degree of conversion increases up to aðt ST Þ % 0.68 (ST stands for stress transfer), see Fig. 5. From t ST , temperature remains constant but the strain (only due to chemical contraction) decreases up to 2175 me while a reaches it maximum value, about 0.846. Two main results are highlighted: the stress transfer to the fiber occurs for a fairly advanced level of polymerization, and the value of e chem experienced by the FBG is well below that expected. In fact, it is deduced from [24] and the estimated degree of conversion, a e chem value of approximately 23110 me for a ranging between 0.68 and 0.84, significantly greater than that experienced by the FBG. Finally, the glass transition temperature given by the model is about 125 C at the end of the reaction. This value is very similar to the one estimated by the change in the slope of the cooling curve in Fig. 4 (T g SM , section Curing Cycles). Let us now turn to analyze the curing cycle of the resin inside the RM mold. Figure 6 illustrates the change in temperature, mechanical strain, and degree of conversion with time for the heating phase and the cure plateau. In the present experiment, the initial degree of conversion has been identified to 6.5% for a perfect match between experimental and calculated temperatures. The main tendencies are the same as previously observed for the SM mold: during the heating phase and a part of the cure plateau, the strain remains almost equal to zero while the degree of conversion increases up to aðt ST Þ, estimated in this experience at about 0.63. From t 5 t ST to t 5 100 min, the strain evolves from 0 to 2140 me while a evolves from 0.63 to 0.884. The stress transfer to the fiber occurs for a similar degree of conversion as in SM experiment, and the chemical strain measured by the FBG is also well below that expected, estimated at 24770 me from [24] for this range of degree of conversion. Finally, the model shows that T g is close to the curing temperature, this is why we do not detect any slope change at the beginning of the cooling curve in the RM mold.

Finally, the curing inside the SI mold is presented. The initial degree of conversion was also set to 6.5% for a perfect match between experimental and calculated temperatures [20]. Figure 7 shows the change in temperature, mechanical strain, and degree of conversion with time for the heating phase, temperature overshoot, and cure plateau that occurred in the SI mold. The strain remains zero until t ST and aðt ST Þ % 0.66. From t ST to t Tmax corresponding to the maximum temperature of the exothermal peak, where aðt Tmax Þ % 0.95, a decrease of the strain of about 2300 me is measured, but the increase in the temperature (about 25 C) implies an underestimation of the measurement of this chemical strain by the FBG. The stress transfer to the fiber occurs for a similar degree of conversion as previously and the chemical strain measured by the FBG is also well below that expected, 26100 me, estimated from [24] for this range of degree of conversion. From t Tmax to t 5 80 min, the temperature decreases strongly and results in a large decrease in strain (about 26360 me) mainly conducted by the thermal strain because a changes only from 0.95 to 1. During this cooling, the temperature of the resin crosses the glass transition temperature T g SI , and an inflexion point is seen in Fig. 6 on the strain curve as well as in Fig. 3. At the end of this last experiment, the resin is fully cured and vitrified.

DISCUSSION

Quasistatic Momentum Balance

Let us consider an infinite unidirectional medium made of three components: a mold, a fiber, and a resin sample. The quantities related to each component are denoted by a superscript, respectively, m, f, and r. The different cross-sections are, respectively, S m , S f , and S r .

The total strain rate of the medium is assumed to be uniform in each material and noted _ e. Thus, no warpage occurs. The total strain rate is decomposed into a mechanical part, a thermal part for each material, and a supplementary chemical part for the resin, Eq. 5.

The momentum balance inside the medium implies Eq. 6. The three materials are supposed to behave as linear elastic bodies, the Young's modulus E r of the resin evolves with time because of reticulation and temperature changes, whereas the ones of the fiber, E f , and the mold, E m , are supposed to be constant (Eq. 7). For the sake of simplicity, the viscous behavior of the resin above its glass transition is not accounted for what is obviously a strong assumption.

S m _ r m 1S f _ r f 1S r _ r r 50 ( 6) _ r m 5E m _ e m mech _ r f 5E f _ e f mech _ r r 5E r _ e r mech 1 _ E r e r mech 8 > > < > > : (7)
Using Eq. 7, the integration of Eq. 6 between two arbitrary instants t 0 and t 1 leads to Eq. 8.

S m E m e m mech t 1 ð Þ2e m mech t 0 ð Þ À Á 1S f E f e f mech t 1 ð Þ2e f mech t 0 ð Þ 1S r E r e r mech  à t 1 t 0 50 (8)
If t 0 is chosen as the onset of stress transfer from resin to fiber (which can be the gel time, as reported in [11,12]), the mechanical strain of the fiber is then zero up to t 0 .

Using Eq. 5 and assuming a constant coefficient of thermal expansion for the mold and the fiber, Eq. 8 can be rewritten as Eq. 9.

S m E m 1S f E f À Á e t 1 2e t 0 ð Þ2 S m E m CTE m 1S f E f CTE f À Á T 1 2T 0 ð Þ
1S r E r;t 1 e r;t 1 mech 2E r;t 0 e r;t 0 mech À Á 50

Using Eq. 5 and taking the temperature T 0 , at time t 0 , as the reference temperature for the thermal dilatation, Eq. 9 gives Eq. 10, which shows the complex interaction involved in the present case between the different materials.

S m E m 1S f E f À Á e t 1 2e t 0 ð Þ2 S m E m CTE m 1S f E f CTE f À Á T 1 2T 0 ð Þ 1S r E r;t 1 e t 1 2 ð t 1 t 0 CTE r T;a ð Þ dT dt dt2e r;t 1 chem 2E r;t 0 e t 0 2e r;t 0 chem À Á 50 (10)
Obviously, a key point arises from the fact that both the CTE and the Young's modulus of the resin are functions of temperature and degree of conversion. Two particular cases are discussed further: cooling once the reticulation has ended (SI experiment) and isothermal cure (RM and SM experiments).

Coefficient of Thermal Expansion of RTM6 at the End of the Experiments

During the cooling to room temperature, and as discussed previously in section Curing Cycling, the effect of the mold is negligible in SM and SI experiments, as well as in RM experiment once debonding has occurred. In addition, the modulus of the resin for T significantly lower than T g can be assumed to be almost constant as it only depends on temperature once the reticulation has ended (from 2.76 GPa at room temperature to 1.9 GPa at 150 C [26]). Finally, at the end of the experiments, S f E f can be neglected compared to S r E r (and not only S f << S r ). In this particular case, the strain measured by the FBG is directly related to the thermal contraction of the resin (Eq. 11).

e t 1 2e t 0 5

ð t 1 t 0 CTE r T; a end ð Þ dT dt dt (11)
Hobbiebrunken et al. [26] reported a linear dependency of the CTE of RTM6 to temperature for a fully cured resin. Assuming this linearity, the evolution of the resin CTE with respect to temperature is identified using the FBG data at the end of the SI experiment. The time t 0 is taken at the end of the experiment. Figure 8 shows that this assumption allows to perfectly fit the experimental data from the FBG, compared to the assumption of a constant CTE. The linear evolution with temperature of the CTE of fully cured RTM6, identified from the SI experiment data, is compared to the experimental data reported by Hobbiebrunken et al. [26], Fig. 9. Despite a few overestimated slope of the curve, our results are very close to those of reference 26], with an error lower than 5% between 25 and 150 C.

The same identification was carried out for the SM and RM experiment, showing a variation of the CTE of the resin as a function of the degree of conversion of the resin. It is demonstrated that in-situ measurement using FBG sensor can be processed to assess the CTE of vitrified thermosets for various degree of cure.

Isothermal Chemical Shrinkage Monitoring

In an isothermal case, Eq. 10 can be reduced to Eq. 12.

S m E m 1S f E f À Á e t 1 2e t 0 ð Þ1S r ðE r;t 1 e t 1 2e r;t 1 chem À Á 2E r;t 0 e t 0 2e r;t 0 chem À Á Þ50

In Eq. 12, the strain variation e t 1 2e t 0 is given by the FBG and e r;t 1 chem 2e r;t 0 chem can be evaluated by the model using shrinkage data from the literature [24]. The former is one order of magnitude lower than the latter in both RM and SM experiments, see section Strain Induced during the Curing. This difference can be attributed on the one hand to the mold stiffness, and on the other hand, to the resin modulus variation, which remains the main unknown. In both experiments, the cure-induced strain of the fiber is almost the same, as well as the actual chemical shrinkage, despite a significant difference in the mold contribution (S m E m =S r is 60 times higher in the RM mold). This assertion results from the underlying assumptions leading to Eq. 10, that does not take into account the viscous relaxation in the resin. The mold does not interact with the resin in the vicinity of the fiber. As shown previously in section Coefficient of Thermal Expansion of RTM6 at the End of the Experiments, the assumption of Eq. 5, uniform strain in the fiber, the mold, and the resin, was valid for the fully cured resin. It is no longer valid here because the resin is not a solid, changing from a viscous liquid to a viscoelastic gel.

SUMMARY AND CONCLUSIONS

Curing experiments of RTM6 resin were conducted in several types of molds, each one instrumented with FBG and thermocouple, to assess the cure-induced strains. A thermokinetic model allowed to evaluate the degree of conversion of the resin throughout the cure. At the onset of stress transfer to the optical fiber, values of the degree of conversion between 0.63 and 0.68 are reported. This is slightly higher than the degree of conversion at gel point at 150 C reported in [24], that is, 0.59, despite very different experimental approaches. Nevertheless, this seems logical as the resin must be rigid enough to deform the optical fiber, what is not the case just after the gel point where the Young's modulus is close to zero. Whatever the curing experiment, the FBG deforms under the chemical shrinkage of the resin. Nevertheless, the amplitude of the FBG strain is one order of magnitude lower as the expected chemical strain. This can be attributed to the viscoelastic behavior of the partially cured resin above its glass transition, allowing stress relaxation and limiting the interactions with the mold as well as with the optical fiber. The FBG signal is thus difficult to link to the actual chemical shrinkage. However, once the resin is cured, it is demonstrated that the FBG strain provides directly the CTE of the resin, taking into account its linear change with temperature. The effect of the degree of conversion on the CTE value was also emphasized from the different experiments.

  à r =e à z is independent of the temperature change because each strain component is proportional to DT. It is interesting to analyze the change in B by varying the elastic modulus of the resin and the external radius of the resin cylinder. The other parameters m r and CTE r are assumed to be constant and set, respectively, to 0.4 and 60 3 10 26 K 21 (see section Curing Cycles) corresponding to the particular case of the cooling of a vitrified resin at the end of the curing process. Results are shown in Fig. 1a for E r in the range [0.1-10] GPa and R e in the range [62.5-1000] mm.

FIG. 1 .

 1 FIG. 1. (a) Biaxiality ratio and (b) axial strain error committed by assuming a constant B ratio, function of R e for various resin elastic modulus. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FIG. 3 .

 3 FIG. 3. Bragg wavelength shift throughout the curing cycle of the resin for the three involved molds. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FIG. 4 .

 4 FIG. 4. Mechanical strain variations versus temperature throughout the curing cycle of the resin for the three involved molds. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary. com.] FIG. 5. Temperature, mechanical strain, and degree of conversion versus time for the temperature increase and the cure plateau at 138 C of the RTM6 curing cycle in the SM mold. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FIG. 6 .

 6 FIG. 6. Temperature, mechanical strain, and degree of conversion versus time for the temperature increase and the cure plateau at 148 C of the RTM6 curing cycle in the RM mold. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

_ e5 _ e m 5_ e m mech 1_ e m therm _ e f 5_ e f mech 1_ e f therm _ e r 5_ e r mech 1_ e r therm 1_ e r chem 8

 8

FIG. 8 .

 8 FIG. 8. Comparison of the thermal strain of fully cured RTM6 assuming a linear evolution of the CTE with temperature and a constant CTE value to the experimental data from SI experiment. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary. com.]

FIG. 9 .

 9 FIG. 9. CTE of RTM6 as a function of the temperature evaluated at the end of the SI, RM, and SM experiments compared to the experimental data reported by Hobbiebrunken et al. [26]. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary. com.]

TABLE 1 .

 1 Material parameters of optical fiber SMF-28 V R[START_REF] Tao | [END_REF].

		Elastic modulus	Poisson's	Coefficient of thermal
	Radius [R (mm)]	[E (Gpa)]	ratio (m)	expansion [CTE (K 21 )]
	62.5	70	0.17	0.55 3 10 26 K 21

ACKNOWLEDGMENT

The authors acknowledge C. Vergnes for his major contribution to the experimental work.


R Gantois
email: renaud.gantois@mines-albi.fr

A Cantarel
email: barthur.cantarel@iut-tarbes.fr

G Dusserre
email: cgilles.dusserre@mines-albi.frd

J.-N Félices

F Schmidt

Mold Filling Simulation of Resin Transfer Molding Combining BEM and Level Set Method

Keywords: Numerical simulation, Boundary Element Method, Level Set, Liquid Composite Molding, Darcy's law

Liquid Composite Molding (LCM) is a popular manufacturing process used in many industries. In Resin Transfer Molding (RTM), the liquid resin flows through the fibrous preform placed in a mold. Numerical simulation of the filling stage is a useful tool in mold design. In this paper the implemented method is based on coupling a Boundary Element Method (BEM) with a Level Set tracking. The present contribution is a two-dimensional approach, decoupled from kinetics, thermal analysis and reinforcement deformation occurring during the flow. Applications are presented and tested, including a flow close to industrial conditions.

Introduction

Liquid Composite Molding (LCM) is used in a wide range of industries (e.g. automotive, aeronautics, spatial, marine, defense, goods, etc.). Since the apparition of Resin Transfer Molding (RTM), a great number of variants has been developed in order to improve the forming conditions, lower production costs and reduce processing time. Among them can be mentioned Vaccuum Assisted Resin Transfer Molding (VARTM), Structural Reaction Injection Molding (SRIM) and Liquid Resin Infusion (LRI). The main feature in common is a liquid resin flow through the fibrous preform. A typical RTM setup [1] is depicted on Figure 1. The process can be divided into three major steps a) the dry plies are cut and stacked, b) the reinforcement is preformed and put in place in the mold, c) the mold is filled using the liquid resin, d) the part is cured and demolded.

In industry, a condition to meet the mechanical specifications is the entire impregnation of the preform. For years, engineers have worked using the trial and error method in order to optimize the filling conditions. An acceptable design is often reached after iterations because it involves the control of a large numbers of processing parameters. Considering the tool costs, it is clear that numerical simulation of flow is of a great interest, in particular for complex parts. Motivated by this need, numerical modeling has become an extensive field of investigation for the last three decade. Several models based on continuum mechanics and heat transfer has been studied at different scales. Simulations at fibre or tow scale are generally conducted using Stokes' flow.

Due to the complexity of the fibrous network, Representative Elementary Volumes are usually considered. Simulations at part scale are frequently adopted in industry considering an average flow modeled using Darcy's law. As the resin flows, a tracking technique is employed to follow the moving front in the mold. It can be performed using Lagrangian or Eulerian approaches. The first approach offers an interesting framework because the moving mesh explicitly describes the boundary [2][3][4]. Previous works [3] showed that the technique is highly suitable in combination with a Boundary Element Method. However, it becomes difficult to implement as soon as topological changes appears, which is frequently encountered in mold filling simulation. Eulerian approaches are clearly more consistent for that task. In such a framework, an auxiliary field is used to capture the front. In Volume Of Fluid (VOF) methods [5,6], it is defined at each cells as the fluid fraction. Front is implicitly located through connecting the partially filled cells. Some enhancements aim to approximate the front shape using reconstruction algorithms (e.g. Simple Line Interface Calculation (SLIC) or Piecewise Linear Interface Construction (PLIC)). As the field is discontinuous, a major drawback is that it requires a large number of cells to reach an acceptable level of accuracy, which tends to dramatically increase the computational time. In Level Set methods [7,8] the auxiliary field is a continuous function, usually chosen as the signed distance to the front. Thanks to the regularity of the field, an accurate shape of the front can be computed. The present paper focus on a technique combining a Boundary Element Method (BEM) [2,3,8,9] with a level set technique, following the approach developed in reference [8]. Level set theory indicates that the front motion is governed by the boundary velocity field. As a result the use of BEM together with level set techniques is straightforward to compute the filling pattern. Besides, the versatility of boundary mesh makes easier the remeshing procedure comparing to Finite Element Methods. The first part of this paper considers the governing equations. The second part describes the implemented model. The last part covers some applications.

Governing equations

Modeling Assumptions. Impregnation is modeled as a Darcy's flow. Fibrous preform is regarded as a porous medium. As a first approach, we focus on the flow front numerical modeling. Hydro-mechanical couplings are neglected, so that the medium is assumed to be not deformable and of a constant permeability. It should be noted that BEM technique allows coupling preform deformation and resin flowing but this is not addressed in this paper. The flow is considered as isothermal and resin viscosity is assumed to be constant, neglecting chemical reactions effects. Again, coupling heat balance equation including reactive source term remains possible but has not been investigated. Consequently, governing equations reduce to Darcy's law and incompressibity equation

[ ]      >= < ∇ ∇ - >= < (2) 0 . (1) ' v p K v ! ! ! ! µ where > < v ! is
     >= < ∇ ∇ - >= < (4) 0 . (3) ' v p K v e ! ! ! ! µ where e K is given as 2 1 K K K e = , 1
K and 2 K being the principal permeabilities. The isotropic equivalent transformation [3,10] is given by (5) Outline. The algorithm is implemented using Matlab environment. At the beginning, preprocessing imports a standard mesh file, performs the isotropic equivalent transformation, assigns material and processing data, and locates the injection ports. Next, the filling routine based on a Level Set formulation advances the front. At the end, a post-processing program is run in order to visualize the filling pattern once the real domain is recovered. The filling stage is divided into a finite number of quasi steady states. At each step time, governing equations are solved using a constant Boundary Element Method. The front is advanced by feeding a level set solver with the BEM-computed velocities. The process is repeated until the mold is completely filled. (10) Rewritten in a matrix form using H and G , Eq. 10 can be transformed into

1 0 0 1 2 1               =       y x K K K y

Boundary Element

∫ ∑ ∫ ∑ Γ = Γ = Γ = Γ + j d p q d q p p c N j j N j j i i
ij 1 1 G q H p N j j ij N j j ∑ ∑ = = = (11)
where

∫ Γ Γ + = j d q H ij * ij
F AX = (12)
Where X is the vector of N unknowns, A is a 2 N matrix obtained by reordering H and G . F is the known vector computed from the boundary conditions, H and G matrix.

Level Set method. The moving front is captured using the signed-distance function φ , defined so that the zero level set corresponds to the interface [7,13]

{ } 0 ) , , ( / ) , ( ) ( 2 = ∈ = Γ t y x IR y x t φ (13)
For each point under consideration, distance is signed negative if located in the impregnated area and positive otherwise. Two types of mesh are utilized. The first one is a background grid. It is generated by a standard triangle element mesher using the CAD of the cavity. The other mesh is a computational boundary mesh. It bounds the fluid phase at each step time. Level set equations govern the evolution of the signed-distance function. To combine with BEM, we use a formulation involving a propagating interface with a velocity in its normal direction. It is given as follows

     = = = ∇ + ∂ ∂ (15) ) 0 , , ( (14)
0 0 φ φ φ φ t y x F t !
where F is the extended normal velocity at time t, built in order to coincide with the real velocity of the front (it is extrapolated outside without physical meaning). The filling stage is initialized by choosing 0 φ as the signed-distance to the inlet gates. Next, values for each point are updated using an Euler scheme

(16) ) , , ( ) , , ( t F t y x t t y x ∆ ∇ - = ∆ + φ φ φ ! Where t
∆ is the time step, adjusted to match its upper limit according to Courant-Friedrich-Levy (CFL) conditions. F and φ ∇ ! are computed at time t.

For the sake of efficiency, φ can be updated on few nodes around the front in a "narrow band" [7] instead of the entire grid. By the next, front is updated using a boundary element mesher developed during this work. It is based on interpolating the zero level set on the grid through a Marching Triangles algorithm.

It should be noted that the method directly handles topological changes, involving merging or dividing fronts. However it does not ensure that the resin remains inside the mold. Resin-mold contact is implemented using the fixed level set describing the mold walls.

Applications

Isotropic radial injection. The numerical method, was first validated by considering the academic case of an isotropic radial injection. The mold chosen for that study is a 0.25 m 2 square plate. The resin enters the part from the center and flows through the preform describing circular patterns. Four vents located on the corners ensure the complete filling. Material and processing parameters are given in Table 1. Table 1 Material and processing parameters

K [m 2 ] µ [Pa.s] ρ [kg.m -3 ] ε [-] 0 p [Pa] f p [Pa] 0 r [m]
1.10 -9 0.1 1150 0.5 2.10 5 1.10 5 2.5.10 -3

Neglecting gravity, the front motion in the free flow stage is given by the following analytical solution [3,10] (17)

) ( 4 1 1 ln 2 . 2 0 0 0 2 0 r t p p K r r r r f f f εµ - = +         -                
where f r and 0 r are the radii of the moving front and the inlet gate, 0 p and f p the inlet and outlet pressures. The mold is meshed using 2500 triangle elements. CPU time is around 20s on a 2.26 GHz / 1.93 Go of RAM laptop. The predicted filling time is 98s. Figure 2 shows the comparison between the implemented the analytical results at time 10.7s and 63.3s. The agreement is fair, with an error of 2.1%. Anisotropic injection involving complex shapes. The model was also compared with a more realistic case. A test was performed consisting of impregnating a 1x1 rib knit fabric made of glass fibers [14] using a canola oil of known viscosity. The preform was placed under a transparent flexible bag and filled using two inlet gates under vaccum conditions (infusion). Front is pulled using two vaccum lines on the upper and lower sides. The dimensions of the mold are given in Fig. 4. Some inserts were placed to assess the contact algorithm. Material and processing parameters are given in Table 2. CPU time is 200s for a model involving 3693 elements. Data was compared at 0.5s, 3s, 5s and 7s. Figure 6 indicates that agreement is fair at any time. In particular, a front merging is accurately predicted.

Conclusion

The presented method is based on combining BEM with Level Set method. Velocity computation is reduced to boundary, which tends to reduce the computational time. An important feature of the method is that the front shape is accurately computed which is more difficult to achieve in standard VOF techniques. Comparing with a Lagrangian approach, the Level Set front tracking facilitates the treatment of topological changes, typically appearing when fronts merge into each other.

Possible extensions of this work may include other physics occurring during the flow (e.g. hydromechanical couplings, kinetics). It should be noted that the implemented model includes a gravity term, and can be easily modified to take into account other body forces (inertia effect due to rotation for example). The technique remains similar in three dimensions, considering the resin flow as surfaces propagating.

  Figure 1 RTM process

(

  the physical coordinates.Finally, the resin velocity v ! is given by dividing the average velocity by the reinforcement porosity ε, which turns the system of Eqs. 3 and 4 into Neumann condition)where p (resp. q ) is the prescribed value of pressure (rep. normal derivative). In this paper, 0 = q on the mold walls (non penetration condition), p being the prescribed value of pressure at inlet and outlet.

2 N

 2 matrix. Finally Eq. 11 is reordered to take the form of a linear system as

Figure 2 Figure 3

 23 Figure 2 Filling patterns in isotropic radial injection: comparison with analytical solutionThe same simulation was conducted by activating the gravity contribution. A comparison between the implemented method and a FEM/VOF-based code (PAM-RTM™) is given Figure3. The same mesh was used for both simulations. As expected, the flow is slightly deported in the lower part of the plate. The predicted filling time is 101s for the implemented model and 89s for PAM-RTM™. The difference is less than 14%, which shows a good qualitative adequacy.

Figure 4 1 K [m 2 ] 2 KFigure 4

 4124 Figure 4 Dimensions of the mold

  and g denote respectively the resin specific mass and gravity. Combining the previous equations and transforming the coordinates in the isotropic equivalent domain e

	p	'	=	p	+	gz ρ
	where, ρ Ω of permeability

the macroscopic velocity, ] [ K the permeability tensor, µ the liquid resin viscosity and ' p the acting pressure. That pressure includes a gravity term. It is given by e K , we obtain

  Method. For clarity, the subscript e referring to the isotropic equivalent domain is omitted by the next. Laplace's equation is multiplied by the Green function and θ is the internal angle of the corner in radians. For a two-

	Where i p is the value of the pressure at a point i M located on the boundary, * q the pressure
	gradient associated with
	dimensional domain,	* p and * q are given as follows
	p	*	=	2 1 π	  	1 r	  	and	q	*	=	2 1 π	2 . n r r -	(9)
	where	r	=	M	i	P	, i M the point of application of the Dirac delta function and P the point under
	consideration. Meshing the boundary into N constant boundary elements and applying Eq. 8 leads
	to
											*								*
						1										1
					i c	i p	+	Γ ∫	pq	*	Γ d	=	∫ Γ	p	*	q	Γ d	with	i c	=	π θ 2	(8)

* p .

Integrated twice by parts over the calculation domain and using Green's theorem leads to Somigliana's equation

[3,11,12]

* p

Acknowledgments

This work was carried out in the framework of FUSCOMP project coordinated by DAHER-Socata. The authors gratefully acknowledge B. Cosson and S. Leroux for their contributions during this work.