Ali Daouadji

Pierre-Yves Hicher

Afif Rahma

An elastoplastic model for granular

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The mechanical behaviour of granular materials is dependent on the properties of the grains of which they are constituted. Hence, it is natural to search for correlations between the parameters of the constitutive model of a supposedly continuous medium and the properties of the grains (form, size, granulometry. . . ). In order for these correlations to have meaning it would be necessary to introduce a system of understanding from the discontinuous to the continuous medium. This methodological work has notably been developed by Biarez and Favre (1977) and [START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF]. It has been applied for determining the parameters of an elastoplastic model. Correlations between the parameters of the grain, on one hand, and the parameters of the elastoplastic model, on the other hand, have been proposed and validated in the case of sands which possess diverse grain size distributions [START_REF] Hicher | Micro-macro correlations for granular media. Application to the modelling of sands[END_REF].

The purpose of this article is to re-analyse, in the spirit of this approach, the influence of the granulometric characteristics on the behaviour of the grain assembly, particularly in the case where the above would evolve in the course of mechanical loading under the effect of grain ruptures. The aim is to determine the mechanical characteristics of the medium considered as continuous which would become integrated within the constitutive model by using complementary knowledge about the properties of the discontinuous medium.

We shall therefore analyse the evolution of behaviour under the effect of granulometric change. We shall then present a method that takes into account progressive grain ruptures through a systematic actualisation of the parameters of the constitutive model, by making them dependent on the evolution of the physical properties of the discontinuous medium.

Methodology for classifying parameters: from a discontinuous to a continuous medium

If we consider a grain assembly in which only the grains are treated as continuous media, the equilibrium of the assembly submitted to exterior forces could be calculated by the traditional method in continuum mechanics [START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF]. The solution to the problem provides mechanical properties of the granular material considered as a continuous medium. By use of this method, we can deduce that the mechanical properties of the granular medium will depend on the following two groups:

• the mechanical properties of the grains themselves,

• boundary conditions of a solid formed by grains which depend on the geometry of grains and on the geometry of their arrangement. The geometry of grains may be defined by the parameters of size, form, and grain size distribution. The geometry of the arrangement consists of a parameter which expresses compacity and parameters which describe its anisotropy. In all cases where grain rupture may be considered as negligible, the irreversible deformations are due only to a change in the arrangement of particles and not to a change of the particles themselves. The first group of parameters will therefore be constant, only the parameters of the arrangement will evolve.

Under these conditions, we can distinguish two classes of parameters: Class 1: constant parameters independent of the arrangement also called nature parameters.

(a) parameters of the constitutive law of the grains:

• of the grains themselves, e.g., the linear elasticity E g ,

• of the contact between the grains: interparticulate friction angle ϕ g , bonds between the grains C g ;

(b) parameters of the geometry of the grains:

• size, grain size distribution. . . ,

• shape: angularity. . . ,

• surface state.

Class II: variable parameters describing the arrangement geometry • the isotropic part: void ratio (e), compacity I D ,

• anisotropy: its description requires in general a tensorial representation.

This analysis of the discontinuous / continuous medium passage allows us to establish a framework of expression for the correlations between parameters. It will hence be possible to obtain correlations either within the same class, or among classes by using the following scheme:

General

Constitutive law Boundary Solution + + ⇒ equations of grains conditions balance (geometry of (parameters of equations . . . arrangement) continuous medium)

(1) and in the absence of grain rupture:

Nature (I) + Compacity (II) ⇒ Rheology of equivalent continuous material. (2)
The above relations allow us to build a framework for the correlations among the parameters representative of the grains and the parameters of a constitutive law of the continuous medium [START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF] and [START_REF] Hicher | Micro-macro correlations for granular media. Application to the modelling of sands[END_REF].

Influence of grain ruptures on the behaviour of granular materials

Grain breakage is an important subject in Soil Mechanics. [START_REF] Fukumoto | Particle breakage characteristics of granular soils[END_REF] qualified it as a subject of major importance as the soil undergoes significant changes due to this phenomenon. There are many practical problems related to grain breakage such as the stability of earth dams, the bearing capacity of piles. . . . We may cite several investigations on this subject: [START_REF] Roberts | The compressibility of sand[END_REF][START_REF] Marsal | Large scale testing of rockfill materials[END_REF][START_REF] Lee | Compressibility and crushing of granular soils in anisotropic triaxial compression[END_REF][START_REF] Lee | Drained strength characteristics of sands[END_REF][START_REF] Vesic | Behavior of granular materials under high stresses[END_REF]. . . and more recently [START_REF] Colliat-Dangus | Comportement des matériaux granulaires sous fortes contraintes. Influence de la nature minéralogique du matériau étudié[END_REF][START_REF] Fukumoto | A grading equation for decomposite granite soil[END_REF][START_REF] Fukumoto | Particle breakage characteristics of granular soils[END_REF][START_REF] Hagerty | One-dimensional high-pressure compression of granular media[END_REF][START_REF] Coop | The Behaviour of Granular Soils at Elevated Stresses[END_REF][START_REF] Kim | Etude expérimentale du comportement mécanique des matériaux granulaires sous fortes contraintes[END_REF]Yamamuro andLade, 1996a, 1996b;[START_REF] Lade | Significance of particle crushing in granular materials[END_REF][START_REF] Dowell | The fractal crushing of granular materials[END_REF][START_REF] Leung | The role of particle breakage in pile creep in sand[END_REF][START_REF] Dowell | On the micromechanics of crushable aggregates[END_REF].

In the preceding chapter we emphasised the hypothesis of invariance of the parameters which characterise the geometry of the grains. When significant grain ruptures occur during mechanical loading, this hypothesis is no longer valid. Grain rupture may be classified according to three modes as showing in figure 1 [START_REF] Guyon | Du sac de billes au tas de sable[END_REF]:

• fracture: a grain breaks into smaller grains of similar sizes,

• attrition: a grain breaks into one grain of a slightly smaller size and several much smaller ones,

• abrasion: the result is that the granulometry remains almost constant but with a production of fine particles. Hence, the rupture of sharp angles may be classed in the mode of rupture by attrition, whereas the rupture of micro-asperities during the sliding of the grains may be classed in the mode of abrasion rupture. This may occur in the absence of fracture or attrition notably, in the case of cyclic tests of small stress intensity.

In this study, we have concentrated on the two first modes of rupture. They both lead to a significant evolution of the grain size distribution curve. This phenomenon brings about an evolution of the mechanical properties of the granular material, which we may summarise as follows: Grain ruptures on an isotropic or one-dimensional compression path produce an increase of compressibility, which is represented by a slope increase in the e-log p plane (e is the void ratio and p the mean effective stress). In the case of identical grains such as an assembly of glass balls of the same diameter, the compressibility increase is very abrupt since the ruptures occur simultaneously. The phenomenon is more progressive in sand on account of the grain's size and shape differences, which produces a less homogeneous division between the contact forces. In figure 2 we can see the superposition (in the e-log p plane) of several one-dimensional compression tests on a sand for stresses up to 100 MPa. We obtain first for very loose samples a compression index (slope Cc in the e-log p plane) in the range of 0.15 to 0.20. The value of this slope increases for values of the mean effective stress p > 1 MPa as a result of grain crushing. We obtain then a slope of apparent consolidation whose value can double or triple in relation to standard values. This slope increase depends clearly on the amplitude of grain rupture [START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF].

On a deviatoric path, for example a compression triaxial test, we observe first of all a significantly bigger contractancy. This may be explained by the displacement of the perfect plasticity relation in the e-log p plane. The explanation is schematically given in figure 3 and figure 4 from [START_REF] Biarez | Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires[END_REF]. We can see in figure 3 that the values of the minimal (e min ) and maximal (e max ) void ratio decrease when the grain size distribution (C u ) increases (until C u ≈ 10). For a given shape, we obtain correlations between the position of the critical state line and nature parameters (e min , e max , grain angularity, C u = d 60 /d 10 , where d n represents the diameter corresponding to n% of material being smaller by weight) are drawn ( figure 4). One can see that an increase of C u due to grain ruptures will produce a slip of the critical state line towards smaller void ratios. In order to reach the critical state during a triaxial test, a bigger volume change will be needed in case of significant grain breakage. This phenomenon will also result in a change of the stress-strain relationship and in particular in an increase of the axial strain corresponding to the maximum strength. Therefore a significant change takes place in the global behaviour of granular media due to grain breakage, which needs to be taken into account in a constitutive model.

This may be verified by comparing the results obtained from a sand under small stresses (negligible ruptures of grains) and under elevated stresses. We observe a very strong increase in the volume change. In the adimensional plane q/Mp-ε 1 ( figure 5), where q is the second invariant of the stress tensor and M a constant of the material (value of q/p at large deformations), the test curves under elevated stresses are shifted towards the right in relation to the network which corresponds to the behaviour of a loose contractant material without grain breakage (characteristics of materials tested by Le [START_REF] Long | Contribution à l'étude des propriétés mécaniques des sols sous fortes pressions[END_REF], [START_REF] Colliat-Dangus | Comportement des matériaux granulaires sous fortes contraintes. Influence de la nature minéralogique du matériau étudié[END_REF], [START_REF] Ziani | Contribution à l'étude du compactage des sols: cas particulier du comportement des sables très peu denses[END_REF] and [START_REF] Kim | Etude expérimentale du comportement mécanique des matériaux granulaires sous fortes contraintes[END_REF] are given in table I ).

Parameters influencing grain rupture

We shall briefly analyse and classify the parameters which are likely to influence grain rupture. Two groups of parameters may be extracted:

• nature parameters, that is to say parameters related to characterising the discontinuous medium. These are parameters of class I, as defined above.

• mechanical parameters, that is to say parameters characterising the mechanical state of a medium supposed to be continuous: stress and strain tensors.

Influence of nature parameters

Several studies [START_REF] Vesic | Behavior of granular materials under high stresses[END_REF][START_REF] Colliat-Dangus | Comportement des matériaux granulaires sous fortes contraintes. Influence de la nature minéralogique du matériau étudié[END_REF][START_REF] Bard | Comportement des matériaux granulaires secs et à liant hydrocarboné[END_REF][START_REF] Kim | Etude expérimentale du comportement mécanique des matériaux granulaires sous fortes contraintes[END_REF]. . . have clearly shown the influence of nature parameters of the grains: the mechanical resistance of the grains, the form, size, and grain size distribution.

The mechanical resistance of grains may be linked, on one hand, to its mineralogical nature and on the other hand to its state of alteration. It therefore depends on the history of how the natural deposits were formed, the history of its fabrication (crushing. . . ) for the artificial mixtures. The influence of the shape may be linked to the preceding point in that the statistical distribution of faults (grain joints, fissures. . . ) depends on the shape of the particles, the smallest being more resistant than the large ones (scale effect).

Grain breakage augments with their angularity. This may be attributed to a greater fragility of the contact points of small curvature radius.

Another important factor is the grain size distribution and different studies appear to agree that a badlygraded grain size distribution enhances grain rupture. This effect may be explained by a better distribution of intergranular forces in a well-graded medium. A recent experimental study realised by [START_REF] Hicher | Experimental evidence and modelling of grain breakage influence on mechanical behaviour of granular media[END_REF] provides proof of this influence with the help of two mixtures made from crushed granite. Tested in the same mechanical conditions and for identical nature parameters except for the value of d 60 /d 10 , the well graded mixture G1 (d 60 /d 10 = 10) presented very slight evolution in its grain size distribution in contrast to the badly graded mixture G2 (d 60 /d 10 = 2). Some examples of evolution on a one-dimensional compression path are presented in figure 6. It is interesting to note that it is the totality of the grain size distribution which has evolved: all the grain sizes are concerned by the rupture phenomenon and not only the largest among them.

Influence of stress-strain history

It is evident that the intensity of applied stresses is a preponderant factor on the development of grain rupture. The level of stress from which they appear in a significant manner depends on the preceding parameters. However, the stress path also plays an important role. Experimental studies on this subject allow us to conclude that the stress deviator plays a major role and that the triaxial tests, for example, generate more ruptures than the isotropic or one-dimensional compression tests of the same average stress (see figure 7) [START_REF] Hicher | Experimental evidence and modelling of grain breakage influence on mechanical behaviour of granular media[END_REF].

This may very well be explained by the relative displacements of particles which encourage grain rupture. Figure 7 shows this phenomenon in the case of crushed granite as seen above. We observe that for a triaxial test of a lateral stress constant and equal to 15 MPa, there is an evolution of the grain size distribution between the axial strains of 15% and 30%, while the axial stress has attained a practically constant value. These results permit us to illustrate the effect of mechanical parameters associated with the continuous medium: the level and path of stress, the amplitude of strains during mechanical loading.

A brief review of existing relations

Most of the investigations cited in section 3 are experimental ones. They permit us to understand the influence of grain crushing on the behaviour of granular materials as well as the parameters influencing this breakage. Some authors have proposed relations that we can schematically divide into four categories:

(i) factors which allow the degree of grain breakage to be quantified: [START_REF] Marsal | Large scale testing of rockfill materials[END_REF], [START_REF] Lee | Compressibility and crushing of granular soils in anisotropic triaxial compression[END_REF], [START_REF] Hardin | Crushing of soil particles[END_REF] and [START_REF] Lade | Significance of particle crushing in granular materials[END_REF]; (ii) relations to determine the evolution of the grain size distribution: [START_REF] Fukumoto | A grading equation for decomposite granite soil[END_REF][START_REF] Fukumoto | Particle breakage characteristics of granular soils[END_REF], [START_REF] Dowell | The fractal crushing of granular materials[END_REF] and [START_REF] Dowell | On the micromechanics of crushable aggregates[END_REF]; (iii) relations between the maximal friction angle and the friction angle at critical state or the interparticulate friction angle: [START_REF] Vesic | Behavior of granular materials under high stresses[END_REF], [START_REF] Koerner | Effect of particle characteristics on soil strength[END_REF], [START_REF] Barton | Shear strength of rockfill[END_REF], [START_REF] Bolton | The strength and dilatancy of sands[END_REF], [START_REF] Naylor | Prediction of construction of Beliche Dam[END_REF], [START_REF] Maksimovic | A family of nonlinear failures envelopes for non-cemented soils and rock discontinuities[END_REF]; (iv) models to provide the evolution of the stress-strain relationship during loading: [START_REF] Loret | Modelization of sand behaviour over a wide stress range[END_REF], [START_REF] Simonini | Analysis of behavior of sand surrounding piles tip[END_REF], [START_REF] Dowell | The fractal crushing of granular materials[END_REF] and [START_REF] Dowell | On the micromechanics of crushable aggregates[END_REF].

Among the different grain breakage factors proposed (i) Hardin's factor (1985) is the most widely used. More recently, [START_REF] Lade | Significance of particle crushing in granular materials[END_REF] have proposed a new interesting factor based on the effective grain size. Unfortunately, these factors are computed after the tests and do not give information on the grain size distribution during loading nor on the stress-strain relationship.

It might be interesting to know the grain size distribution after grain crushing occurs (ii). For example, in order to evaluate, albeit with a gross approximation, the permeability k of the filter of earth dams, we can use the well-known Hazen's relation:

k = 100(d 10 ) 2
where k is given in cm/s and d 10 in cm.

(3) [START_REF] Fukumoto | Particle breakage characteristics of granular soils[END_REF] proposed a relation allowing us to plot the grain size distribution after testing, regardless of the type of test performed (one-dimensional or triaxial compression tests, gyratory test or compression test) or of the stress level or the number of blows applied. He assumed that the percentage in weight of a given size of grain can be expressed by a geometric progression. However, it is necessary to obtain the grain size distribution after testing in order to determine the parameters of the relation. [START_REF] Dowell | The fractal crushing of granular materials[END_REF] and [START_REF] Dowell | On the micromechanics of crushable aggregates[END_REF] have also proposed a relation between the percentage in weight passing through a given sieve and the grain size. For that, they assumed that the grain size distribution has a fractal character. However, the resulting curves are only qualitative. [START_REF] Vesic | Behavior of granular materials under high stresses[END_REF] have plotted the secant friction angle as a function of the mean stress for several sands (iii). This angle corresponds to the maximal circle in the Mohr plane and is given by (see next section for notations):

φ s = arcsin σ 1 -σ 3 σ 1 + σ 3 . ( 4
)
The value of this secant angle decreases with the increase of the mean stress and tends to a minimal value. [START_REF] Vesic | Behavior of granular materials under high stresses[END_REF] have proposed that the measured shear strength can be expressed as the strength due to sliding friction + dilatancy effects + crushing and rearranging effects. [START_REF] Koerner | Effect of particle characteristics on soil strength[END_REF] has proposed the same expression in terms of friction angles. [START_REF] Bolton | The strength and dilatancy of sands[END_REF] stated that the friction angle can be decomposed in two terms: the friction angle at the critical state and the dilatancy angle. Moreover, he proposed, for a given type of test, a relation between the maximum angle of friction and the maximum dilatancy rate. These relations are widely used. The relation proposed by [START_REF] Naylor | Prediction of construction of Beliche Dam[END_REF] connects the angle of friction to the confining pressure. [START_REF] Maksimovic | A family of nonlinear failures envelopes for non-cemented soils and rock discontinuities[END_REF] also related the friction angle to the state of stress.

The last category (iv) concerns the models that can provide the evolution of the stress-strain relationship. [START_REF] Simonini | Analysis of behavior of sand surrounding piles tip[END_REF] proposed a model composed of two yield surfaces. A first nonlinear Mohr-Coulomb criterion is used when non-significant grain crushing occurs and a second criterion is proposed by Zienkiewicz and Humpheson (1977) (cited by the author) when it is not. In both cases, the flow rules are nonassociated. Simonini also used Bolton's relation (1986) and computes Hardin's factor to quantify the degree of grain crushing. It is obvious that, in the transition zone, the model cannot estimate accurately the state of stresses and strains. The results obtained with the plot of the constant Hardin's factor are interesting even though no results in terms of stress or strain are given. In his elastoplastic model, [START_REF] Loret | Modelization of sand behaviour over a wide stress range[END_REF] has modified the well-known relation proposed by [START_REF] Roscoe | Yield of clays in states wetter than in critical[END_REF] and [START_REF] Schofield | Critical State Soil Mechanics[END_REF]:

q dε p d + p dε p v = Mp dε p d (5)
widely used for cohesionless materials or clays. The right-hand side of relation ( 5) was identified as the internal friction dissipation. Loret added a term to the right-hand side that corresponds to the breakdown work. Thus, M does not remain constant but is a function of the mean stress p. [START_REF] Dowell | The fractal crushing of granular materials[END_REF] and [START_REF] Dowell | On the micromechanics of crushable aggregates[END_REF] proposed a relation between the void ratio and the logarithm of the stress which can be used for one-dimensional compression tests. They also added, to the right hand side of the relation ( 5), a term that took into account the increase of the surface due to grain breakage and the surface energy.

As we can see, very little has been proposed concerning the modelling of grain breakage influence on the stress-strain relationship of granular materials. We present in the following an elastoplastic model developed for the modelling of soil behaviour. Inside the framework of this model, we have introduced a special treatment concerning grain breakage influence in order to enlarge the capabilities of the model.

An elastoplastic constitutive model for granular media

The constitutive model we used is an elastoplastic model which combines four mechanisms: three deviatoric mechanisms of plane strain and an isotropic mechanism. We shall present it here succinctly in the case of monotonic loading. For more information, our readers may refer to the following studies: Aubry et al. (1982) and [START_REF] Hujeux | Une loi de comportement pour le chargement cyclique des sols[END_REF].

Notations

Assuming that σ , ε and ε are respectively the tensors of effective stress, strain and strain rate, we can write the following:

p = 1 3 Tr(σ ), ε v = Tr(ε), s = σ -pI, e = ε - ε v 3 I, ( 6a
) q = 3 2 s ij s ij 1/2 , εd = 2 3 ėij ėij 1/2 .
In the k-plane, we can define effective stress, strain and strain rate tensors respectively σ k , ε k and εk . If ii represents the component "ii" of the tensor (no summation on i) then:

p k = 1 2 Tr(σ k ) = 1 2 (σ ii + σ jj ), (ε v ) k = Tr(ε k ) = ε ii + ε jj , q k = 1 4 (σ ii -σ jj ) 2 + σ 2 ij 1/2 , γk = (ε ii -εjj ) 2 + 4ε 2 ij 1/2 , ( 6b
)
s k = s k 1 s k 2 = σ ii -σ jj 2 σ ij ,
p a is the atmospheric pressure, I is the unit tensor.

As in the theory of elastoplasticity, we assume that the strain rate tensor ( ε) can be decomposed into an elastic strain rate tensor ( εe ) and an irrecoverable strain rate tensor ( εp ): ε = εe + εp .

(7)

Elastic part

The elastic part of the model is isotropic, non-linear and depends on the mean stress:

εe v = 1 K(p) ṗ, (8)
εe d = 1 3G(p) q, (9)
K(p) = K a p a p p a n , G(p) = G a p a p p a n , ( 10
)
where K and G are respectively the bulk and shear modules and n is a constant (0 n 1).

Plastic part

We assume that the plastic part of the deformation is governed by three deviatoric plane strain mechanisms in three orthogonal planes and a purely volumetric mechanism.

Deviatoric mechanisms

The equations of each mechanism k (k = 1, 2, 3) are identical. We assume that the mechanism k, in the ij -plane, produce plane strain and strain rates.

The yield surface, for each k-plane, is given by:

f k (p, p k , s k , p c , r k ) = sk -r k , ( 11
)
where the components of sk are:

sk1 = s k 1 F k (p, p k , p c ) , sk2 = s k 2 F k (p, p k , p c ) (12)
and F k (p, p k , p c ), a factor of normalisation, is given by:

F k (p, p k , p c ) = sin φp k 1 -b Log p p c , ( 13
)
r k is an internal variable allowing a progressive mobilisation of the mechanism k:

r k = r el k + γ p k 1 + γ p k . ( 14
)
Its evolution law is given by:

ṙk = λ p k (1 -r k ) 2 a(r k ) . ( 15
)
p c is the critical pressure and is related to the plastic volumetric strain ε p v by the following equation:

p c = p co exp βε p v ( 16
)
where 1/β is the slope of the critical state line in the ε v p -log p plane and

a(r k ) = a cyc + (a m -a cyc )α k , ( 17
)
a m , a cyc , b, r el k , p co , φ and α k are constant parameters of the model. In the k-plane, the plastic deviatoric strain rate tensor ėp k and the plastic volumetric strain rate (ε p v ) k are defined as:

ėp k 1 = εp ii k - εp jj k = λ p k d k 1 , ėp k 2 = 2 εp ij k = λ p k d k 2 , ( 18
)
εp v k = εp ii k + εp jj k = λ p k v k ,
where λ p k is the plastic multiplier and k gives the direction of plastic strains, assuming that this tensor can be decomposed into a deviatoric and a volumetric part.

We make the hypothesis that, in the normalised deviatoric plane, the law is associated, which defines d k :

d k = ∂ s k f k = s k s k . ( 19
)
The plastic potential relationship permitting contractancy and dilatancy is given by:

v k = dε p v k dε p d k I = α k sin φ - s k : d k p k I.
(20)

Isotropic mechanism

Experimental results exhibit plastic strains during an isotropic compression. Because they are not activated under this stress path, the deviatoric mechanisms cannot generate purely volumetric strains. An isotropic mechanism is thus introduced.

f 4 (p, p c , r 4 ) = | p| -r 4 , (21)
where

p = p F 4 (p c ) (22)
and F 4 (p c ) is a factor of normalisation given by:

F 4 (p c ) = dp c , ( 23
)
r 4 is the degree of mobilisation of this mechanism. Its value increases continuously from r el 4 (elastic) to 1 and its evolution law is as hyperbolic form given by:

ṙ4 = λ 4 p (1 -r 4 ) 2 c p a p c . ( 24
)
The mechanism generates plastic volumetric strains:

εp ij 4 = λ p 4 v 4 δ ij 3 . ( 25
)
c and d are constant parameters of the model.

During loading, the four mechanisms can be activated and are coupled by the density hardening p c as follows:

ε p v = 4 m=1 ε p v m . ( 26
)
The equation of compatibility, which permits us to determine the plastic multipliers, is given for each of these mechanisms by:

ḟ σ, ε p v , r = ∂ σ f : σ + ∂ ε p v f εp v + ∂ r f ṙ. ( 27
)
Several authors [START_REF] Lee | Drained strength characteristics of sands[END_REF], [START_REF] Vesic | Behavior of granular materials under high stresses[END_REF], [START_REF] Ramamurthy | Behavior of coarse-grained soils under high stresses[END_REF], Colliat-Dangus (1986), [START_REF] Lade | Significance of particle crushing in granular materials[END_REF], Yamamuro andLade (1996a, 1996b)) observed that the compressibility of granular medium increases with grain breakage. Thus, the dilatancy decreases and tends to vanish. [START_REF] Lee | Drained strength characteristics of sands[END_REF] have observed: "The stress-strain volume change characteristic of dense sand at high confining pressures are not unlike those of low pressures". We believe that the concept of critical state in sands can then be used (see also [START_REF] Been | The critical state of sands[END_REF].

The methodology for connecting the parameters of the discontinuous medium and the equivalent fictitious continuous medium has been applied in order to propose links, by means of correlations, between the parameters of the elastoplastic model and the nature and compacity parameters of granular media. The approach was developed by [START_REF] Hicher | Micro-macro correlations for granular media. Application to the modelling of sands[END_REF], who used a database called Modelisol [START_REF] Favre | Modelisol: a database for reliability analyses in geotechnique[END_REF] which contained experimental results on 21 different sands. They first assumed that the previous studies gave enough elements for the statistical determination of elastic parameters E, v, n and perfect plasticity (or critical state) parameters: φ, β, d, p co (Biarez and Favre, 1977;[START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF]. For p co , the initial density influence was introduced by means of the density index I D = e max -e 0 e max -e min taking into account the values of the two nature parameters e max and e min as well as the value of the initial void ratio e 0 .

They concentrated their study on the other parameters, specific to Hujeux's model: r el , a, b, α which appear in the constitutive equations ( 11) to ( 27). For this purpose, the principal component analysis method was chosen. Using the methodology of connections between parameters of discontinuous medium described above, the principal component analysis expresses the fact that the parameters describing the discontinuous medium can be presented in factorial space either by parameters describing the nature and the compacity of the granular material (FSPD) or by the parameters of the constitutive model (FSPC). We can say that the space FSPD is the solution, by means of the correlations and statistical links of the space FSPC. Thus the values of the constitutive model parameters are the result of the projection of FSPD on FSPC.

It was therefore possible to determine the variables which were strongly correlated in order to propose, by means of multiple regression, mathematical links between model parameters and nature parameters.

The following results were obtained: if Y i is a parameter of the elastoplastic model, Y i can be written:

Y i = g i (y j ), j = 1, . . . , n, ( 28
)
where y j are the parameters describing the properties of the granular medium (nature and compacity parameters). In fact, only E and p co are depending on the compacity, the other parameters are constant for a given medium and therefore, depend only on the nature parameters. At the present stage, these correlation models can be considered as very valuable, not for the final determination of the model parameters, but more likely for the definition of an initial set of parameters, realistic enough to permit an optimal optimisation process to be undertaken.

Accounting for grain breakage in the constitutive model

The first step was to propose a method for quantifying the amount of grain breakage. Several proposals which consider generally the evolution of different parameters representative of the grain size distribution have already been made. Having tested some of them [START_REF] Kim | Etude expérimentale du comportement mécanique des matériaux granulaires sous fortes contraintes[END_REF], we decided to use a parameter which could be representative of the whole evolution of the grain size distribution in agreement with [START_REF] Fukumoto | Particle breakage characteristics of granular soils[END_REF]. For this reason, we select the surface S intercepted by the initial and the actual grain size distribution curves (see figure 8).

As we saw above, the evolution of S occurred for two kind of parameters: nature and mechanical parameters. A convenient way of taking into account both stress and strain is to consider the plastic work:

W p = σ ij dε p ij . ( 29
)
Experimental results show that triaxial tests produce more grain breakage than isotropic or one-dimensional compression tests. In the first expression, as axial and radial strains have a different sign (the first is negative while the second is positive), the plastic work, for a given stress and axial strain, is smaller for a triaxial test than for an one-dimensional compression test. By introducing the absolute value of the strain increment, we can correct this point. So we decided to use the following expression:

W * p = σ ij dε p ij . ( 30
)
The whole surface S can be given by a function of W * p , which depends on nature parameters and can be written:

S = ψ ph (W * p ), ( 31
)
ψ ph is a function depending only on the initial granular material. With reference to Hujeux's model, we can now introduce this function to take into account the material changes due to grain breakage.

The equation ( 28) was proposed assuming that the nature parameters were constant along any stress and strain path. Should this hypothesis no longer be valid, as in the case of significant grain ruptures, the equations have to be modified in order to follow these changes.

We can assume that y j are modified by grain ruptures in the following manner:

y j = h j (S). ( 32
)
The model parameters can therefore be written:

Y i = g i h j (S) . ( 33
)
A study of different factors of influence has shown that the principal parameters to be taken into account were those which control the position of the perfect plasticity relation, that is for the model p co and β which control hardening in density by the relation:

p c = p co exp βε p v . ( 34
)
Let us reconsider the correlations in figure 4. We can observe that the grain ruptures, by provoking an increase of d 60 /d 10 , will cause a shift of the initial relation of perfect plasticity towards the smaller void ratios.

Under these conditions, the critical pressure p co corresponding to a reference void ratio e 0 diminishes while the slopes of the lines are almost equal. We therefore propose the hypothesis that the modulus of plastic compressibility β remains constant during the whole loading even if the physical and mechanical properties evolve as shown in figure 4.

The calculation of the value of p co obtained from experimental tests made on quartz sand as well as on carbonate sand shows that its diminution is correlated with that of calculated plastic work. The value of p co carried up to its maximum value varies between 1 and 0.

We can then write:

p co = p coi 1 -χ(S) . ( 35
)
The most appropriate form for the relation between p co and W * p is the hyperbolic form ( figure 9). The variable χ(S) has the following analytical form [START_REF] Hicher | Experimental evidence and modelling of grain breakage influence on mechanical behaviour of granular media[END_REF] and [START_REF] Daouadji | Modélisation de l'influence de la rupture des grains sur le comportement des matériaux granulaires[END_REF]:

χ(S) = 1 - p co p coi = W * p B + W * p . ( 36
)
Hardening in density will then be written as follows:

p c = p coi . 1 - W * p B + W * p . exp β.ε p v (37)
with p co : value of the critical reference pressure after loading (grain crushing may occur), p coi : value of the critical reference pressure before loading (without grain crushing), B: parameter taking into account the effect of grain ruptures which depend on the initial properties (shape, size, distribution. . .) of the grains.

Parameter B allows us to characterise different materials while taking into account the individual characteristics of grains (resistance, form. . .) and of their grain size distribution.

We have also taken into account the influence of the loading path as it has been stated above. The use of a multimechanism model that distinguishes between deviatoric and isotropic paths allows us to take into account this influence and support the choice of this model. We have therefore taken into consideration the plastic work corresponding to each of these mechanisms. Thus, if one mechanism is more activated than another, the plastic work corresponding to this mechanism is higher than that of other mechanisms. This is the case notably for the simulation of a triaxial test where the deviatoric mechanisms intervene in a more important way than does the isotropic mechanism.

In re-writing equation ( 37) for each mechanism, we obtain:

p ck = p coi . 1 - W * p k B k + W * p k . exp β.ε p v . ( 38
)
Parameter B k takes on different values for deviatoric mechanisms and for isotropic mechanism. We note that the coupling between different mechanisms is always provided by the plastic volume change. The plastic multipliers are obtained by resolving for each mechanism the equation ( 27).

Validation of the model

In order to validate our approach, we have simulated the behaviour of different kinds of sands [START_REF] Colliat-Dangus | Comportement des matériaux granulaires sous fortes contraintes. Influence de la nature minéralogique du matériau étudié[END_REF]). The first one is a quartz sand coming from the Hostun quarry in France. This is a fine sand of angular grains and of uniform grain size distribution. The second sand is of marine origin with a very high percentage of calcite (98%), and containing a great amount of shell debris. The grains are rounded and have an intra-granular porosity of 9%. The grain size distribution for the second sand is more spread out than it is for the first sand. Table II recapitulates the information concerning these two sands.

The individual resistance of the sand grains varies very much. It is calculated in running a simple compression test on a grain, in assimilating it to a sphere of diameter d. The resistance to rupture by indirect traction is given by Jaeger (1967), [START_REF] Billam | Some aspects of the behaviour of granular materials at high pressures[END_REF], [START_REF] Colliat-Dangus | Comportement des matériaux granulaires sous fortes contraintes. Influence de la nature minéralogique du matériau étudié[END_REF] and [START_REF] Dowell | On the micromechanics of crushable aggregates[END_REF]:

R t = γ * P r d 2 (39)
with P r : load provoking rupture, ( * * ) For an effective mean stress p' 0 = 2.2 MPa.

d: diameter of the sphere, γ : coefficient of proportionality dependent on the grain. These tests were made on ten grains of the same size and on all grain sizes. The values of resistance to indirect traction are presented in figure 10 as are the dispersions obtained for each grain size. The tests performed by [START_REF] Billam | Some aspects of the behaviour of granular materials at high pressures[END_REF] as well as those by [START_REF] Ramamurthy | Behavior of coarse-grained soils under high stresses[END_REF] are also presented. For example, the value of resistance to indirect traction for a diameter corresponding to d 50 is 20 MPa for the calcareous sand, whereas for the quartz sand, it is 200 MPa.

We have simulated laboratory tests on isotropic and triaxial paths in a large range of stresses. We have tried to reproduce the influence of grain ruptures which were produced during the isotropic phase as well as during the shearing phase. The determination of the parameters of the model was made from tests performed under small mean stresses. We have used the correlations proposed by [START_REF] Hicher | Micro-macro correlations for granular media. Application to the modelling of sands[END_REF] in order to determine a set of initial parameters which has been adjusted on triaxial test results at small stresses. All the parameters, without exception, were maintained constant for all the tests and are given in table III.

We can therefore take into account the effect of the grain rupture with the help of two parameters B k (deviatoric mechanism) and B 4 (isotropic mechanism). Figure 11 illustrates the influence of parameter B k on a drained triaxial path at lateral stresses equal to 1 MPa on the calcareous sand. We can note first of all that the set of parameters defined for the small stresses (50 KPa) and which characterises the initial material can no longer propose a good simulation for this level of mean stress. The evolution of the grain size distribution curve ( figure 12) shows effectively a significant change in the grain sizes. We also notice that an appropriate value of parameter B k (obtained by curve fitting) allows us to simulate well the experimental results for both the stress-strain relationship (q-ε 1 plane) and the volume change (ε 1 -ε v plane).

When we increase the amount of grain breakage, the peak in the stress-strain curve given by the model without grain breakage vanishes, the material becomes more ductile and the axial strain corresponding to the maximun value of the stress deviator increases. We can note in figure 11 that the dilatancy phenomenon vanishes with the increase of breakage. The volumetric strain reaches the value noted during the test. Let us point out that at the end of the simulation the volumetric strain increases from 2% (without breakage) to 10% (with breakage).

The same can be said for the isotropic stress paths thanks to the choice of the parameter B 4 . In figures 13 and 14 we present the isotropic consolidation curves up to a stress value of 15 MPa. The value of B 4 obtained for the quartzic sand and for the calcareous sand are respectively 1.5 MPa and 2 KPa. We notice that in the case of the Hostun sand very little breakage occurs during this phase whereas a great number of grains are crushed in the case of the calcareous sand (see the grain size distribution curves in figure 12). In the first case, the volumetric strain obtained is 8% while for the calcareous sand, this value rises to 30%.

This increase in the consolidation slope is here directly tied to the quantity of grain breakage obtained during loading. For the Hostun sand, its value is close to the values usually obtained for granular materials C c ≈ 0.1-0.2 [START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF]. On the other hand, for the calcareous sand, it is clearly higher C c = 0.48. This high value is a consequence of very great grain breakage [START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF] and [START_REF] Daouadji | Modelling of grain breakage influence on mechanical behavior of sands[END_REF]. We can state that given a well-chosen value B 4 , the modified model is capable of reproducing the experimental results.

It is therefore necessary to simulate the isotropic consolidation phase preceding shearing since there may be the presence of significant grain breakage and the triaxial tests cannot be simulated in a satisfactory manner without taking into account this phase.

We present in figure 15 the experimental results and the simulations of four triaxial tests performed on calcareous sand with lateral stresses between 50 KPa and 1 MPa (B k = 1.5 KPa). We can state that the modified model is capable of taking into account the behaviour of granular materials at both small stresses where it is strongly dilatant as well as at higher stresses under which the materials become very contractant due to the effect of numerous grain ruptures provoking a significant variation of the grain size distribution.

Simulations of fine Hostun sand are shown in figure 16. Simulations without breakage lead to results in the q-ε 1 plane which approximate the experimental curves. On the other hand, values obtained in the ε v -ε 1 plane for simulations under confining pressure greater than 3 MPa are much smaller than values obtained during testing. For a confining pressure under 2 MPa, little grain breakage occurs as we can see from the evolution of the percentage of fine material ( figure 12). With grain breakage, we can simulate more correctly triaxial tests up to 15 MPa. This test is presented in figure 17 with (B k = 0.7 MPa) and without grain breakage. In the case of grain breakage, we observe that for 20% of axial strain, we can reduce the deviatoric stress for more than 5 MPa and increase significantly the volumetric strain. Contrary to tests on calcareous sand which require little energy in order to induce grain breakage, tests on fine Hostun sand show that several MPa in confinement are needed in order to induce the noticeable influence of grain breakage.

Conclusion

In establishing the correlations between the characteristics of the discontinuous medium and the mechanical behaviour of the equivalent continuous medium, we find it possible to study the role of grain breakage, particularly of an elasto-fragile type breakage.

The parameter which best explains this phenomenon is the evolution of the grain size distribution represented by the value of the coefficient of uniformity d 60 /d 10 , which exerts a strong influence on the position of the isotropic consolidation line and of the perfect plastic one in the e-Logp plane. In particular, an increase of the coefficient of uniformity corresponds to a slipping of these relations towards the small void ratios. Grain breakage, in provoking an increase of d 60 /d 10 , consequently increases, on the one hand, the compressibility of granular material on a consolidation path (isotropic or one-dimensional compression), and on the other hand, the contractancy on a deviatoric path. This increase in contractancy is accompanied by an increase in ductility.

The grain breakage influence has been introduced in an elastoplastic model by making the critical state line dependent on the evolution of the grain size distribution. This evolution is computed by means of a function of the plastic work, which depends on the initial nature parameters of the granular material.

In modifying the model in this manner, we may simulate the behaviour of granular materials under a large range of stresses when there is an evolution of the physical and mechanical characteristics of grains.

A primary validation of this approach has been done on two different types of sands, a quartzic and calcareous one and on various loading paths with different stress intensities. The results we obtained reproduced with a sufficient accuracy the results of the laboratory tests.

Figure 1 .

 1 Figure 1. Different modes of grain breakage.

Figure 2 .

 2 Figure 2. Grain ruptures in one-dimensional compression tests for poorly graded sands.

Figure 3 .

 3 Figure 3. Evolution of the minimal and the maximal void ratio as a function of grain size distribution and the angularity of the grains.

Figure 4 .

 4 Figure 4. Perfect plasticity relation as a function of the grain size distribution.

Figure 5 .

 5 Figure 5. Grain breakage influence on stress-strain relationship on Hostun sand and crushed granite in comparison with standard behaviours.

Figure 6 .

 6 Figure 6. Evolution of the grain size distribution curve in a one-dimensional compression test. Influence of the initial distribution.

Figure 7 .

 7 Figure 7. Influence of stress and strain paths on grain breakage.

Figure 8 .

 8 Figure 8. Definition of the surface S.

Figure 9 .

 9 Figure 9. Correlation between the critical reference pressure and plastic work for Hostun sand and a calcareous sand.

Figure 10 .

 10 Figure 10. Tensile strength of individual particles as a function of their diameter.

  For an effective mean stress p' 0 = 0.05 MPa.

Figure 11 .

 11 Figure 11. Influence of grain breakage amplitude in a numerical simulation of a triaxial test on calcareous sand at p' 0 = 1 MPa.

Figure 12 .

 12 Figure 12. Grain size distribution curves before and after tests on Hostun and calcareous sands.

Figure 13 .

 13 Figure 13. Isotropic test on Hostun sand. Comparison between numerical and experimental results (p' up to 15 MPa).

Figure 14 .

 14 Figure 14. Isotropic test on calcareous sand. Comparison between numerical and experimental results (p' up to 15 MPa).

Figure 15 .

 15 Figure 15. Triaxial test on calcareous sand. Comparison between numerical and experimental results (0.05 p' 0 1 MPa).

Figure 16 .

 16 Figure 16. Triaxial test on Hostun sand. Comparison between numerical and experimental results (3 p' 0 5 MPa).

Table I .

 I Characteristics of materials used to determine a network (without grain breakage) compared to tests under elevated stresses (with grain breakage). : e 0 = 0.890 and I D = 32% for σ 3 = 300 KPa. * * : e 0 > e max as w = 5%.

		Le Long	Colliat-Dangus	Ziani	Kim
	Nature	Quartzic sand	Quartzic sand	Quartzic sand	Granite
		(Hostun)	(Hostun)	(Hostun)	G1	G2
	C u = d 60 /d 10	1.62	1.70	/	10	2
	e min	0.580	0.656	0.620	0.430	0.430
	e max	0.866	1.000	0.840	0.800	1.230
	e 0	0.680	0.673	1.190	0.467	0.870
	I D (%)	65	95 *	(very loose) * *	90	90

*

Table II .

 II Characteristics concerning the Hostun's sand and the calcareous sand.

	Sand	Mineralogy	e min	e max	D r (%)	C u	d 50 (mm)	Shape	R t (MPa)
	Hostun	Quartzic	0.656	1.00	90	1.7	0.32	angular	200
	Calcareous	Calcareous	1.014	1.670	71	2.8	0.17	rounded	20

Table III .

 III Values of the parameters of the model for the calcareous sand and the quartzic Hostun's sand.

	Value of the parameters	Calcareous sand	Hostun's sand
	Elastic parameters
	K	215	354
	G	165	265
	N	0.5	0.5
	Parameters of the critical state
	φ	39	31
	β	30	10
	p co	1.5 ( * )	2.2 ( * * )

Figure 17. Triaxial test on Hostun sand. Comparison between numerical and experimental results (p' 0 = 15 MPa).


Philippe Lorong

David Ryckelynck

J Yvonnet

Ph F Chinesta

D Lorong

Ryckelynck

The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces

Keywords: Natural neighbor interpolation, C-NEM, Constrained Voronoi diagram, Moving interface, Stefan problem

come

Introduction

Phase boundaries represent material interfaces across which several fields may exhibit sharp gradients, and even discontinuities. A wide range of numerical methods have been developed for treating these problems according to the pertinent physics and assumptions about the interface [START_REF] Lewis | Finite element simulation of metal casting[END_REF]. When a sharp interface is considered, its motion is governed by the jump in the temperature gradient normal to the phase boundary and is accompanied by latent heat effects (Stefan condition). In order to satisfy these conditions the most common approach lies in explicitly tracking the interface motion. Within the interface tracking approach two main alternatives exist: the moving mesh methods and the mixed Eulerian-Lagrangian methods. Moving finite element mesh methods conform element boundaries to the interface as it evolves. Although these methods are very accurate, they are limited by severe mesh distortion. Thus, frequent remeshing is needed, with the associated field projections between successive meshes. Moreover, remeshing is, even today, a delicate task in three dimensions. To alleviate remeshing efforts a number of Eulerian-Lagrangian methods have been developed recently that track the interface while solving the equations on a fixed grid [START_REF] Udaykumar | Computation of solid-liquid phase fronts in the sharp interface limit on fixed grids[END_REF]. Many of these methods effectively smear the discontinuity over a few grid cells, and are therefore not capable of representing the true discontinuity across the interface.

A new approach for representing localized behaviours has recently emerged in the field of the finite element method, known as the partition of unity method [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF]. The main idea is to extend the classical approximation considering the product of the standard shape functions and local enrichment functions. The extended finite element method (X-FEM) is a variation on this technique. Recently, the X-FEM has been coupled to the Level Set Method [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flows[END_REF] to represent interface [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF]. In this way, the discontinuity evolution can be properly represented on a fixed background mesh, just by adding an appropriate enrichment in the functional approximation in the elements that are intersected by the moving discontinuity [START_REF] Ji | A hybrid finite element/level set method for modeling phase transformations[END_REF]. However, when the material in which the interface is moving, is subjected to large displacements, an When the background mesh evolves, remeshing will also be required to avoid large mesh distortions. In the context of moving phase boundary problems, it has been shown in [START_REF] Ji | A hybrid finite element/level set method for modeling phase transformations[END_REF] that due to the C 0 finite element shape functions, errors are introduced when the moving interface crosses the edge of the elements.

To alleviate dependence on the mesh and to provide smoother shape functions, the use of meshfree or meshless methods has been investigated. Meshless methods discretize a continuum body by a finite number of particles (or nodes) and the field of interest is interpolated under these nodes without the aid of an explicit mesh. Many meshless methods have been proposed, including the smooth particle hydrodynamics (SPH) [START_REF] Lucy | A numerical approach to the testing of fusion process[END_REF], the Radial Basis Function methods [START_REF] Kansa | Multiquadratics-a scattered data approximation scheme with applications to computational fluid dynamics. I-Surface approximations and partial derivative estimates[END_REF][START_REF] Kansa | Multiquadratics-a scattered data approximation scheme with applications to computational fluid dynamics. II-Solutions to hyperbolic, parabolic and elliptic partial differential equation[END_REF] the diffuse element method (DEM) [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF], the element free Galerkin (EFG) [START_REF] Belytschko | Element-free Galerkin methods, Internat[END_REF], the reproducing kernel particle method (RKPM) [START_REF] Liu | Reproducing kernel particle methods[END_REF], the HP clouds [START_REF] Duarte | An H-p adaptative method using clouds, Comput[END_REF] and the partition of unity method (PUM) [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF].

The introduction of moving discontinuities in these meshless methods can present difficulties for the following reasons:

(i) The quality of the approximation as well as the conditioning of the global system is pathologically dependent on the size of the support of the shape functions; (ii) Imposition of essential boundary conditions needs special treatment; (iii) Integration is not accurate enough; and (iv) The physical discontinuity across the interfaces must be introduced accurately.

In order to overcome these different problems, we propose the use of the constrained natural element method (C-NEM) [START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi constraint-Méthode des éléments C-Naturels[END_REF][START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM)[END_REF] for treating thermal models involving moving interfaces. This approach is an extension of the natural element method [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF] in which both trial and test functions are constructed on the basis of the Voronoi based interpolants [START_REF] Hiyoshi | Improving continuity of Voronoi-based interpolation over Delaunay spheres[END_REF][START_REF] Sibson | A vector identity for the Dirichlet tessellations[END_REF]. These interpolants satisfy the Kronecker delta property and their support is defined by the union of the Delaunay spheres passing through the visible nodes. In the C-NEM, the introduction of a visibility criterion and its related constrained Voronoi diagram preserves the appealing properties of the NEM in any geometry (convex or not) and allows the introduction of material discontinuities.

Furthermore, the C-NEM is strongly related to the finite element method through the following points: (a) an underlying structure is needed (the constrained Delaunay triangulation/constrained Voronoi diagram); (b) the external boundary conditions can be imposed directly, due to the interpolant character and the strict linearity of the shape functions along any kind of boundary or interface. Nevertheless, in the C-NEM, both interpolation and integration are not constructed on the basis of the Delaunay triangles, but on its dual, the Voronoi diagram (which can be extended to constrained Delaunay triangulation [START_REF] Seidel | Constrained Delaunay triangulations and Voronoi diagrams with obstacles[END_REF] for non-convex domains or when they involve fixed or moving discontinuity), which provides smoother approximation and allows to proceed on the basis of very distorted Delaunay triangles. In a recent paper [START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM)[END_REF] the C-NEM has been successfully applied in some problems involving non-convex domains and cracks discontinuities [START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM)[END_REF].

The layout of this paper is as follow: in Section 2, we introduce a simple mathematical model of a thermal problem involving a moving interface (Stefan problem [START_REF] Vuik | Some Historical Notes about the Stefan Problem[END_REF]). In Section 3, the C-NEM is summarized, and it will be applied in Section 4 to discretize the weak formulation of the nonlinear Stefan problem. A numerical benchmark is presented in Section 5, which allows us to test the accuracy of the proposed technique by comparing the numerical results with the exact solution of the problem.

Problem formulation

Let Ω ∈ 2 be a bounded domain and T the temperature field. On the domain boundary Γ ≡ ∂Ω the temperature or the thermal flux is prescribed. We will denote by Γ 1 the domain boundary where the temperature is known T (x ∈ Γ 1 , t) = T (x, t) and by Γ 2 the domain boundary where the heat flux q is imposed. The thermal model is defined in the time interval [0, t max ]. The initial temperature T (x, t = 0) = T 0 , where T 0 is assumed to be higher than the material melting temperature T m . At time, t = 0, a part of the domain boundary Γ 1 is suddenly exposed to a temperature T 1 < T m . A moving solidification front Γ I is then generated, whose position evolves in time, i.e., Γ I (t), dividing the domain Ω in two regions Ω 1 (t) (containing the solid phase at time t ) and Ω 2 (t) (which contains the liquid phase) as shown in Fig. 1. For the sake of simplicity we will consider, from now on, a homogeneous and isotropic thermal model in both phases.

The heat transfer model is defined in each phase, neglecting volumetric source terms as well as the motion in the liquid phase induced by thermal gradients, by:

c 1 ∂T (x,t) ∂t = ∇ • (k 1 ∇T ) in Ω 1 (t) c 2 ∂T (x,t) ∂t = ∇ • (k 2 ∇T ) in Ω 2 (t) (1)
where c 1 and c 2 are the volumetric heat capacities of both phases, and k 1 and k 2 are the respective thermal conductivities. The associated initial and boundary conditions are:

T (x, t = 0) = T 0 ∀x ∈ Ω T (x, t) = T (x, t) ∀x ∈ Γ 1 , ∀t ∈ [0, t max ] -k∇T (x, t) • n = q(x, t) ∀x ∈ Γ 2 , ∀t ∈ [0, t max ] (2)
The evolution of the interface Γ I (t) is described by a Stefan condition:

V x ∈ Γ I (t) = |[q]| L n I (x) ( 3
)
where V is the interface velocity, L is the volumetric latent heat of fusion, n I (x) is the normal vector to the interface at point x which is assumed to point into the liquid phase, and |[q]| the thermal flux jump across the interface Γ I (t), i.e.,

|[q]| = k 1 ∇T Γ - I (t) -k 2 ∇T Γ + I (t) n I (4)
The additional constraint prescribed on the interface Γ I (t) is:

T (x, t) = T m ; ∀x ∈ Γ I (t) (5)
where T m is the melting temperature.

The constrained natural element method (C-NEM)

In this section, the utility of the C-NEM to describe moving interfaces and discontinuities in a fixed cloud of nodes is discussed. After a brief review of the Voronoi-based interpolants, we introduce the constrained Voronoi diagram which is used for computing the shape functions in any domain.

Natural neighbor interpolation

We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. For a more in-depth discussion on the Sibson interpolant and its application for solving secondorder partial differential equations, the interested reader can refer to Sambridge et al. [START_REF] Sambridge | Geophysical parameterization and interpolation of irregular data using natural neighbors[END_REF], and Sukumar et al. [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. The NEM interpolant is constructed on the basis of the Voronoi diagram. The Delaunay tessellation is the topological dual of the Voronoi diagram.

Consider a set of nodes S = {n 1 , n 2 , . . . , n N } in 2 . The Voronoi diagram is the subdivision of 2 into regions T i (Voronoi cells) defined by:

T i = x ∈ 2 : d(x, x i ) < d(x, x j ), ∀j = i , ∀i (6)
The Sibson coordinates of x with respect to a natural neighbor n i (see Fig. 2) is defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area (volume in 3D) of the Voronoi cell related to point x:

φ i (x) = Area(af ghe) Area(abcde) (7)
If the point x coincides with the node n i , i.e., x = x i , φ i (x i ) = 1, and all other shape functions are zero, i.e., φ j (x i ) = δ ij (δ ij being the Kronecker delta). The properties of positivity, interpolation, and partition of unity are then verified [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]:

   0 φ i (x) 1 φ i (x j ) = δ ij n i=1 φ i (x) = 1 (8)
The natural neighbor shape functions also satisfy the local coordinate property [START_REF] Sibson | A vector identity for the Dirichlet tessellations[END_REF], namely:

x = n i=1 φ i (x)x i (9)
which combined with Eq. ( 8) implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness). Natural neighbor shape functions are C ∞ at any point except at the nodes, where they are only C 0 , and on the boundary of the Delaunay circles (spheres in 3D) where they are only C 1 , because of the discontinuity in the neighbors nodes across these boundaries. Hiyoshi and Sugihara [START_REF] Hiyoshi | Improving continuity of Voronoi-based interpolation over Delaunay spheres[END_REF] have shown that the Sibson interpolant belongs to a more general class of Voronoi-based interpolants, called kth order standard coordinates, proving that the interpolant generated by the kth order standard coordinates have C k continuity on the Delaunay circles (spheres) boundaries. In this context, the Sibsonian and non-Sibsonian (Laplace) coordinates [START_REF] Belikov | The non-Sibsonian interpolation: a new method of interpolation of the values of a function on an arbitrary set of points[END_REF] results to be the standard coordinates of order 1 and 0, respectively.

Another important property of this interpolant is the ability to reproduce linear functions over the boundary of convex domains. The proof can be found in Sukumar et al. [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. An illustration is depicted in Fig. 2(b): as the areas associated to points on the boundary become infinite, the contribution of internal points vanish in the limit when the point approaches the convex boundary, and the shape functions associated with nodes n 1 and n 2 become linear on the segment (n 1n 2 ). This is not true in the case of non-convex boundaries, and the next section focuses on an approach to circumvent this difficulty.

Consider an interpolation scheme for a scalar function T (x) : Ω ⊂ 2 → , in the form:

T h (x) = n i=1 φ i (x)T i ( 10
)
where T i are the nodal temperatures at the n natural neighbor nodes, and φ i (x) are the shape functions associated with each neighbor node. It is noted that Eq. ( 10) defines a local interpolation scheme. Thus, the trial and test functions used in the discretization of the variational formulation describing the thermal problem treated in this paper take the form of Eq. (10).

The constrained natural element method

Constrained Voronoi diagram

It was proved in [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF][START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM)[END_REF] and [START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbor Galerkin methods[END_REF] that spurious influences between "non-visible" nodes and loss of linearity in the interpolation along boundaries of non-convex domains appear in the framework of the NEM. In order to avoid this drawback and to recover all properties of the method for any geometry (including non-convex domains containing cracks or involving field discontinuities), a visibility criterion is introduced in order to restrict influent nodes among natural neighbors. The computation of the shape functions is done on the basis of the so-called constrained (or extended) Voronoi diagram (CVD), introduced by Seidel in [START_REF] Seidel | Constrained Delaunay triangulations and Voronoi diagrams with obstacles[END_REF]. Imagine 2 is a sheet of paper Σ 0 , with the points of the set of nodes S and the line segments defining the boundary in a set L, drawn on it. For each l i ∈ L, we cut Σ 0 open along l i and glue another sheet Σ i , which also cut open along l i . The gluing is done around l i such that every traveler who crosses l i switches from Σ 0 to Σ i and vice versa. A schematic view of the particular gluing necessary to achieve that effect is illustrated in Fig. 3(a). We know what it means for two points on the primary sheet to be visible from each other. For other pairs we need a more general definition. For i = 0, points x 0 ∈ Σ 0 and y i ∈ Σ i are visible if xy crosses l i , and l i is the first constraining line segment crossed if we traverse xy in the direction from x to y.

In Fig. 3(b), the intersection between the CVD and the domain closure is depicted. The resulting diagram is composed of cells T C i , one for each node n i , such that any point x inside T C i is closer to n i than to any other node n j visible from x. We call this the constrained Voronoi cells, which are defined formally by: where Γ is the domain boundary, composed with segments l i ∈ L and S a→b denotes the segment between the points a and b.

T C i = x ∈ n : d(x, x i ) < d(x, x j ), ∀j = i, S x→n i ∩ Γ = ∅, S x→n j ∩ Γ = ∅ (11)
The constrained Delaunay triangulation does not always exist in 3D without adding new nodes, as shown in [START_REF] Schönhardt | Über die Zerlegung von Dreieckspolyedern in Tetraeder[END_REF]. Nevertheless, some techniques for constructing 3D constrained Delaunay tessellations are available and provided in [START_REF] Shewchuck | Tetrahedral mesh generation by Delaunay refinement[END_REF][START_REF] Shewchuck | Sweep algorithms for constructing higherdimensional constrained Delaunay triangulations[END_REF] by addition of Steiner points.

The constrained natural element approximation

In order to solve partial differential equations defined on non convex domains, or to reproduce functional discontinuities, we consider the following approximation for both the trial and the test functions:

T h (x) = V i=1 φ C i (x)T i ( 12
)
where V is the number of natural neighbors visible from point x and φ C i is the constrained natural neighbor shape function related to the ith node at point x. The computation of the C-n-n (constrained natural neighbor) shape functions is similar to the natural neighbor shape function, when one proceed using the constrained Voronoi diagram introduced previously. It was shown in [START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi constraint-Méthode des éléments C-Naturels[END_REF] and [START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM)[END_REF] that the use of the constrained Voronoi diagram does not affect the properties of the NEM interpolation, allowing the extension of the linearity of the shape functions on the convex domains boundaries, to any geometry, convex or not.

The ability of the C-NEM for treating problems involving cracks has been illustrated in [START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM)[END_REF]. In the present paper, we focus on its application in the context of a moving interface defining two domains with different thermal properties. Thus, defining at time t two CVD (constrained Voronoi diagrams) of Ω 1 (t) and Ω 2 (t), both with respect to the interface Γ I (t), it can be proved that the interpolated temperature field is C 1 everywhere, except at the nodes and on the interface Γ I (t) where it is only C 0 . Thus, this interpolation seems to be appropriate to simulate the Stefan problem considered in this paper.

To illustrate this behavior, we consider the situation depicted in Fig. 4, where the point x moves from Ω 1 to Ω 2 . If x is in Ω 1 , the interpolated field is constructed from Eq. ( 12) using the neighbor visible nodes from point x (Γ I is assumed opaque). If x is on Γ I , according to the previous discussion, the interpolated field is strictly linear because it only depends on the two neighbor nodes located on Γ I . Finally, when x is in Ω 2 , the interpolated field is defined using the visible neighbor and visible nodes from point x (Γ I being opaque). The continuity of the interpolated field is then guaranteed, but a discontinuity appears in the field derivatives, because of a sudden change in the neighbor nodes across the interface. We can then reproduce the temperature field continuity, as well as the expected flux discontinuity on the interface.

C-NEM discretization

The weak formulation associated with Eq. ( 1) results:

Find T ∈ H 1 (Ω) with T = T on Γ 1 such that: Ω c ∂T ∂t δT dΩ = - Ω k∇T • ∇ δT dΩ + Γ I (t) |[q]| δT dΓ ∀δT ∈ H 1 0 (Ω) ( 13
)
where H 1 (Ω) and H 1 0 (Ω) are the usual Sobolev functional spaces. Substituting the trial and test functions (both approximated in the C-NEM framework) in the above equation and using the arbitrariness of the field δT , the following system of equations is obtained:

C Ṫ + KT = F ( 14
)
where T is the vector containing the unknown nodal temperatures. We consider the solution on the time interval [0, t max ], partitioned into steps as [t n , t n+1 ] and the generalized trapezoidal time stepping algorithm characterized by the parameter α:

∂T n+1 ∂t = T n+1 -T n -(1 -α) t ∂T n ∂t α t ( 15
)
which leads to:

(C n+1 + α tK n+1 )T n+1 = F n+1 T n , |[q]| n+1 ( 16
)
with:

C n+1 = Ω n+1
N t cN dΩ [START_REF] Merle | Solving thermal and phase change with the extended finite element method[END_REF] where

Ω n+1 = Ω n+1 1 ∪ Ω n+1 2 ; K n+1 = Ω n+1
B t kB dΩ [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF] and

F n+1 = C n+1 T n + (1 -α) t Ω n+1 N t c ∂T n ∂t dΩ + α t Γ n+1 I N t |[q]| n+1 dΓ ( 19
)
where N is the vector containing the nodal shape functions:

N = { φ 1 φ 2 . . . φ N }
and B is the matrix containing the shape functions derivatives:

B = φ 1,x φ 2,x . . . φ N,x φ 1,y φ 2,y . . .

φ N,y

The stabilized conforming nodal integration proposed by Chen et al. in [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] is employed for the numerical integration of K (see our former work [START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM)[END_REF] for more details). A lumped mass matrix C is computed making use of the constrained Voronoi cells areas as nodal weights.

The iteration procedure is defined as: Knowing T n and |[q]| n at time t n , the non-linear problem associated with Eq. ( 16) results in finding T n+1 and |[q]| n+1 such that Eqs. [START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbor Galerkin methods[END_REF] and [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] are satisfied. For this purpose we proceed as follows:

(1) Compute the interface velocity V n (x) using Eq. ( 3) and update the interface position at time t n+1 using the forward Euler formula:

x n+1 J = x n J + tV n x n J ( 20
)
where x J are the nodes defining the interface. (2) Update locally the constrained Voronoi diagram and the shape functions associated with integration points in the interface neighborhood. Then, we compute C n+1 and K n+1 . (3) Solve Eq. ( 16) using a Newton-Raphson procedure where the tangent matrix is computed numerically. (4) Repeat while t n+1 < t max .

An alternative scheme using the Latin method [START_REF] Ladevèze | Non-Linear Computational Structural Mechanics[END_REF] in the extended finite element framework can be found in Merle and Dolbow [START_REF] Merle | Solving thermal and phase change with the extended finite element method[END_REF].

Numerical examples

Unidirectional solidification of a semi-infinite solid

In this section, we illustrate the potential of the proposed technique in simulating a two-phase Stefan problem. The problem is essentially one-dimensional, but we solve it here in two dimensions to underscore the appealing features of the method.

The Stefan problem models the one-dimensional freezing of a semi-infinite domain (x 0). The initial temperature T 0 is assumed constant in the whole domain, being higher than the melting temperature T m . At time t = 0 the temperature at the left boundary x = 0 is suddenly prescribed to a value T 1 lower than the melting point, originating a solidification front that progresses from the boundary x = 0 in the x direction. The exact flow front position x f (t) is given by:

x f (t) = 2λ β s t ( 21
)
where β s = k s /c s is the thermal diffusivity of the solid phase, and the constant λ satisfies the following relationship:

e -λ 2 erf (λ) = k l √ η(T 0 -T m )e -ηλ 2 k s (T m -T 1 )erfc(λ √ η ) + λL √ π c s (T m -T 1 ) (22)
with η = β s /β l being the ratio of the thermal diffusivities and where k l represents the liquid phase conductivity. The temperature field in the solid phase 0 x x f (t) is then:

T (x, t) = T 1 + T m -T 1 erf (λ) erf x 2 √ β s t (23)
and in the liquid phase x x f (t):

T (x, t) = T 0 - T 0 -T m erfc(λ √ η) erfc x 2 √ β l t (24)
where erf and erfc are the error function and complementary error function, respectively.

In the present investigation, we use the water-saturated sand thermal properties provided in [START_REF] Lynch | Continuously deforming finite elements for the solution of parabolic problems, with and without phase change[END_REF] that are listed in Table 1. T 1 and T 0 were set to 263.15 K and 277.15 K, respectively (λ = 0.3073). We simulate the evolution of the temperature field in Ω = [0, 0.01] × [0, 0.005] m. In order to use the infinite domain solution as reference solution, the temperature at x = 0.01 m is prescribed to its expected value according to Eq. [START_REF] Sibson | A vector identity for the Dirichlet tessellations[END_REF].

In a first test, we consider in the domain Ω a 20 × 10 uniform grid and t = 2 s. In the following, a comparison between C-NEM and FEM has been performed. For this purpose, the FE shape functions have been computed on the basis of the constrained Delaunay triangulation, that can be obtained by connecting the natural neighbors in the constrained Voronoi diagram. Fig. 5 compares the computed interface position and the exact one. Excellent accuracy can be noticed for the C-NEM solution, as depicted in Fig. 6, where the error in the front position is represented. Fig. 7 shows the temperature profile at different times. The C-NEM accuracy is found greater than the FEM solution related to the same Delaunay triangles. In order to assess the thermal flux jump across the interface, the computed solution |[q h (t)]| is compared to the analytical solution |[q ex (t)]| given by: in the temperature gradient is accurately obtained, being in excellent agreement with the exact solution.

q ex (t) = Lλβ s √ β s t (25)
In a second test, we consider the domain Ω containing 200 nodes distributed at random. The purpose of the present test is to investigate the meshless feature of the technique, in which due to its meshless character no geometrical restrictions concerning the relative nodal positions are involved. Thus, neither the background nodal distribution nor the relative position of the nodes defining the moving interface with respect to the background nodes, induce a lack of accuracy when high distortions in the Delaunay mesh, used to compute the Voronoi diagram, takes place. This is the main difference between the proposed strategy and the standard finite element method whose accuracy depends significantly on the geometrical quality of the mesh. Moreover, this test approaches the situations encountered when the material is also moving, inducing highly irregular nodal densities and high background mesh distortions. Fig. 10 depicts the cloud of nodes and the interface position as well as the associated constrained Voronoi cells. Remarkably, despite the very ir-regular nodal distribution and density, we can notice that the interface computed in the C-NEM simulation is not distorted as it moves through the domain, unlike the FEM solution, from which we can conclude that the accuracy in the C-NEM is not significantly affected by the regularity of the nodal distribution. In Fig. 7 some temperature profiles along the line y = 0.0025 m are depicted, from which excellent accuracy can be inferred for the C-NEM solution.

In order to examine the convergence behaviour of the present approach, we conduct a series of calculations on increasingly refined meshes. The results for the relative error in the average front position for t = 50 s are provided in Fig. 11, which proves its convergence.

Solidification from a line heat sink

In this problem, a domain Ω = [-0.01, 0.01] × [-0.01, 0.01] m initially entirely in a liquid state (T 0 > T m ) is exposed to a continuous line heat sink located at the point (0, 0), originating the nucleation of an axisymmetric phase boundary separating the liquid and the solid phases whose radius increases in time. This problem was investigated by Ji et al. in [START_REF] Ji | A hybrid finite element/level set method for modeling phase transformations[END_REF] in the context of the X-FEM method. The analytical solution for this problem defined in the unbounded 2 can be found in [START_REF] Carslaw | Conduction of Heat in Solids[END_REF]. The radius of the front is given by: R f = 2λ β s t (26) The temperature in the solid region r < R f is given by:

T (r, t) = T m + Q 4πk s Ei - r 2 4β s t -Ei -λ 2 (27)
and in the liquid region r > R f by:

T (r, t) = T 0 - T 0 -T m Ei(λ 2 η) Ei - r 2 4β l t ( 28
)
In the above, λ is the root of the equation:

Q 4π e λ 2 = λ 2 β s L - k l (T 0 -T m ) Ei(-λ 2 η) ( 29
)
with Q the heat sink intensity, Ei the exponential integral and η the ratio of thermal diffusivities:

η = β s β l ( 30
)
In our problem, only the upper quadrant Ω = [0, 0.01] × [0, 0.01] of the domain is modeled, due to the symmetry. By imposing the exact solution on the external boundaries, we may model an axisymmetric problem in a square domain. It is noteworthy that the above solution is singular for r = 0.

To overcome this difficulty, we move the node located at the origin to (+h/2, +h/2), h being the nodal distance between Simulations with regular and randomly distributed nodes are performed. In the first simulation, the domain contains 144 nodes on a uniform grid, and 8 additional nodes on the interface, which match the exact solution for t 0 = 1 s. Fig. 12 shows the evolution of the Delaunay triangles depicted at times t = 1.15 s, t = 2.65 s and t = 5.5 s. The triangulation is constrained by the moving interface. The dual constrained Voronoi cells are depicted at time t = 9.25 s.

In Fig. 13, the computed front is compared with the analytical solution. Higher accuracy is obtained using C-NEM in comparison with the FEM solution. Thus, an excel- lent agreement with the analytical solution is obtained (see Figs. 14 and15).

Conclusion

In this paper, the salient features of the C-NEM method are used for treating thermal problems involving moving interfaces. In the C-NEM framework, the Delaunay triangulation is used for its convenience in the frequent connectivity update during the motion of the interface. The use of a constrained Delaunay triangulation allows to conform the boundary with the interface and to accurately reproduce the jump in the heat flux across the interface. The traditional C 0 finite element shape functions are replaced in the Galerkin scheme by the smoother natural neighbor (Sibson) shape functions, which are constructed on the basis of the dual to the Constrained Voronoi diagram (dual to the constrained Delaunay triangulation). A nodal integration is also performed on the basis of the Voronoi cells by using this dual structure instead of the triangles structure. This technique allows the use of very distorted Delaunay triangles without loss of accuracy, and avoids the numerical problems associated with the discontinuity across the interface. The enrichment of the natural neighbor interpolation with discontinuous shape functions related to a level-set description in the context of the partition of unity is a work in progress.

Fig. 1 .

 1 Fig. 1. Two phases problem.

Fig. 2 .

 2 Fig. 2. Construction of the Sibson shape functions.

Fig. 3 .

 3 Fig. 3. Interpretation of the constrained Voronoi diagram. (a) Schematic view of the CVD. (b) Constrained Delaunay triangulation with respect to segment l i and intersection between the dual CVD and the domain closure.

Fig. 4 .

 4 Fig. 4. Reproducing discontinuous derivatives using the constrained Voronoi diagram.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Computed FEM and C-NEM front position versus the exact solution using a 20 × 10 regular grid.

Fig. 8

 8 Fig.8shows a comparison between the FEM and C-NEM solution for Ω = [0, 0.01] × [0, 0.005] m containing 200 nodes, uniformly or randomly distributed. we can notice from Fig.8that C-NEM solution provides a higher accuracy compared to the FEM solution. Fig.9depicts the error between the C-NEM computed flux jump and the exact solution[START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF] for different refined meshes, proving the convergence of the approach. We point out that the discontinuity

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Temperature profiles along the line y = 0.0025 m using a regular grid.

Fig. 10 .

 10 Fig. 10. Computed interface position using an irregular cloud of nodes: (a) Cloud of nodes and interface position; (b) Constrained Voronoi cells.

Fig. 11 .

 11 Fig. 11. Relative error in the average front position for t = 50 s.

Fig. 12 .

 12 Fig. 12. Evolution of the constrained Delaunay triangulation at times t = 1.15 s, t = 2.65 s and t = 5.5 s, and constrained Voronoi cells at time t = 9.25 s, using an unstructured nodal distribution.

Fig. 13 .

 13 Fig. 13. Comparison between computed and exact analytical interface position: (a) t = 1.15 s, (b) t = 2.8 s, (c) t = 4.3 s, (d) t = 5.65 s, (e) t = 7.15 s, (f) t = 8.65 s.

Fig. 14 .

 14 Fig. 14. Interface radius position.

Fig. 15 .

 15 Fig. 15. Error in the interface location prediction.

Table 1

 1 Thermal properties of the water saturated sand

	Properties	Solid	Liquid
	Volumetric heat capacity (kg•m -1 •s -2 •K -1 )	2.051 × 10 6 2.595 × 10 6
	Thermal conductivity	4.019	2.888
	(kg•m -1 •s -3 •K -1 )
	Melting temperature (K)		273.15
	Volumetric latent heat of fusion
	(kg•m -1 •s -2 )		8 .038 × 10 7