H Chraïbi

M Prat
email: prat@imft.fr

O Chapuis

Influence of contact angle on slow evaporation in 2D porous media

Keywords:

We study numerically the influence of contact angle on slow evaporation in 2D model porous media. For sufficiently low contact angles, the drying pattern is fractal and can be predicted by a simple model combining the invasion percolation model with the computation of the diffusive transport in the gas phase. The overall drying time is minimum in this regime and is independent of contact angle over a large range of contact angles up to the beginning of a transition zone. As the contact angle increases in the transition region, cooperative smoothing mechanisms of the interface become important and the width of the liquid gas interface fingers that form during the evaporation process increases. The mean overall drying time increases in the transition region up to an upper bound which is reached at a critical contact angle θ c . The increase in the drying time in the transition region is explained in relation with the diffusional screening phenomenon associated with the Laplace equation governing the vapor transport in the gas phase. Above θ c the drying pattern is characterized by a flat traveling front and the mean overall drying time becomes independent of the contact angle.

Drying time fluctuations are studied and are found to be important below θ c , i.e. when the pattern is fractal. The fluctuations are of the same order of magnitude regardless of the value of contact angle in this range. The fluctuations are found to die out abruptly at θ c as the liquid gas interface becomes a flat front.

I. INTRODUCTION

Prediction of evaporation from porous media is of interest for many environmental and industrial applications, such as the water exchange between soil and atmosphere, drying of many products or the recovery of volatile hydrocarbons from underground oil reservoirs, to name only a few. In these applications, the liquid is generally wetting and drying can then be regarded as a drainage process, where the evaporating wetting liquid is replaced by a nonwetting gas. This case has been the subject of many studies in recent years within the framework of pore-network models and invasion percolation theory, e.g. [1], [START_REF] Prat | [END_REF] and references therein. However, there are also important applications in which evaporation occurs in a hydrophobic porous medium. For example, buildings stones can be rendered hydrophobic in order to limit the salt weathering hazards due to the salt crystal formations resulting from evaporation, [START_REF] Goudie | Salt weathering hazards[END_REF]. More generally, rendering the porous medium hydrophobic is often considered for reducing the evaporation from porous surfaces or soils. PEFCs (polymer electrolyte fuel cells) is another example in which evaporation, as well as other two-phase flow processes, takes place in partially hydrophobic (teflonized) porous (fibrous) layers. A proper understanding and modeling of two phase flows, including evaporation, in teflonized systems is needed in relation with the problem of water management of PEFCs, [START_REF] Wang | Handbook of Fuel Cells[END_REF]. Yet, it is often considered that the wettability condition may change in time and this may be the cause of performance degradation. Hence it is important to study the impact of wettability change on evaporation. Here we explore the case where the wettability condition, i.e. the contact angle, is uniform throughout the porous medium. The study of situations implying spatial variations of contact angle is left for future works. As shown in [START_REF] Cieplak | [END_REF] for the simpler process of mechanical quasi-static displacement, the invasion pattern changes from a fractal pattern to a compact one as the contact angle rises. The compact pattern is obtained for contact angles greater than a critical angle θ c . We are especially interested in the transition region right below θ c where the invasion pattern changes significantly. We expect a similar pattern transition for the slow evaporation process considered in this paper.

Experimentally, it is not obvious to study in details the effect of contact angle on drying pattern in the transition zone since it is no easy to impose at will a given value of θ. As shown in Figure 1, a partial illustration on the effect of contact angle can, however, be obtained from a simple drying experiment of a model porous made of a monolayer of 1mm glass beads randomly distributed and sandwiched between two glass plates. Without additional treatment, this leads to a contact angle θ ∼ 30-40° (measured in the liquid phase). To obtain a much higher contact angle, the glass beads and the glass plates are rendered hydrophobic by silanization, a process which leads to a contact angle of the order of 105°-107°, [6]. The model is initially fully saturated by pure water. Vapor escapes through the four open lateral edges of model. The 8.8 cm long and 7.5 cm wide model is placed horizontally in a small transparent Plexiglas chamber of controlled temperature (22 ± 1 o C ). The relative humidity in the chamber is stabilized using a LiCl saturated solution (RH =12 %) and the evolution of phase distribution within the porous medium is recorded using a CCD camera set above the chamber. More details on this experiment can be found in [START_REF] Sghaier | Prat submitted to Transp[END_REF]. For θ ∼ 30-40°, one can observe in Figure 1 the highly ramified capillary fingerings and trapped liquid clusters typical of invasion percolation patterns as expected. In contrast, the invasion pattern is much more compact with no trapping for θ ∼105°-107° and resembles the faceted pattern that will be shown bellow for θ=99° (see Figure 7). In [START_REF] Chapuis | [END_REF], we analyzed and compared the two regimes flanking the transition region. In the present effort, we concentrate on the transition region. This is made possible by taking into account in the simulations all the local mechanisms involved in the growth of the liquid-gas interface.

Throughout the text the wettability of the system is characterized by the contact angle θ measured in the dense phase, i.e. liquid water in the present context. This contact angle can be interpreted as the receding contact angle since menisci recedes during invasion of pores as a result of evaporation. As shown in a previous study, [START_REF] Chapuis | [END_REF], the change in wettability from a hydrophilic materials (θ << 90°) to a hydrophobic one (θ >>90°) leads to a major change in terms of invasion pattern and therefore also in terms of evaporation rate. As in [START_REF] Chapuis | [END_REF], we study this change and its impact on drying rates in the quasi-static limit, i.e. when the pressure evolution in the liquid phase is only due to capillary effects. The study is conducted for a simple two dimensional (2D) model porous medium constructed by placing disks of random radii on a square lattice under the assumption of dilute concentration of water vapor in the gas phase (which is acceptable for evaporation at temperatures close to the ambient temperature).

For convenience, the liquid is referred to as "water" throughout the text. It is obvious that the results are general and apply to other volatile liquids (under the condition of slow evaporation).

As shown in previous works, e.g. [START_REF] Chapuis | [END_REF], [9], [10] and references therein, liquid films can be a major transport mechanism in the drying of porous materials. Liquid films in drying refer to "thick" films trapped by capillary effects in corners, grain contacts or surface roughness in the pores (as opposed to the thin films adsorbed by solid surface forces). For the not too small contact angles of interest in this study, it is, however, reasonable to neglect the effect of films, see [START_REF] Chapuis | [END_REF] for more details.

The paper is organized as follows. In Sec.II the model of drying developed to simulate the evaporation process in the two dimensional (2D) model porous medium is presented. The influence of contact angle on drying is discussed in Sec. III through the analysis of evolution of mean overall drying time and drying time statistical fluctuations with the contact angle.

We close in Sec. IV by offering some concluding remarks.

II MODEL

A. Model porous medium, transport mechanisms and external boundary condition

For simplicity, we consider a situation where heat transfer can be neglected. This corresponds typically to evaporation of a not too volatile liquid, e.g. water, at the ambient temperature. As in [START_REF] Cieplak | [END_REF], we use a 2D array of disks with random radii (Figure 2). As shown in [START_REF] Cieplak | [END_REF], a distinguishing advantage of this model is to allow for a full solution of the interface shape. In particular, it is possible to study in detail the influence of contact angle on the invasion pattern resulting from the evaporation process. In this model, the liquid-gas interface consists of circular arcs connecting disks. One can refer to [START_REF] Cieplak | [END_REF] for more details on the advantage of this model compared to the more traditional pore network models made of interconnected channels. In our case, the disks are placed on a square lattice with lattice constant a. There are L disks per row and L rows on the lattice. Disk radii are distributed randomly according to an uniform distribution law in the range [ ]

max min r , r
, see details in Table I regarding the particular systems considered in this paper. The porosity ε in Table I is the fraction of area not covered by disks.

Initially, the pores of the model porous medium are saturated with liquid, which is then allowed to evaporate isothermally. As depicted in Figure 2, the three upper sides of the square domain are sealed, i.e. they represent the impervious surfaces. The fourth lowermost side of the domain is open to air for drying. On this side, we use a standard mass-flux boundary condition of the form

) P P ( h j , v o , v ∞ - = (1)
where j is the evaporation flux density, h is the mass transfer coefficient, o v P , is the vapor partial pressure at the porous medium surface, and ∞ , v P is the vapor partial pressure in the surrounding air. The mass transfer coefficient h is expressed as

( ) δ T M R D h v = ( 2
)
where D is the diffusion coefficient of water vapor, R is the universal gas constant, M v is the vapor molecular weight, T is temperature, and δ is the external transfer length-scale (we took δ = a throughout this paper). Referring to the classical convective drying situation, δ can be seen as the average thickness of the mass boundary layer developing at the porous medium surface.

At the surface of the menisci, the vapor partial pressure is the saturation vapor partial pressure P vs and the difference in partial pressure

) P P ( , v vs ∞
can be regarded as the mass transfer potential driving the evaporation process.

B Fluid interface geometry and interface growth mechanisms

As in [START_REF] Cieplak | [END_REF], the interface between the gas phase and the liquid phase consists of sequence of arcs between pairs of disks. As sketched in Figure 3, each arc intersects both disks at the proper contact angle θ. At the beginning of the process, we assume that the medium is saturated with liquid and the interface is formed by L stable arcs attached all along the first row of disks, see Figure 2. Then as a result of evaporation, the curvature of the arcs increases until one arc becomes instable. This leads to the invasion of adjacent pore as depicted in Figure 3 and new arcs are positioned at the entrance of the new interfacial pores. Then this process is repeated until full invasion of the system. Hence, in this quasi-static limit, it is assumed that the invasion is characterized by a succession of stable configurations, with very rapid motion of the interface between two successive stable configurations.

As discussed in [START_REF] Cieplak | [END_REF], three basic types of local instability may be distinguished:

1) "burst" when the capillary pressure between the two fluids is such that it becomes impossible to find a stable arc between two disks intercepting them with the right angle θ.

2) "touch" the arc connecting two disks intersects during its growth another disk at the wrong θ.

3) "overlap" = menisci coalescence, two neighboring arcs on the interface intersect.

Figure 4 shows the relative importance of various growth mechanisms as a function of contact angle θ in our system to reach breakthrough, i.e. when the invading gas phase reaches the top edge. As can be seen from Figure 4, the interface growth is dominated by bursts for sufficiently small contact angles whereas the overlap mechanism is dominant for sufficiently high contact angles. As discussed in details in [START_REF] Cieplak | [END_REF], this change in the probabilities of each local growth mechanism is responsible for major changes in the invasion pattern, from invasion percolation patterns to a compact pattern with no trapping.

To construct the evaporation algorithm, we begin by considering the simpler case of the quasi-static mechanical displacement of the fluid in place by the invading fluid. Contrary to [START_REF] Cieplak | [END_REF] where the interface growth was driven by small changes in the pressure applied to the system, the procedure is slightly different in our case and can be considered equivalent to imposing a very small flow rate. In particular, only one pore is invaded between two successive stable configurations of the interface. Figure 5 shows the various cases of pore invasion considered in our model. A pore occupied by the liquid phase and adjacent to a pore occupied by the invading gas phase is defined as an interfacial pore. We determine for each interfacial pore the local instability mechanism (burst, touch or overlap) that is encountered first during the growth of the interfacial arcs and denote by R(i) the corresponding curvature radius of the arc(s) associated with the considered interfacial pore, i.e. the curvature radius at the onset of the local instability. The instability radii R(i) are obtained by computing the evolution of the curvature radius of each interfacial arc, as it moves between the two disks it connects. The computation is analogous to that presented in [START_REF] Cieplak | [END_REF]. For this reason, the details, which are presented in [11], are omitted here. We can then define as pore invasion potential

a ) i ( R ) i ( P = (3)
which is inversely proportional to the pressure difference between the two fluids needed for the invasion of the pore (=

) i ( R
γ where γ is the surface tension). Under these circumstances, the invasion algorithm becomes analogous to the invasion percolation algorithm, [12], i.e. the interfacial pore of greatest potential is invaded at each step of the invasion. The fundamental difference lies in the computation of invasion potential (3), which takes into account all possible local instability mechanisms and not only the burst mechanism as in the invasion percolation algorithm. If one adds the rule that the trapped pores, i.e. the liquid pores that become completely surrounded by the invading phase, cannot be invaded, one obtains an algorithm describing the quasi-static displacement of one phase by the other under the condition of a very small flow rate of injected fluid, i.e. when the capillary forces dominate the displacement.

As first proposed in [13], the algorithm describing slow evaporation can be constructed by combining the mechanical quasi-static displacement algorithm described above and the computation of diffusive transport of the water vapor in the gas phase. This will make the process time dependent through the computation of evaporation rate, the (variable) time step being of the order of ρ l V p /F, where ρ l ,V p and F are the liquid density, average volume of one pore and evaporation mass flux at the boundary of liquid cluster to which belongs the considered pore, respectively. More precisely the (quasi-)steady diffusion equation . 0 = ∆ v P is solved in the part of the pore space occupied by the gas phase with the boundary conditions: P v = P vs on the interfacial arcs, and Eq.( 1) at the entrance of pore located along the first row of model porous medium, see [13] for more details regarding the computation of P v thanks to a finite volume type discretization technique, which has been straightforwardly adapted here to the particular case of disks array.

For the sake of completeness, the drying algorithm is summarized here: (1) every liquid cluster present in the network is identified, (2) the interfacial pore of greatest invasion potential is identified for each cluster, the volume of liquid contained at time t in this pore is sc V , (3) the evaporation flux F c at the boundary of each cluster is computed from the finite volume computation of the liquid vapor partial pressure in the gas phase (obtained from the numerical solution of equation .

0 = ∆ v P
), (4) for each cluster, the time t c required to evaporate the amount of liquid contained in the interfacial pore identified in step (2) is computed:

c sc c F V t l ρ =
, (5) the element among the elements selected in step (2) eventually invaded is that corresponding to t cmin = min (t c ), ( 6) the phase distribution within the network is updated, which includes the partial evaporation of liquid contained in the interfacial pores

selected in step (2) with min min ) ( ) ( c c c t F t V t t V - = + l l l l ρ ρ
(where l ρ is the liquid density and l V the volume of liquid contained in the interfacial pore) except for the interfacial pore selected in step (5) which becomes completely saturated by the gas phase. The procedure can be repeated up to full drying or stopped at some intermediate stage.

III RESULTS

As mentioned in the introduction the two regimes limiting the transition zone were studied in [START_REF] Chapuis | [END_REF], where it was shown that the dimensionless average overall drying time of a hydrophilic system was system size dependent approaching exponentially a limit for large size system. However, the dependence with the size is relatively weak. Furthermore, the statistical fluctuations of the drying time were found to be of the same order of magnitude for all system sizes investigated. Therefore, we can concentrate here on small systems without loss of generality and do not explore again the (weak) influence of system size. In passing, it can be noted that the computational time of drying for one realization is significantly greater here due to the consideration of all local instability mechanisms of the interface. This makes more difficult the numerical study since fewer realizations can be computed for a similar computational effort.

A. Influence of contact angle on overall average drying time and drying pattern

The overall drying time t o is the time needed to fully evaporate the liquid contained in the model porous medium. Figure 6 displays the evolution of overall drying time as a function of contact angle for the four systems considered. To obtain the results shown in Figure 6, we have considered water at a temperature of 20°C as working fluid and taken a = 1mm. Unless otherwise mentioned, the results presented in this paper were obtained considering 100 realizations of a 25x25 disk array.

Interestingly, the average overall drying time is independent of contact angle for sufficiently low or high contact angles (as shown by the plateaus on the left hand side and the right hand side of the curves in Figure 6). These plateaus correspond to the regimes studied in [START_REF] Chapuis | [END_REF], i.e. the invasion percolation regime for sufficiently low contact angle when the growth of the interface takes place in the largest local constriction available along the interface of a liquid cluster (this regime corresponds to bursts as dominant local instability mechanisms) or, for sufficiently high contact angle, when the pattern is characterized by a flat traveling front (menisci overlaps are then the dominant local instability mechanism). These two regimes are termed IP (invasion percolation) regime and FF (flat invasion front) regime in the following.

The results shown in Figure 6 for the IP and FF regimes are consistent with the ones obtained in [START_REF] Chapuis | [END_REF] for a pore network model of interconnected channels of rectangular cross-section, i.e.

drying in the hydrophilic model (IP regime) is faster than in the hydrophobic one (FF regime).

The ratio between the lowest and largest average drying times is slightly greater than 0.6 as shown in the inset of Figure 6 for the disk array whereas a ratio slightly higher of 0.7 was obtained with the pore network model in ref. [START_REF] Chapuis | [END_REF] for a network of comparable size.

Interestingly, Figure 6 shows that the transition between the IP and FF drying times occurs over a relatively narrow range of contact angles. As shown in Figures 4 and6, the greater the porosity, the wider the transition zone.

In the analysis of the pressure driven quasi-static invasion process considered in [START_REF] Cieplak | [END_REF], a critical contact angle θ c was introduced to characterize the transition between uniform flooding and fractal invasion patterns. Noting that the contact angle is defined in the invading phase in [START_REF] Cieplak | [END_REF] and not in the invaded phase as in our case, see Figure 2, the results reported in [START_REF] Cieplak | [END_REF] indicates that uniform flooding occurs above θ c (using our definition for the contact angle). As it will be made clear from the consideration of drying patterns and drying time fluctuations, θ c corresponds to the contact angle marking the end of transition region on the right hand side in Figure 6, i.e. the contact angle marking the beginning of the FF regime. Approximate values of θ c deduced from Figure 6 are reported in Table I for four porosities. Consistently with the results of [START_REF] Cieplak | [END_REF], θ c appears to increase monotically with porosity. Whereas θ c increases with porosity, the contact angle marking the beginning of the transition region is found to be independent of porosity and close to 90°, as can be seen from Figure 6.

It is interesting to observe that the simple drying model proposed in [13], which is based on a the IP algorithm for modeling the capillary effects, is sufficient to predict the average overall drying time over the full range of contact angle [0, 90°] , provided that liquid film effects can be neglected, see [9] and references therein for more details on the effect of films. Referring to both Figures 4 and6, it can be seen that the average overall drying time begins to be affected only when the probability of overlap mechanism becomes greater than the burst mechanism probability.

Figure 7 shows examples of drying patterns for system D (ε = 0.7). The same realization of the 25x25 disk array is considered for four different values of the contact angle, namely θ=80° (invasion percolation (IP) regime), θ=94° (transition region), θ=99°(transition region) and θ=120° (flat traveling front = FF regime). The IP regime (θ=80°) is characterized by ramified fractal patterns and the occurrence of many disconnected liquid clusters whereas the FF regime (θ=120°in Figure 7) is characterized by an almost flat traveling invasion front with no trapping. As shown in [START_REF] Cieplak | [END_REF] for the case of the quasi-static fluid invasion and illustrated in Figure 7 for the case of the evaporation problem considered here, the trapping phenomenon progressively disappears as the contact angle increases in the transition region. The growth of the interface is increasingly smoother as the contact angle is increased in the transition region and this is due to the increasingly cooperative nature of invasion as the overlap mechanism becomes the dominant local instability mechanism. As explored in detail in [START_REF] Cieplak | [END_REF], larger and larger segments of the interface move forward coherently as the contact angle increases in the transition region and this leads to the "faceted" interface shown in Figure 7 for θ=99°. For a sufficiently high contact angle (above θ c with the convention for the contact angle used here) , the entire interface advances coherently and this leads to the FF regime.

The fact that the average drying time increases with the contact angle in the transition region can be understood, at least qualitatively, from the evolution of the drying patterns shown in Figure 7 and the results reported in [START_REF] Cieplak | [END_REF]. If one ignores the external mass transfer resistance (~ h -1 from Eq.( 1)), which is independent of the contact angle, the evaporation flux roughly varies as d -1 where d is the average distance between the evaporation front inside the porous medium and the exposed surface (bottom side in Figure 7) of the porous medium. Hence the deeper the evaporation front inside the porous medium, the lower the evaporation rate. When the invasion pattern is flat and compact (θ=120° in Figure 7), the evaporation front coincides with the flat invasion interface. For a given number of invaded pores (or equivalently a given overall liquid saturation), this corresponds to the deepest evaporation front position that is possible to obtain since there is no "holes" in the domain occupied by the remaining liquid.

Hence this type of invasion pattern leads to the largest drying time. In the other cases (transition region and IP regime), the evaporation front does not coincide anymore with the interface between the two fluids and this leads to faster drying time. Note that the equation governing the transport of water vapor in the gas phase is the Laplace equation .

0 = ∆ v P
It is well known that the transport governed by the Laplace equation is characterized by the phenomenon of diffusional screening when it takes place from an irregular interface, e.g.

[14] and references therein. Here this means that the local evaporation flux density is not uniform along the interface but strongly varies from one place to another. For example in the case of the IP patterns shown in Figure 7 (θ=80°), the evaporation flux density is quasi-null at the boundary of the fjords located deep inside the porous medium due to diffusional screening and the evaporation flux density is in fact only significant at the menisci located along the outermost boundary of the liquid / gas region shown in Figure 7.

A simple manner of illustrating the diffusional screening effect is to shown the invasion pattern together with the vapor partial pressure field v P . This is shown in Figure 7. As expected the most marked screening effect is obtained for the IP regime where the interface is the most ramified and the numerous liquid clusters, close to each other, contribute to the screening phenomenon. From the vapor isoconcentration lines shown in Figure 7, it is obvious that the "effective" evaporation front is much less deep inside the porous model for the IP regime than for the FF regime (for the same overall saturation in both systems). Hence the faster drying observed in the IP regime is a consequence of the diffusional screening effect together with the fact that bursts are the local dominant instability mechanisms in this regime (which lead to the IP highly ramified structure). When the growth of the interface becomes more coherent, the screening effect is less effective as the ramified fjords of the IP regime transform into large gulfs. This is illustrated in Figure 7 (θ = 99°). As a result, the effective evaporation front is located deeper in the porous medium on the average and this leads to greater average overall drying time compared to the IP regime. This can be expressed in terms of the average finger width w studied in [START_REF] Cieplak | [END_REF]. The finger width can be interpreted here as the mean distance between two menisci. In the IP regime, w is on the order of pore size and the screening effect is very effective. In the transition region, w (the size of gulfs) increases (up to divergence at θ c ) with θ and the screening effect becomes less and less effective.

This can be characterized quantitatively through the consideration of the equivalent flat evaporation front, which can be defined as follows. Suppose that j θc (t) is the evaporation flux density at time t for the FF regime, i.e. for θ ≥ θ c . j θc (t) can be expressed as 4), the equivalent flat front position for any flux density j θ can be defined as

j θc (t) ) P P ( D ) t ( d D RT M , v vs pm c ext v ∞ - θ -         + δ = 1 (4)
        - - = θ ∞ e v , v vs mp eq G RT M j ) P P ( D d ( 5
)
Figure 8 shows the evolution of d eq (θ) /d θc (averaged over 100 realizations) as a function of the overall mean saturation (the overall saturation is the volume fraction of the pore space occupied by the liquid). As can be seen from Figure 8, the diffusional screening is very effective for θ < θ c , which explains the faster drying as discussed above

To conclude this section, it can be observed that the influence of contact angle on average drying time can be still more marked when an intermediate low overall saturation is sought. This is shown in Figure 9 where

θ S t
is the average drying time needed to reach the mean overall saturation > < S when the contact angle is θ. As can be seen from Figure 9, the evolution of the relative difference in drying time as a function of mean overall saturation is not monotonic. The maximum difference is reached at a low saturation close to 0.5 for our computations. The decrease in the drying difference for lower saturations is associated with the drying of the last row of pores, see [START_REF] Chapuis | [END_REF] for some details on this very last period of drying.

B. Influence of contact angle on statistical fluctuations

As shown in [START_REF] Chapuis | [END_REF], the drying time is subject to significant statistical fluctuations in the IP regime whereas the drying time variability was nearly null for the FF regime, which is obvious from the invasion pattern. More precisely, it was found that the distribution of overall drying time was nearly Gaussian with σ / IP t 0 ≈ 0.1 for a 25 x 25 network, where σ is the standard deviation of overall drying time over many realizations. Hence it is interesting to study how the magnitude of statistical fluctuations varies in the transition zone. Figure 10 shows the evolutions of σ / 0 t in the transition zone. As can be seen the ratio σ / 0 t increases in the transition zone up to the sudden decrease when the invasion pattern corresponds to the FF regime, i.e. when there is no more spatial fluctuations in the growth of the interface. Hence it is remarkable to observe that the statistical fluctuations of the drying time are important up to θ c , which is consistent with the existence of random invasion patterns below θ c . The behavior of σ / 0 t near θ c is shown in Figure 11 and is consistent with a linear variation of the form

σ/ 0 t ) / . ( c θ θ - ∝ 1
. Figure 12 shows the evolution of PDF (probability density function) of overall drying time in the transition region. As mentioned earlier, the computational time of one realization is significantly greater than for the simpler drying models used in [START_REF] Cieplak | [END_REF]. As a result, the PDFs shown in Figure 12 were obtained considering 10000 realizations of a 15x15 model porous medium and not 100000 as in [START_REF] Chapuis | [END_REF] (recalling that only the IP regime was considered in [START_REF] Chapuis | [END_REF]). Despite the more limited number of realization considered here, a clear evolution of the overall drying time PDF is observed through the transition region. For θ ≥ θ c , the drying time fluctuations die out and the PDF is a dirac delta function (not shown in Figure 12). For the IP regime (θ = 80° in Figure 12), the PDF is consistent with the nearly Gaussian distribution obtained in [START_REF] Chapuis | [END_REF] for a 2D square network of interconnect channels. As noted in [START_REF] Chapuis | [END_REF], the drying time distribution is in fact slightly dissymmetric and skewed to the left for the IP regime. As can be seen from Figure 12, this dissymmetry of the PDF significantly increases with the contact angle in the transition region. This can be also shown from the values of the skewness reported in Table II. Note the tail on the left hand side for θ = 99° and the complete dissymmetry for θ = 102°, i.e. as θ c is approached from below. Although probably difficult to obtain experimentally, the dissymmetry of the drying time PDF is a clear signature of the transition zone, i.e. of the changes in the local interface growth mechanisms.

These results are for the overall drying time and it can be surmised that the effect of statistical fluctuations can be still greater for intermediate average drying time. This is illustrated in Figure 13 which shows, for several values of contact angle, the evolution of the standard deviation σ s of the overall liquid saturation over 100 realisations of a 25 x 25 disks array as a function of the average saturation S over the 100 realizations, whereas the insert in Figure 13 shows how S varies with time t . As can be seen from Figure 13, the maximum magnitude of σ s is similar for all values of contact angle below θ c and close to 0.08. However, the average saturation at which this maximum is observed increases with the contact angle.

It can be concluded that the drying time fluctuations are significant below θ c and die out at θ c . However, one should recall here that the liquid films have been ignored in the analysis. As stated in the introduction, it is reasonable to neglect the effect of films for the relatively high contact angles of interest in this study, see [START_REF] Chapuis | [END_REF] for more details. Interestingly, the results presented in [START_REF] Chapuis | [END_REF] about the influence of liquid film in the IP regime indicates that the statistical fluctuations are strongly dampen when the films can develop over a significant distance toward the surface of porous medium ahead of the liquid cluster region. Although the details remain to be elucidated, this effect is due to an enhanced screening effect of the heterogeneity of liquid phase distribution in the liquid cluster region. This suggests that the drying time statistical fluctuations and related fluctuations are significant only over a limited range of contact angles. The upper bound of this range is well defined and is given by θ c . The lower bound is less well defined since the results presented in [START_REF] Chapuis | [END_REF] (noting in passing that the "critical" angle θ c used in [START_REF] Chapuis | [END_REF] should not confused with the angle θ c considered in the present study) suggests that the damping of fluctuations due to films is not abrupt but takes place over a (albeit narrow) range of contact angle.

IV CONCLUSIONS

In this paper, we have studied the influence of contact angle on slow evaporation in 2D model porous media. For sufficiently low contact angles, the drying pattern is fractal and can be predicted by a simple model combining the invasion percolation model with the computation of the diffusive transport in the gas phase. The overall drying time is minimum in this regime and is independent of contact angle over a large range up to the beginning of a transition zone.

For the cases considered in this study, the transition zone spans over about 10° in contact Liquid films were neglected in the analysis. From previous studies, e.g. [START_REF] Chapuis | [END_REF], [9] it is expected that films do not affect significantly the invasion patterns but do affect the mean drying time and drying time fluctuations. However, this should no affect the analysis of the transition region presented in this study since film effects are expected to be non negligible only for sufficiently low contact angles, i.e. for contact angles expected to be significantly lower than θ c .

Our study was restricted to a square lattice of discs with random radii. Disks on a triangular lattice were also considered in [START_REF] Cieplak | [END_REF] and the variations of invasion pattern with θ c exhibited the same behavior as the square lattice. Thus, evaporation in triangular lattices as well as in other similar random systems is expected to exhibit the same trends as the system considered in the present study. This is illustrated in part by the experimental results shown in Figure 1 as well as those reported in [15] for a 2D network of channels of rectangular cross section.

As pointed out in [START_REF] Cieplak | [END_REF], a transition region is expected to be observed also in 3D. According to [START_REF] Cieplak | [END_REF], the critical contact angle θ c is, however, expected to be lower in 3D. The main difficulty with a 3D model is that an accurate description of the evolution of the liquid-gas interface shape becomes much more involved. Despite this difficulty, it is certainly desirable to consider the 3D case, which is obviously more representative of real porous media, since it is well known that there are significant differences in the structure of the phase distribution during drying (for contact angle below θ c ) compared to the 2D case owing to the fact that each phase can form a percolating cluster in 3D, e.g. [16] . To develop a 3D pore network model a evaporation for any value of contact angle, one option is to use the approximate parametric models of pore capillary entry pressure developed in previous works on two-phase flows , e.g. [17]. Work in this direction is in progress.

Also, only capillary forces were considered in the study and the increase in the drying time in the transition region and for θ ≥ θ c is entirely associated with the change in the drying pattern due to the change in the local instability mechanisms controlling the growth of the interface.

Under other circumstances, in particular when the drying pattern is not only controlled by the capillary forces but by the competition between capillary forces and gravity or viscous forces, the increase in drying time with the contact angle can be due to a completely different mechanism, i.e. a change in the invasion pattern due to the reduction in capillary forces compared to gravity or viscous forces (and therefore not to a change in the interface local growth mechanisms). As discussed for instance in [18], it is well know that gravity forces as well as viscous forces can contribute to stabilize the drying patterns, i.e. to the formation of a drying front. This should be kept in mind when analyzing experimental data, especially when θ c is close to 90°. For example, the change in the saturation profiles observed in the drying experiments analyzed in [19] between a random packing of 250µm beads with receding contact angles of 0° and a similar packing but with a receding contact angle of 84° (referred to as "hydrophobic" beads in [19]) could be interpreted at first glance as consistent with the transition from IP pattern to flat traveling front analyzed in the present study. However, as pointed out in [19], the influence of gravity effects is important, compared to capillary forces, with the "hydrophobic" beads. We note that the front is not sharp in this experiment with the "hydrophobic" beads and the saturation profiles are in fact quite consistent with a 3D version of the gravity stabilized front discussed in [18]. Hence, it can be surmised that θ c > 84° for the system studied in [19]. Naturally, both effects, i.e. the change in the local growth mechanisms and the change in the competition between capillary forces and gravity or viscous forces can be responsible for the increase of drying time with the contact angle in real systems.

List of Tables and Figures

.

  where d θc , D ext and D pm are the position of the flat front (= distance between the front and the porous medium surface), the effective diffusion coefficient in the external boundary layer and the porous medium effective diffusion coefficient respectively. Defining the external mass transfer resistance as G e = ext D δ , it can be seen from Eq. (4) that j θc (t) Hence from the numerical data for θ ≥ θ c , we can plot j θc -1 as a function of values of G e and D pm . Then, according to Eq.(

  angle variation. As the contact angle increases in the transition region, cooperative smoothing mechanisms of the interface become important and the width of the liquid gas interface gulfs that form during the evaporation process increases. As a result, the diffusional screening phenomenon becomes less effective and the mean overall drying time increases up to an upper bound which is reached at a critical contact angle θ c . Above θ c the drying pattern is characterized by a flat traveling front and the mean overall drying time becomes independent of the contact angle. Drying is found to be almost 65% slower for contact angles above θ c compared to the minimum drying time observed for contact angles below the contact angle marking the beginning of the transition region (~ 90° in our case). Below θ c the drying time fluctuations induced by the spatial fluctuations in the liquid phase distribution associated with each realization are significant. The drying time distribution standard deviation is of the same order of magnitude regardless of the value of contact angle in this range. However, there is a significant change in the drying time PDF. The PDF is nearly Gaussian at the beginning of the transition zone and becomes more and more skewed to the left as the contact angle increases in the transition zone. The fluctuations die out abruptly at θ c as the liquid gas interface becomes a flat front.

Figure 1

 1 Figure 1 Example of phase distribution during drying in a model porous medium made of a monolayer of 1 mm glass beads randomly distributed between two glass plates. Vapor escapes from the four sides of model porous medium: a) θ ∼ 30-40° b) θ ∼ 105°-107°. The black areas correspond to liquid saturated regions. Glass beads forming the porous medium are visible in the regions (in grey) invaded by the gas phase.

Figure 2 Figure 3

 23 Figure 2 Model porous medium formed by a regular array of disks of random diameter. Distribution of fluids at t = 0 (invading fluid in grey, displaced fluid in white). The interface is visible as a series of arcs joining the bottom row of disks. Vapor escapes from bottom edge of system

Figure 4

 4 Figure 4 Probability (in %) of interface local growth mechanisms as a function of contact angle in the displaced fluid.

Figure 5 Figure 6

 56 Figure 5 Growth of interfacial arcs during the invasion. The arrows indicate the meniscus displacement direction. Local invasion events: a) burst, b) touch, c) and d) arc coalescence (overlap)

Figure 7

 7 Figure 7 Example of drying patterns (liquid phase in dark gray, gas phase in light gray) in the transition region for one realization of a 50x50 disks array (system D). Value of contact angle

Figure 8

 8 Figure 8 Equivalent flat front position as a function of mean overall saturation for various values of contact angle.

Figure 9

 9 Figure 9 Evolution of ) ( t t ) ( t c o S c S θ -θ θ as a function of average overall saturation. ) ( t c θ 0

Figure 10

 10 Figure 10 Standard deviation σ of overall drying time as a function of contact angle in the transition region for different porosity Figure 11 Evolution of overall drying time standard deviation near θ c (~101°). The inset shows the variation of σ over the transition zone

Figure 13

 13 Figure 13 Evolution of standard deviation of overall liquid saturation σ S as a function of average saturation S . The inset shows the evolution of S as a function of time.

  Figure 9

Figure 12 Figure 13

 1213 Figure 12.

Table I

 I Parameters of the system studied (the lattice spacing a is used as reference length

	System	Range of radii	scale)	Porosity	θ
	A	0.32-0.349		0.63	98
	B	0.29-0.349		0.65	99
	C	0.27-0.349		0.67	100
	D	0.23-0.35		0.70	101
	Table II Standard deviation (σ), skewness (sk) and kurtosis (β) for the drying time
	Contact angle θ (°)	distributions shown in Figure 12. σ/ ) ( t θ 0 sk	β
	80	0.145		-0.243	2.515
	94	0.166		-0.339	2.556
	99	0.179		-0.705	2.947
	102	0.100		-2.487	10.31

ACKNOWLEDGEMENTS

Financial support from GIP ANR : (project ANR-06-BLAN-0119-01 « Intensifilm ») is gratefully acknowledged. We thank N.Sghaier for permission to present the experimental patterns shown in Figure 1.


OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible.

Study of the behavior of bi-axially loaded composite laminates

damien.busca@enit.fr, marina.fazzini@enit.fr, bernard.lorrain@enit.fr, sebastien.mistou@nimitech.fr, moussa.karama@enit.fr, marie.laetitia.pastor@iuttarbes.fr

Keywords: Composites, bi-axial, thermography, Digital Image Correlation

Composite materials are increasingly used in the industry. Combining light weight and high resistance, they represent a good alternative to metals for weight saving in structures, such as airplanes. The behavior of those materials is intensively studied by researchers all over the world. Since a structure is usually submitted to complex loading, the study of the bi-axial behavior of those materials is a key issue to the understanding of damage progression of in service composite structures. Very few experiments have been done on bi-axial loading. The particular shape of the specimen and the absence of standard specimen make the task complex and the study remains open.

A specific cruciform specimen has been designed to respond to some criteria to allow performing reliable testing, based on already existing geometries [START_REF] Escárpita | Modified cruciform specimen for bi-axial testing of fibre reinforced composites[END_REF][START_REF] Smits | Design of a cruciform specimen for biaxial testing of fiber reinforced composite laminates[END_REF]. The design has been made using a finite element model. The design criteria are:

 Stresses focused in the central area of the specimen  Fracture in the central area  Homogeneity of stresses in the central area

The specimen is then manufactured using the liquid resin infusion process, with a specific mold and a specific strategy to match with all the requirements. The material tested is a woven balanced carbon fabric with an epoxy resin. The specimen has a reduced section in the central area.

Damage characterization of materials is a complex procedure which needs to know a lot of data to be performed correctly. Digital image correlation (DIC) is an optical full field measurement method which is increasingly used [START_REF] Brault | Bi-axial fatigue analysis by stereocorrelation measurement[END_REF][START_REF] Lamkanfi | Strain distribution in cruciform specimens subjected to biaxial loading conditions -Part 1 -Two-dimentional versus 3-dimentional finite element model[END_REF]. When strain gages allow an accurate measurement of the strains on a small surface, DIC allows the measurement of the displacements all over the surface of a specimen. The accuracy of the system used is 5 μm. It is also possible to calculate the strains all over the surface.

Damage is also characterized by a local increase of the temperature of the material, due to energy release caused by fracture of the matrix or the fibers. These increases are detectable by infrared thermography, which is also a full field optical measuring method.

Combining those two monitoring systems, a specific monitoring setup has been designed, to get the most information during a bi-axial static tensile test. Two cameras are set on the opposite sides of the specimen to perform DIC and have access to the displacement fields on both sides of the specimen. A thermic camera is also used to get the temperature evolution on the surface of the specimen until complete fracture. Finally, a high-speed camera is used to get successive pictures during the fracture of the specimen at a rate of 3000 images per second to study accurately the fracture kinematic. The specimen is tested on a bi-axial tensile machine. The test consists in a bi-axial tension with the same load on both arms. The load is applied continuously until complete fracture of the specimen.

Finally, X-ray tomography is performed on the broken specimen to determine the fracture profile. X-ray tomography is an efficient way to perform non-destructive control by obtaining a 3 dimensional picture of the specimen, including the inside. It is then pssible to detect delamination, as well as fibers and matrix cracks.

All the data obtained are correlated to attempt to identify the damage mechanisms of those materials under bi-axial loading. A finite element approach is then attempted with a progressive damage model to match with the experience.

  D. Busca 1 , M. Fazzini 1 , B. Lorrain 1 , S. Mistou 1,2 , M. Karama 1

	1 : UNIVERSITE DE TOULOUSE, INP,
	Ecole Nationale d'Ingénieurs de Tarbes, Laboratoire Génie de Production, 47 avenue
	d'Azereix, BP 1629 65016 TARBES CEDEX, FRANCE
	2 : Nimitech Etudes, Parc Industriel de la Haute Bigorre, 1 Avenue des victimes du 11
	juin 1944, 65200 BAGNERES DE BIGORRE, France


T Fourcade

O Dalverny

J Alexis

C Seguineau

J M Desmarres

Parametric identification of elastic-plastic constitutive laws using spherical indentation

Keywords: Parametric identification, Finite Element Modelling, Spherical tip, instrumented indentation, pile-up

In this paper we present a methodology to characterize elastic-plastic constitutive law of metallic materials using spherical indentation tests and parametric identification coupled with finite elements model (FEM) simulation. This procedure was applied to identify mechanical properties of aluminium, copper and titanium and the identified models show differences with experimental uniaxial tensile tests results that do not exceed 10%. A sensitivity study according to a 2 Design of Experience (DoE) was achieved to determine which data that can be extracted from pile-up is the most relevant to use in order to enhance the identification procedure. It appears that the maximum pile-up height seems to be the best suited for this purpose.

Introduction

During the lasts decades, extensive research has been devoted to the characterization of elasticplastic properties for a wide range of materials using indentation techniques. Two main approaches can be noticed in the literature. The first one consists in the resolution of an analytical direct problem [1,2,3] and the second one consists in the resolution of an inverse problem [4,5,[START_REF] Bolzon | NanoScience and Technology[END_REF][START_REF] Dowhan | [END_REF]8,9] using or not Finite Elements Model (FEM). The aim of the study is to develop a methodology to characterize elastic-plastic constitutive law of metallic materials using spherical indentation tests and parametric identification coupled with FEM simulation.

In the first part of the paper, the identification procedure is described and the results of the identification for three materials are compared to those obtained through uniaxial tensile test in order to validate its accuracy and robustness.

The second part is dedicated to discuss the possibility of adding experimental information from the pile-up after unloading as a further information to use into the identification procedure.

2.

Identification of constitutive laws using spherical indentation

Tested materials and experimental procedure

The three tested materials are aluminium 1050-O, temper hard copper and mill annealed titanium T40. The titanium has a hexagonal close-packed (hcp) lattice whereas both the other materials have a face centred cubic (fcc) structure.

• Standard characterization

The identification results have been compared to those of a well-known normalized test. For each materials, at least 7 tensile tests have been performed according to the ISO NF EN 6892-1 and the results are summarized in Table 1. All the uncertainties on mean values given in this paper are calculated with a confidence interval of 95% using the Student law. The elastic-plastic behaviour of the three material cannot be represented by a unique constitutive law. Both the aluminium and the copper can be described using Ramberg-Osgood expression [10] which allows a soft elastic-plastic transition:

= + 0,002 , ( 1 )
This model involves three independent parameters: the modulus of elasticity , the conventional yield stress . and the hardening parameter !.

Because of its linear hardening, the titanium is better described by a modified version of the Rasmussen model [11,[START_REF] Fourcade | Etude des propriétés mécaniques de matériaux métalliques en couches minces[END_REF] which is based on the model of Ramberg-Osgood:

= " + 0,002 , si ≤ , & , + 0.25 ' ( ) * + . + 0.002 si > , ( 2 )
Until reaches . , the Rasmussen model is the same as the Ramberg-Osgood expression. When stress has reach yield stress, the transformed stress is defined as:

= -. ( 3 )
In order to insure a continuity between both the parts of the curve, the . modulus is calculated as the tangent modulus of the stress-strain curve when = . . It is then obtained from equation ( 1 ) as -/-| 0 . :

. = 12 . ' 3 . 4 ) ( 4 )
The transformed ultimate stress is calculated according to equation ( 3 ) as :

= -. ( 5
)
The Rasmussen model is then constituted of 5 independent parameters.

• Indentation tests: The indentation tests have been performed using a 1mm diameter spherical tip. The maximum load is 150 N for the aluminium and 300N for both the copper and the titanium and ten imprints were performed for each material. During the test, the load 5 has a function of the penetration ℎ of the tip has been continuously monitored. Figure 1 shows the mean indentation curve of the indentation tests performed on aluminium. The error bar represent the uncertainties on the mean value with a confidence interval of 95% using Student Law.

Identification algorithm and numerical model

The identification of the parameters values of the behaviour law was achieved using the results from the indentation tests and corresponds to the resolution of an inverse identification problem. A usual methods of resolving such a problem is summarized in Figure 2. The considered constitutive laws are the Ramberg-Osgood one for the aluminium and the copper and the Rasmussen one for the titanium. The modulus of elasticity is usually calculated with the earlier stages of the unloading slope [START_REF] Oliver | [END_REF] and we consider it as known during the identification. The values used in the identification for the three materials are the ones calculated from the tensile tests shown in Table 1. Thus, the number of parameters to evaluate with the identification for the Ramberg-Osgood law and the Rasmussen law are respectively 2 and 4.

The inputs of the problem are the experimental response from the indentation tests, the . 5 GKR is the mean experimental indentation curve, 5 S is the simulated one, 5 TUV S and 5 *WK GKR are respectively the maximum of the experimental and simulated indentation curve. F G and F * are weights and have respectively a value of 1 and 2 in order to maximize the weight of the value at maximum depth, which is the one with the smallest uncertainty. Specimen roughness, presence of surface hardening and tip defect may have a significant influence on the earlier stages of the loading curve. In order to minimize such effects, the 10 first percent of the loading curve are not taken into account in the calculation of the difference vector. The cost function is then defined as the Euclidean norm of the difference vector:

X><? = Y=><?Y ( 7 )
The identification is led by minimising the cost function using the Levenberg-Marquardt algorithm from the scipy.optimize python library. As a descent algorithm, the Levenberg-Marquardt algorithm is sensitive to local minima and thus the choice of the initial parameters value is critical for the algorithm not to converge towards the global minimum. The initial parameters are then chosen close to the global minimum using a Monte-Carlo stochastic random method. The number of draws is function of the number of parameters Z to evaluate as 2 C with a minimum of 16 draws. The validity domain of the parameters are given in Table 2. The parameters corresponding to the minimum value of the cost function are then used as initial parameters for the Levenberg-Marquardt algorithm. The numerical model of the indentation test has been achieved using Abaqus. It consists in an axisymmetric model [4,14,15,[START_REF] Dowhan | [END_REF]16] with a refined "ALE adaptative meshing" of the area under the tip in order to minimize distortion of the elements. The elements size is linked to the maximum depth of indentation. The elastic-plastic behaviour of the material is a multi-linear model with an isotropic hardening. The specimen is made of CAX3 elements under and near the tip and the furthest are infinites CINAX4 elements in order to take into account the far displacement field even in a "local" problem such as indentation tests. The tip is made of CAX4R and CAX3 elements and its mechanical behaviour is linear with a modulus of elasticity and Poisson coefficient are respectively equal to 240GPa and 0.33. A "Surface-to-surface" method is used to model the contact between the tip and the surface and the contact is considered frictionless [4,14]. The shape of the simulated indentation curve are similar to the experimental ones and no numerical oscillation can be noticed on the curve.

Identification results and discussion

• Results : The random Monte-Carlo procedure allows to get a set of parameter close enough to the global minimum. Table 3 shows that the Levenberg-Marquardt procedure allowed to minimize the cost function for the three materials. It can be noticed that the final values of the cost function are similar for the three materials. The final simulated indentation curve is very close from the experimental one as it can be seen in Figure 3-a. The parameters of the constitutive law for the three materials are also given in this table. • Discussion: The mean relative differences between identified constitutive law and the experimental curve for the aluminium, the copper and the titanium are respectively equal to 4.6%, 6.8% and 3%. These results show that the identified laws are in good agreement with experimental tensile test results. Nevertheless, the experimental investigations of the imprint's shape shows a slight difference between experimental and simulated pile-up. The maximum height differs from about 25%. In the next section we use pile-up information to enhance the identification procedure.

Pile-up influence on the identification

For this study, we focused on the aluminium because Ramberg-Osgood law is simpler than the Rasmussen one and unlike the copper it shows great elongation at break and notable hardening. Some authors already use data extracted from pile-up after the unloading as an additional information in order to approximate more accurately the elastic-plastic behaviour like topography of the surface around the imprint [17,[START_REF] Bolzon | NanoScience and Technology[END_REF] or the maximum pile-up height as additional data [18,16].

In order to measure various quantity from the pile-up, the topography of the imprint have been measured after unloading using an optical interferometer. Optical technology measurement allows to get a three dimensional information around the imprint, but because of the great slope into the imprint, the information cannot be measured here as shown in Figure 4. Moreover, the determination of the initial plan is also difficult due to the roughness of the surface. The simpler quantity to measure which is representative of the pile-up is its the maximum height ℎ *WK . The pile-up volume seems to be also representative of the pile-up, but the missing part in the middle of the imprint and the determination of the initial plan makes difficult a precise measurement of the pile-up volume. In order to overcome this difficulty we defined 4 different strategies for measuring the volume of pile-up which are summarized in Figure 5. For the totaltotal volume \ ]] , all the volume is measured and the non-measured part into the imprint is neglected. The total-partial volume \ ]H is measured between the maximum height and 80% of the maximum height. The beginning of the measurement of the partial-total volume \ H] is the abscissa of the maximum height and the partial-partial volume \ HH is calculated between the abscissa of the maximum height and 80% of the maximum height.

Sensitivity study

In order to determine which of the pile-up parameters it is efficient to add to solve the identification problem, we performed a numerical sensitivity study according to a 2 design of experience (DoE). The factors of the DoE are the maximum indentation depth, the yield stress

. and the hardening parameter !. Their values are respectively _0.1, 0.2`mm _30, 100`MPa and _3,20`.

The responses are ℎ *WK , \ ]] , \ ]H , \ H] and \ HH , Figure 6 shows the result of the DoE. It can be noticed that the maximum height ℎ *WK is quite sensitive to the material parameters whereas the volumes are very sensitive to the indentation depth but less to the yield stress. Indentation test indeed introduces plastic strain which occurs with a constant volume. Thus, the volume of pile-up is mainly driven by the indentation depth and is linked to the volume of the tip penetrating into the material. The use of pile-up volume as an extra information to resolve our inverse problem is then not relevant. The better way seems to use the maximum height of pile-up.

Identification with maximum pile-up height

An identification has been achieved using the maximum pile-up height. The experimental measure of ℎ *WK has been used and the difference vector has been defined as: =><? = @A 1 ><?, . . . , A C ><?, … , A ><?, D 1 ><?, … , D C ><?, … , D ><?, a 1 ><?, … , a C ><?, … , a ><?E ( 8 ) where . The values of F G , F * and F c are respectively 1, 1 and 2. The cost function X(<) is also defined as the Euclidean norm of =(<).

• Results and discussion: The final value of the cost function is 0.39 which is higher than the one without taking into account ℎ *WK into the cost function. The identified . and ! parameters of the Ramberg-Osgood are respectively 42 and 4.3 and the mean relative difference between the experimental tensile curve and the identified Ramberg-Osgood law is 6.3%. This difference is twice the difference of the identification without taking the pile-up maximum height into account.

The difference between the experimental and simulated ℎ *WK is 7.82 × 10 Lf which is really low. On the other hand, the difference between the experimental and simulated curve is much higher than the one without taking pile-up into account and the shape of the indentation curve are also more different. 4.

Conclusion

In this paper, a methodology to identify elastic-plastic constitutive law of metallic materials using inverse method with spherical indentation and FEM is proposed. Three materials have been tested and the results seems to be in good agreement with the reference tensile tests. Furthermore, topography measurement of the surrounding of the imprint were performed and we investigated which information about pile-up the most relevant information is the most relevant to use in the identification procedure in order to get more accurate results. It appears that the parameter which is the most sensitive to the material maximum height of pile-up whereas the information about pile-up volume are mainly linked to the indentation depth. The maximum height of pile-up has been introduced into the identification procedure and used to identify the mechanical properties of aluminium. However the results of the identification shows that the identified behaviour is further from the tensile tests when taking into account the pile-up maximum height. The comprehension of this phenomenon is still under investigation.

Figure 1 :

 1 Figure 1: Mean load-displacement curve of indentation test on aluminium specimen.

Figure 2 :

 2 Figure 2: Block diagram of the identification of constitutive law from indentation tests.

  the tested material associated with a set of initial parameter and the numerical model of the test. For a set of parameters <, a difference vector =(<) is then calculated as: =><? = @A 1 ><?, … , A C ><?, … , A ><?, D 1 ><?, … , D C ><?, … , D ><?

Figure 3 -Figure 3 :

 33 Figure 3: Experimental and simulated indentation curves for aluminium (a), experiment and identified tensile curves for aluminium (b), copper (c) and titanium (d)

Figure 4 :

 4 Figure 4: Profile of an imprint measured with an optical interferometer on aluminium.

Figure 5 :

 5 Figure 5: The different strategies for measuring the volume of the pile-up.

Figure 6 :

 6 Figure 6: Results of the 2 DoE.

Table 1 :

 1 Tensile tests results.

		Elastic Modulus (GPa)	Yield stress . (MPa)	Ultimate strength (MPa)	Ultimate strain (MPa)
	Aluminium Copper Titanium	69 ± 4 106 ± 9 107 ± 2	59 ± 4 262 ± 39 222 ± 3	120 ± 9 279 ± 39 413 ± 50	0.26 ± 0.03 0.12 ± 0.01 0.42 ± 15

Table 2 :

 2 validity domain for Ramberg-Osgood and Rasmussen models during identification process.

	Ramberg-Osgood	. (MPa) Min Max 10 350	Min 1.1	!	Max 60	Min	[	Max	Min	(MPa) Max
	Rasmussen	10	350	1.1		60	1.1		3	20	850

Table 3 :

 3 Results of the inverse problem resolution for the tree materials.

		63	51	0.10	53	5.7
	Copper	0.65	51	0.07	211	9.8
	Titanium	0.79	29	0.16	173	9.6	370	2.57