Moez Smaoui
email: smaouimoez@voila.fr

Zoubeir Bouaziz
email: zoubeir.bouaziz@enis.rnu.tn

Ali Zghal
email: ali.zghal@esstt.rnu.tn

Maher Baili
email: maher.baili@enit.fr

Gilles Dessein
email: gilles.dessein@enit.fr

Compensation of a ball end tool trajectory in complex surface milling

Keywords: milling, cutting forces, deflection, compensation

This work is consecrated to the minimising of machining errors based on a method for the compensation of the trajectory to be machined in hemispherical milling. This compensation is found to be necessary because of the tool deflection due to the cutting forces. In order to remedy to the machining errors, caused by this deflection, a compensation method has been proposed. The latter is inspired from the mirror method, since the compensated position is going to be determined as being the trajectory reflection, deviated onto the mirror. The advantage of this proposed method is that it takes into account the three deflections d x , d y and d z , respectively to the directions X, Y and Z. After that, two-parallel machinings, separated by a groove and achieved absolutely in the same conditions and with the same tool, are carried out, on the same complex part. The first machining is with compensation, but the second is without compensation.

The coordinates of the two obtained surfaces are recorded by a 3D measuring machine. The comparison of the two-surfaces shows the presence of an important correction of the tool trajectory, and reveals a similarity between the part obtained by simulation and the one conceived in CAM.

Introduction

In high-precision industrial machining environments, the machined surface quality and precision constitute important criteria at the time of machining. This precision may be altered by some machining errors. Among these errors, the one induced by the thermal effect can be mentioned. [START_REF] Huang | !Optimization of fixture design with consideration of thermal deformation in face milling[END_REF] have proposed an effective method of fixture design to reduce the machining errors caused by the cutting heat in face milling. [START_REF] Liu | !A study on optimal compensation cutting for an aspheric surface using the Taguchi method[END_REF] afforded us with detailed processes for machining an aspheric surface which used a four-axis desktop ultraprecision and the optimisation of compensation cutting has used the Taguchi method. Concerning the tool wear [START_REF] Matsumura | !Cutting process of glass with inclined ball end mill[END_REF] have suggested a compensation method for the wear or the deterioration of the tool cutting edge. A large portion of the dimensional errors is caused by positioning errors in the axes of motion on a machine-tool. [START_REF] Fines | !Machine tool positioning error compensation using artificial neural networks[END_REF] have fulfilled the application of artificial neural networks to the problem of calculating these error compensation values. In this context, [START_REF] Wan | !A unified framework of error evaluation and adjustment in machining[END_REF] have proposed an effective model that includes machine-tool, fixture and datum errors, within a unified framework, by means of a differential motion methodology. Furthermore, they present a compensation model of fixture errors, datum errors and the deviation of the tool coordinates.

In the framework of ball-end milling, a machining error compensation, caused by the tool deflection, is proposed. Several methods, allowing the compensation of the tool trajectory, following its deflection, are proposed in the literature. Indeed, [START_REF] Budak | !Peripheral milling conditions for improved dimensional accuracy[END_REF] present an approach aiming to optimise the cutting condition, with the purpose to limit, to the maximum, the machining errors, in order to correct them afterwards. [START_REF] Law | !Error compensation in the end milling of pockets: a methodology[END_REF]Geddam (2001, 2003) try to measure and control the cutting forces or the deflection errors in real-time. Their method requires important means to establish a control system, and carry out the necessary modification on the machine.

Another approach consists in compensating the tool trajectory, with the modification of this trajectory, following the prediction or the measuring of the machining errors. The latter can be estimated with a neural networks [START_REF] Cho | !Integrated machining error compensation method using OMM data and modified PNN algorithm[END_REF] or by an analytic/numerical model of the tool deflection [START_REF] Cho | !Integrated error compensation method using OMM system for profile milling operation[END_REF][START_REF] Suh | !Incorporation of tool deflection in tool path computation: simulation and analysis[END_REF].

In this context, a new method, tending to compensate the tool trajectory, has been proposed. This approach is inspired from the mirror method [START_REF] Dépince | !Active integration of tool deflection effects in end milling. Part 2. Compensation of tool deflection[END_REF][START_REF] Seo | !Integration of tool deflection during the generation of tool path[END_REF], since the compensated position is going to be found out as being the reflection of the trajectory deviated onto the mirror.

Indeed, the same principle will be adopted. However, this method is going to be generalised simultaneously in the three directions X, Y and Z, by taking into account the three flexions d x , d y and d z . The latter are determined according to the cutting forces F X , F Y and F Z .

In this work, the different steps to follow are going to be presented one by one, in accordance with their chronological succession, in order to end, finally, at the effective machining.

After the elaboration of the cutting forces model [START_REF] Smaoui | !Simulation of cutting forces for complex surfaces in ball-end milling[END_REF], all the interest is focused on the study of the tool deflection, due to these forces, with the finite elements method. This method consists in studying the tool deflection, in the three directions under the effect of the forces F X , F Y and F Z applied simultaneously. It also consists in presenting the proposed compensation method. At this level, the problem consists in determining the deflected trajectory, on which the mirror method will be applied. So, the simulation of the deflected trajectory is going to be done, according to the average force between two-consecutive nodes. The whole of nodes form the CAM trajectory. The forces cited above depend on the inclination angle g according to the X axis, measured instantaneously by the tangent.

The deflected trajectory is the vector formed by the CAM trajectory, deviated by the average forces and calculated previously. The choice of the average forces is justified by the inability to indicate, with precision, the force-intensities at the level of each CAM point. These force-intensities can be comprised only between a minimum and a maximum.

A compensation algorithm, inspired from the mirror method, is going to be adopted. It will allow to correct the tool deflected trajectory, all along the part.

The coordinates, which constitute the corrected tool trajectory, are sent afterwards to the CAM software (MasterCam © ), to control the machine.

Finally, two-parallel machinings, separated by a groove and executed absolutely in the same conditions and with the same tool, will be carried out on the same complex part, the first with compensation while the second is done without compensation.

The coordinates of the two obtained surfaces will be taken by a 3D measuring machine. The difference between these two-surfaces corresponds to the deflection average value, occurring in that area. It was ineluctably shown that the compensated surface has improved of the deflection value, which makes it considerably and effectively nearest to the surface conceived in CAM.

Calculation and simulation of the cutting tool deflection error

Cutting forces model

The cutting forces model is based upon the discretisation of the tool, into a series of elementary discs, all along the axial depth. The local radius R(z) of the elementary disc, or of each circumference, can be given under the following form:

* + * + 2 2 R z R R z ? / / (1)
With R the tool radius and !z" the height of each disc (Figure 1). The elementary cutting forces can be calculated at any point P, the middle of each discretised section, in the two-directions. The first is an axial direction, following the Z axis in N z increments. Each one bears the index, i(i = 1, 2, #, N), counted from the point O. The second is an angular discretisation in N s increments, counted from the Y axis. Each increment bears the index j (j = 1, 2, #, N s ). The total cutting force for the j position is:

] ̲ 1 1 ( ) ( ) [ ] ( , , ) sin ( , , ) ( ) f z s b c N X R N Y T z i k Z A dF j K dF j T i j k K f i j k dz dF j K d ? ? Ç Ç È Ù È Ù ? È Ù È Ù È Ù È Ù É Ú É Ú ÂÂ (2)
matrix of the transfer from ( , , , ).

S C R T A ̇if if if to ̇C (O, X, Y, Z) such as:
] ̲ sin sin cos cos sin sin cos sin cos cos cos 0 sin First of all, this model is applied in the case of a plane surface, to end at a generalisation, carried out in the case of concave or convex inclined surfaces [START_REF] Smaoui | !Simulation of cutting forces for complex surfaces in ball-end milling[END_REF].

S C T j d d j d j d d j d j j ̇/ / / Ç È Ù ? / / È Ù È Ù / É Ú (3) N f is
The study of this part milling will be afterwards elaborated, while taking into account the nature of the trajectory to be machined, and of the machining direction. It also takes into account the tool geometrical parameters, such as the axial and the radial depth of cut.

Determination of the specific coefficients

In order to determine the specific coefficients, in the three directions: tangential, radial and axial, respectively noted K R , K T and K A , the cutting forces F X , F Y and F Z must be, first of all, measured. These parameters are bound by a model of cutting forces [START_REF] Milfelner | !Simulation of cutting forces in ball-end milling[END_REF]. These coefficients necessarily vary according to the slope of the trajectory to be machined, and also according to the cutting conditions such as: feed per tooth f zb , the cutting speed V c and the depth of cut: axial A d and radial A r . Indeed, the changing of one of these parameters implies new coefficients.

For this reason, the specific coefficient behaviour is going to be defined according to the inclination angle of the surface to be machined, counted from the horizontal. A series of tests will be carried out, beginning with the machining of a surface of 75° of inclination, and diminishing each time of an angle of 15°, until the horizontal state, as indicated in Figure 2. The adopted trajectory is a zigzag one: the tool carries out several go-and-backs, following X and !X, which respectively correspond to an upward cut and a downward cut.

However, the specific coefficients K R , K T and K A cannot be determined without having recourse to matrix of a reverse passage [T] !1 , which allows the passage from the elementary cutting forces dF X , dF Y and dF Z in the Cartesian coordinate system ̇C (O, X, Y, Z) towards the elementary cutting forces in the local coordinate system ( , , , )

S C R T A ̇if if if such as: } ' ] ̲ } ' 1 , , , , R T A X Y Z dF T dF / ? (5
j d d j d j d d j d j j / / / / Ç Ç Ç È Ù È Ù È Ù ? / / È Ù È Ù È Ù È Ù È Ù È Ù / É Ú É Ú É Ú (6)
The specific coefficients K T , K R and K A are determined from the following equation [START_REF] Milfelner | !Simulation of cutting forces in ball-end milling[END_REF]:

] ̲ * + * + 1 1 1 1 1 1 1 , , sin , , f z f z b X N N N Y i j k R Z T N N N A z i j k F T i j k F K F K K f i j k dz s s d / ? ? ? ? ? ? Ç È Ù È Ù È Ù Ç È Ù È Ù É Ú ? È Ù È Ù Ç É Ú É Ú ÂÂÂ ÂÂÂ (7)
A ball-end cutting tool, 16mm in diameter, having two-teeth (N f = 2) and an angle of helix i 0 = 9°, has been used. The part to be machined is made of steel, type XC38 (C35).

The cutting conditions are defined by an axial depth of cut A d = 1 mm, a radial width of cut A r = 3.873 mm, a feed by tooth f zb = 0.1 mm/tooth, and a spindle speed equal to 4,000 rpm, on a high speed milling machine (HSM). The cutting forces for the extreme angles of inclination, measured by a dynamometer platform (Kistler 9275B). Signals are analysed, due to the software Dynoware © . Figure 3 shows cutting forces for g = 0° and g = 75°.

A magnified length, sampled from the cutting forces strips for this dynamometer surface, corresponding to g = 0, and compared to the simulations obtained by MATLAB © software, is presented in Figure 4.

Figure 4 depicts that there is a correlation between the results obtained in simulation, and those obtained by the dynamometric plate.

The purpose of these series of tests, achieved for the different surfaces, and according to different inclinations, is to determine the specific coefficients variation, according to these slopes. This variation is presented in Figure 5. The variation of the specific coefficients must be introduced in the simulation program. For this reason, a curve of tendency is necessary, in order to find out the theoretical expression of the coefficients K R , K T and K A according to the angle of inclination g.

7
06 5 ! 0.0013 4 0.0899 3 ! 2.7111 2 27.9 500

A K E c c c c c ? / - - - (7)
1 05 5 0.0036 4 ! 0.3341 3 14.867 2 ! 325.78 3600

T K E c c c c c ? / / - - - (8)
3 05 5 0.0062 4 ! 0.5215 3 20.971 2 ! 409.27 3400

R K E c c c c c ? / / - - - (9)

Calculation of the cutting tool deflection

In order to end at the compensation step, it is necessary to know how to determine the tool deflected trajectory, which is the deviated CAM trajectory. For this, it is necessary to simulate the tool deflected trajectory, according to the average force between two-successive nodes.

First, the cutting forces model has been applied on a plane surface. Second, it was compared to an experimental result, which has resulted in a perfect correlation. After, this same model has been generalised, in the case of an inclined surface, and another circular surface, but, by adopting each time a change of reference, depending on the angle of inclination [START_REF] Smaoui | !Simulation of cutting forces for complex surfaces in ball-end milling[END_REF].

The principle of this method consists, first of all, in determining the angle of inclination g, for each segment joining these two-nodes (Figure 6), and in calculating the average force between two-successive nodes Once the cutting forces all along the tool trajectory are determined, the pertinent choice goes towards the finite elements method, in order to determine the tool deformation. It!s a matter of applying a well determined force on the extremity of the tool in the X direction. So, it is necessary to find out its deflection in the three directions X, Y and Z. Then, the force directions is changed in the Y way then in the Z way, and deduce, once more, the tridimensional deflection value according to the force intensity. The tool behaviour is presented further to the force exerted, following the X direction (F = 1,000 N) (Figure7).

The total deflection is the sum of all the deflexions, due to the forces F X , F Y and F Z , applied simultaneously, following X, Y and Z axis, deduced by the finite elements method. It is presented as follows: (10) whole of nodes, constituting the CAM trajectory, are deviated according to the average cutting forces, calculated instantaneously, all along the cutting tool trajectory, swept during the machining. The CAM trajectory is instantaneously replaced by the points, which form the deflected (adopted) trajectory coordinates. It is this same trajectory which is going to be corrected by the mirror method before machining.

X Y Z X Y Z X Y Z F F F X Y Z F F F X Y Z F F F X Y Z dx dx dx dx F F F

Compensation method

In order to compensate this deflection, a new process reveals to be necessary. The method proposed is based upon the tool trajectory modification. It is inspired from the mirror method [START_REF] Hascoët | !Compensation des déformations d"outils pour la génération de trajectoires d"usinage[END_REF]. The latter is an iterative method, conceived for a contour machining, where the deflection is concentrated in the normal direction of the tool [START_REF] Dépince | !Active integration of tool deflection effects in end milling. Part 2. Compensation of tool deflection[END_REF]. In this way, the mirror method is limited to a compensation, following only one direction, and so neglects the tangential direction. In this study, the deflection occurs on three directions X, Y and Z. Hence, the need to a new method, able to remedy to this deflection. The method proposed here, considers the two-trajectories: one is desired (CAM), having the coordinates (x desired , y desired , z desired ), and the other is deflected. A first stage consists in compensating the trajectory of the values dx, dy and dz, in the symmetrically opposite direction, according to the deflected trajectory. Figure 8 presents a simplified diagram of the compensated method for only one direction dx among the three mentioned directions dx, dy and dz. For this new position, a variation of the axial depth of cut A d , of the radial width of cut A r and of the feed by tooth f zb , respectively to the three directions

X F Z X Y Deflection value x f z h
X, Y and Z, will appear. For this reason, the cutting forces and the new corresponding deflections, have to be, once more, calculated. These steps are going to be repeated i-times, in an iterative way, until the deflections will be inferior to a tolerance already fixed at the beginning of this procedure. For this example, the tolerance interval is fixed to 0.005 mm. In this way, it becomes possible to reach the optimal trajectory, which will be imposed to the cutting tool (x imposed , y imposed , z imposed ). This vector will be sent to the MasterCam © software, in order to achieve the part machining, according to the new coordinates.

Application of the compensation method

We want to carry out the milling of a complex workpiece, having six-shapes of surfaces, with a length L = 220 mm and a width l = 78 mm, as indicated in Figure 9.

The same tool used here, is a ball-end tool of a ray equal to 8 mm, having two-teeth (N f = 2) and an angle of helix i 0 = 9°. This tool fits in the workpiece with a feed of a value f zb = 0.1 mm/tooth and an axial depth of cut A d = 1 mm. The two-finish tests: with and without compensation, are going to be realised on the same part. That is why, the two-test surfaces are going to be separated by a groove (12 mm width and 1 mm depth), so the two-machinings will be carried out under the same conditions (Figure 10). The cutting forces are calculated according to the nature of the trajectory part to be machined (flat, inclined!), and following the machining direction (following X or "X).

Figure 11 shows the cutting forces evolution, following the tool-end trajectory. This trajectory is obtained from the coordinates of the points transferred from the software MasterCam © . It is of a zigzag type: the tool carries out several go-and-backs, following X and "X. This respectively corresponds to an upward cut and a downward cut. The cutting forces are drawn in black for the X machining, and in red for the "X machining.

After simulating the cutting forces, the tool deflection is determined according to the average forces, all along the trajectory swept on the surface. The coordinates of the nodes forming the deflected (adopted) trajectory, are gathered in a vector. It is on this trajectory that the correction of the tridimensional trajectory is going to executed with the mirror method. This leads to find the compensated trajectory, which is going to be sent towards the machine (Figure 12), by means of a numerical control file (NC) generated by MasterCam © . The coordinates of the points forming this path, will be sent to MasterCam © . So, the new compensated tool trajectory could be simulated. The two-NC files, containing the compensated and non-compensated trajectories, are sent to the machine, in order to tackle the machining in compensation.

Results

The two-compensated and non-compensated surfaces are measured on a 3D measuring machine (DEA GAMMA RECORD).

In order to confirm the validity and the precision of the compensation program, the compensated trajectory is going to palpate in order to determine the coordinates of the points constituting this trajectory (Figure 13). The coordinates noted down on the part, in the two-machined zones, (with and without compensation), made is possible to deduce a variable and important gain in the different tested zones represented in Figure 14. For the same coordinate in X, The difference in amplitude (Z) for a sampling of the two-compensated and non-compensated surfaces constitutes the amplitude of the correction thus obtained. The latter is listed in Table 1. This geometrical difference, between the two-compensated and non-compensated surfaces, should normally be the nearest possible to the deflection simulated theoretically, in the same area. This area is theoretically located between two very successive nodes.

A second validity measurement has been achieved on the convex zone (zone 4), supposed to be 16 mm in diameter. Four-samplings of coordinates have been carried out following X: (141.8543; 155.1659; 161.9415; 170.4063), which corresponds respectively to a theoretical deflection equal to (!0.0002; 0.1126; 0.1107; 0.0011) in the X direction and (0.0078; 0.0142; 0.0142; 0.0001) in the Z direction.

An average deflection of the order of 0.056 has consequently been remarked. The touch probe has noted down an average diameter of 60.345 mm in the non-compensated area. However, for the compensated zone, the touch probe has recorded a diameter of 60.236 mm, as indicated in the two following diagrams, given by Metrolog © software. It is to be noted that, on one hand, the deflection in the X direction exerts a certain influence on the obtained curvature, since it is more important. On the other hand, and thanks to the recorded values, the radius value of the compensated zone is estimated to 30.173. Nevertheless, the radius value of the other non-compensated zone is of the order of 30.118. So, it becomes possible to deduce from this, that the difference between these two-radii is equal to 0.054. This corresponds to the average value of the deflection occurring in this zone (Figure 15). These results suitably concord with the values expected in simulation. These results unequivocally confirm the results provided in this work and the precision of the methods chosen.

Conclusions

In this paper, a compensation method of the trajectory error is presented in order to be machined, which is based on the cutting tool deflection. The elementary cutting force has been calculated, for each discretised section, by means of the specific cutting force coefficients, in the local spherical reference ( , , , ).

S C R T A ̇if if if

These coefficients have been extracted from the experimental results. Their theoretical expressions, according to the trajectory inclination, are injected in the compensation program. The tool, used here, is 16mm in diameter, and the piece to be machined is made of steel type XC38 (C35).

The tool deflection, caused by these forces, has been studied once their simulation has been realised. This study is a preliminary step which makes it possible to carry out a correction of the trajectory, according to the machining tolerance.

Thus, the results obtained by the finite elements method are going to be used, in order to calculate the tool deflection in the three-directions according to the applied forces F X , F Y and F Z .

After that, a compensation method has been proposed for a change of the tool trajectory. The new coordinates of the points, which constitute the new trajectory, are exported towards the CAM software. This method has been applied in the case of a complex part machining. The simulation results have shown the advantage of the approach proposed since it has allowed to correct the machining errors and ended at an optimal trajectory respecting the machining tolerances.

Finally, two-parallel machinings, separated by a groove, have been achieved, on the same complex part and carried out absolutely in the same conditions, and with the same tool. The first machining corresponds to trajectory with compensation and the second to trajectory without compensation.

The coordinates of the two obtained surfaces are read by a 3D measuring machine. The comparison, between these two-surfaces, confirms the similarity in the results obtained by the two-methods, one with compensation and the other conceived in CAM.

Figure 1

 1 Figure 1 Cutting forces applied to the tool (see online version for colours)

  the number of the tool teeth. K R , K T and K A are respectively the specific coefficients: radial, tangential and axial, empirically defined. d is the angle of width of cut

Figure 2

 2 Figure 2 Different slopes to test (see online version for colours)

Figure 3 Figure 4 Figure 5

 345 Figure 3 Variation of the cutting forces for the different angles of inclination in ̇C (O, X, Y, Z) (see online version for colours)

Figure 6

 6 Figure 6 Average cutting force determined all along the trajectory (see online version for colours)

Figure

  Figure 7 Deflection ( ) X F dx following the X axis for F = 1,000 N (see online version

Figure 8

 8 Figure 8 Compensation method (see online version for colours)

Figure 9

 9 Figure9Workpiece to be machined in 3D (see online version for colours)

Figure 10

 10 Figure 10 Separation of the two-surfaces (see online version for colours)

Figure 11 Figure 12

 1112 Figure 11 Cutting forces simulation in longitudinal milling, (a) FX (b) FY (c) F z (see online version for colours)

Figure 13

 13 Figure 13 Measurement of the geometry obtained by the two-tests (see online version for colours)

Figure 14

 14 Figure 14 Geometrical characteristics of the different tested zones

Figure 15

 15 Figure 15 Tool deflection recorded in the direction of X (see online version for colours)

Table 1

 1 Geometrical correction obtained for the case of a downward cut

	Zone	Value of the correction	Estimated theoretical deflection
	1	!0.002	!0.002
	2	0.014	0.014 et 0.015
	3	0.011	0.007 et 0.014
	4	0.009	Between 0.006 and 0.01
	5	!0.019	Between !0.017 and 0.017
	6	0.003	0.03


Djimedo Kondo
email: djimedo.kondo@univ-lille1.fr

Hélène Welemane
email: helene.welemane@enit.fr

Fabrice Cormery
email: fabrice.cormery@polytech-lille.fr

Basic concepts and models in continuum damage mechanics

Keywords: MOTS-CLÉS : endommagement, poroélasticité, micromécanique, ssures damage, poroelasticity, micromechanics, cracks. 1

Dans cet article, on présente dabord quelques éléments de base de la modélisation macroscopique de lendommagement. On rappelle notamment le cadre thermodynamique des processus irréversibles et son application à la modélisation de lendommagement isotrope. Létude de lanisotropie induite par lendommagement est ensuite traitée en considérant une variable tensorielle dordre 2. Enn, larticle sachève sur une contribution originale dans laquelle on sintéresse à la modélisation des effets unilatéraux liés à la refermeture des microssures. ABSTRACT. In this paper, we present some basic elements of macroscopic modelling of damage. We then recall the general approach of continuum damage based on the thermodynamics of irreversible processes and its application to isotropic damage modelling. The study of damageinduced anisotropy is treated by considering a second order tensorial damage variable. Finally, we present an original macroscopic approach through which is addressed the question of unilateral effects due to the microcracks closure.

Some basic concepts in Continuum Damage Mechanics (CDM)

First of all, let us recall that the main objective of standard CDM is to propose a continuum-mechanics based framework allowing to characterize, represent and model at the macroscopic scale the effects of distributed defects and their growth on the material behavior. This requires the consideration of some basic concepts, some of which are recalled in the present section. To introduce these basic concepts of CDM 1 , it proves useful to consider rst an uniaxial tension experiment.

The classical unixial damage theory

Within the classical approach (see for instance [START_REF] Lemaitre | A course on Damage Mechanics[END_REF]), a very simple measure of the damage amplitude in a given plane is obtained by measuring the area of the intersection of all defects with that plane. For example, based on gure 1, it is readily observed that the effective area of the sample subjected to uniaxial tension is S -S D . S D represents the defects trace in the considered plane. The following positive scalar ω is then commonly considered as a damage variable in the above 1D experiment :

ω = S D S [1]

Figure 1. Cross section of a damaged material

For the undamaged material, S D = 0 and then ω = 0. The damage being related to the growth of defects, ω may grow from 0 to a critical value often taken in litterature equal to 1 which corresponds to an entirely damaged material (effective area S -S D reduced to 0). Instead of the standard uniaxial stress σ = F S , it is convenient to introduce for the damaged material the effective stress :

σ = F S -S D = F S(1 -S D /S) = σ 1 -ω [2]
in which use has been made of [1].

Associated to a strain equivalence principle, the 1D effective stress σ is related to the elastic strain of the material by the uniaxial Hookes law :

σ = Eε [3]
where E is the elastic modulus of the undamaged material. It follows that the constitutive law for the standard stress σ takes the form :

σ = (1 -ω)Eε [4]
For the uniaxial model formulation, equation [4] must be completed by the damage evolution law which can be considered in the form of a dependence between the damage variable ω and the applied load : ω = g(ε)

[5]

A priori, the function g can be identi ed from uniaxial tension test. It must be noted that the relation between ω and ε is valid only in the monotonous loading regime. In an unloading and reloading phase, the damage variable kept its maximum value reached before. A classical way to describe in a uni ed manner these different loading regimes consists in introducing a variable κ which characterizes the maximum level of strain reached in the material before the current time t : κ(t) = max ε(τ ) for τ ≤ t. The damage evolution relation [5] can then be recast in the form : ω = g(κ) [START_REF] Boehler | Application of tensors functions in solids mechanics[END_REF] which remains valid for any kind of loading regime. The complete elastic response of the damaged material is schematized in gure 2. Instead of considering the function g, it is usual to introduce a limit state function :

f (ε, κ) = ε -κ [7]
f ≤ 0; κ ≥ 0; κf = 0 [8]
The condition f ≤ 0 indicates that ε can never be greater than κ, while the second condition means that κ cannot decrease. Besides, the second condition implies that κ can increase only if the current value of the strain is equal to κ.

In summary, the basic elements of the above uniaxial damage theory are as follows : the stress-strain law : σ = E s ε = (1ω)Eε ; this relation appears as a classical Hookes law with a secant modulus E s = (1ω)E associated to the damaged material. A simple determination of the damage variable is then : D = 1 -Es E . It requires however a careful and accurate measure the elastic strain, a damage evolution law wich can be put in the form ω = g(κ) or a a limit state function f ; a rst approach can consist to choose f (ε, κ) = εκ.

A simple isotropic damage theory

A straightforward 3D extension of the previous uniaxial theory is given by the well known Lemaitre-Chaboche [START_REF] Lemaitre | Aspect phénoménologique de la rupture par endommagement[END_REF] 3D elastic damage model. In this model, based on the various concepts presented in section 1, it is postulated that the stiffness D of the damaged material reads :

D = (1 -ω)D 0 [9]
in which D 0 denotes the elastic stiffness of the undamaged material and ω is de ned by [1].

It is readily seen that the generalization of the uniaxial stress-strain constitutive law [4] takes the form :

σ = (1 -ω)D 0 : ε [10]
and the corresponding effective stress tensor σ is given by :

σ = (1 -ω)σ [11]
In the present three dimensional formulation of the isotropic damage model, only the Young modulus is affected, the Poisson ratio remains constant during the damage process. This is clearly a shortcoming of this simple model and will be corrected in section 2.

Let us come now to the 3D loading function f and to the damage evolution. f de nes in the 3D strain space the domain of elasticity whose boundary corresponds to the strain states at which the damage will evolve. An immediate generalization of [7] reads :

f (ε, κ) = ε eq (ε) -κ [12]
in which the equivalent strain ε eq is a norm of ε that needs to be choosen. A rst simple choice ofor this can be :

ε eq = ε : ε [13]
Another choice can be the elastic energy (function of the strain) : ε eq = ε : D 0 : ε.

An isotropic damage modelling of concrete materials

The two norms of ε introduced before for the de nitions of ε eq leads to a symmetric elastic domain in tension and compression. However, several materials (rocks, concrete, ceramics) often show a dissymetric damage surface, the yield value in compression being several times the value in tension. In order to overcome these limitations, Mazars [START_REF] Mazars | Application de la mécanique de lendommagement au comportement non linéaire et à la rupture du béton de structure[END_REF] introduced two damage parameters, ω t associated to a tension mechanism and ω c devoted to the damage under compression. These two parameters, ω t and ω c , are evaluated from two evolution functions, g t and g c , which are assumed to depend both on a unique de nition of the equivalent strain :

ε eq = 3 I=1 < ε I > 2 [14]
where < . > is the Macauley bracket and ε I are the principal strains. < ε I > denotes then the positive part of the principal value ε I of the strain tensor.

The objective of the constitutive model initially proposed by [START_REF] Mazars | Application de la mécanique de lendommagement au comportement non linéaire et à la rupture du béton de structure[END_REF] is to capture the non linear response of geomaterials subjected to loading paths which involves extension. For a general loading path, different from uniaxial tension or compression, the value of damage which enters in the constitutive law is proposed as the following combination :

ω = α t ω t + α c ω c [15]
where α t and α c are taken in the form [START_REF] Pijaudier-Cabot | Damage models for concrete[END_REF] :

α t = 3 I=1 < ε t I >< ε I > ε 2 eq β ; α c = 3 I=1 < ε c I >< ε I > ε 2 eq β [16]
Obviously, in uniaxial tension α t = 1; α c = 0 and then ω = ω t while under uniaxial compression α t = 0; α c = 1 and ω = ω c . An integrated form of the damage evolutuion is (see [START_REF] Pijaudier-Cabot | Damage models for concrete[END_REF]) :

ω t = g t (κ) = 0 if κ ≤ κ 0 1 -(1-At)κ0 κ -A t exp[-B t (κ -κ 0 )] if κ ≥ κ 0 [17]
and

ω c = g c (κ) = 0 if κ ≤ κ 0 1 -(1-Ac)κ0 κ -A c exp[-B c (κ -κ 0 )] if κ ≥ κ 0 [18]
where κ 0 = ft E is the equivalent strain at the begining of the damage process. In this de nition, f t is the tensile strength of the material and E is the modulus of the undamaged material. Let us precise that the total number of the model parameters is 8 : the two elastic coef cients ; κ 0 , A t , B t , A c , B c are materials parameters which can be derived from compressive and tensile tests ; the determination of the constant β requires a shear test (for simplicity β is often taken equal to 1).

An example of simulation of a concrete behavior by means of this model is reported in gure 3.

Thermodynamic framework of CDM : application to isotropic damage

The damage models presented in the preceeding section can be considered in a suitable framework given by the Thermodynamics of Irreversible Processes (TIP). In particular, this TIP framework allows a relatively easy modelling of anisotopic damage as in section 3.

Figure 3. Simulation of the uniaxial (tension and compression) response of a concrete.

From [START_REF] Pijaudier-Cabot | Damage models for concrete[END_REF]

The internal damage variable approach

Generally speaking, damage is de ned as the modi cation of physical properties of materials in relation with the irreversible growth of microdefects. The macroscopic approach requires rst the choice of some appropriate variables physically motivated and mathematically relevant for the description of the effects of microdefects at macroscale. Clearly enough, is desirable that the damage variable has some physical meaning and a mathematical coherence.

A representation of the damage can be obtained from a geometrical characterization of the cracks-like defects. When the damage process is due to the growth of spherical voids as in metallic materials, ω corresponds to the voids volume fraction, i.e. the porosity. In contrast, the non linear behavior of geomaterials may be related to microcracks growth. For a randomly oriented system of similar penny shaped microcracks for which anisotropic aspects can be disregarded, it was shown that the cracks density parameter ω (see [START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF]) can be considered from a micromechanics point of view as a relevant damage variable (see also [START_REF] Dormieux | Micromechanics of damage propagation in !uid-saturated cracked media[END_REF] this issue) :

ω = N a 3 [19]
where N denotes the crack density (number of cracks per unit volume) and a the radius of the penny-shaped cracks (see gure 4 for an example of a representive elementary volume and gure 5 for the geometrical description of a penny-shaped crack).

For non circular microcracks, a generalization of [19] for the representation of the iso-tropic damage has be given by [START_REF] Kachanov | Continuum model of medium with cracks[END_REF] which includes the microcracks area and a shape parameter :

ω = N i ηS 3/2 i [20]
where S i is the decohesion area related to an i-th cracks family and η is a dimensionless cracks-shape parameter. In the case of anisotropic damage, the cracks orientation n may be considered. A very simple representation of the anisotropy of the damage can be obtained by considering the following second order tensor D based on the orientation tensor n ⊗ n (see again [START_REF] Kachanov | Continuum model of medium with cracks[END_REF]) :

D = i ω i n i ⊗ n i [21]
in which ω i and n i are respectively the cracks density parameter and the orientation of a given cracks family i.

Such representation of damage allows to describe only materials symmetry up to orthotropy. The choice of the most appropriate tensorial representation for the damage variable corresponding to microcracking phenomena is still a much-debated question. The reader interrested by this topic can refer to the book [START_REF] Krajcinovic | Damage mechanics[END_REF].

Thermodynamic potential and state laws

The elastic damage behavior is characterized by the existence of a thermodynamic potential as function of the set of state variables. Strain formulation is considered (elastic strain ε as observable variable) with isotropic damage (represented by a single damage variable ω). Assuming linear elasticity at constant damage, the free energy reads :

W (ε, ω) = 1 2 ε : D(ω) : ε [22]
where D(ω) represents the isotropic stiffness tensor for a given damaged level, such that ω = 0 corresponds to the virgin material. Due to the isotropy, the general form of D(ω) is given later by [31]. The potential [22] must be continuously differentiable :

with respect to the strain tensor in order to ensure the existence of the state laws which give : the stress tensor σ(ε, ω)

σ(ε, ω) = ∂W ∂ε (ε, ω), [23]
with respect to the thermodynamic force F ω (ε, ω) associated to the damage variable :

F ω (ε, ω) = - ∂W ∂ω (ε, ω) = - 1 2 ε : D ′ (ω) : ε [24]
where D ′ = ∂D ∂ω . By similarity with the energy release rate in Linear Elastic Fracture Mechanics (see paper by C. Dascalu and by M. Jirasek, this issue), F ω can be interpreted as a damage energy release rate.

Damage surface and damage evolution law

Before introducing the damage yield surface, let us analyze the dissipation D of the material in the presence of damage. One has :

D = σ : ǫ -Ẇ = - ∂W ∂ω ω [25]
Equation [25] suggests that F ω = -∂W ∂ω can be considered as the driving force of the damage process. Following [START_REF] Marigo | Modeling of brittle and fatigue damage for elastic material by growth of microvoids[END_REF] and others, a quite classical form of the damage yield function f is the following one :

f (F ω , ω) = F ω -κ(ω) [26]
κ is a scalar strictly positive function of ω which represents a resistance of the material to the damage propagation. Due to the dependence of κ on ω, the later plays the role of a hardening variable.

As an example, the following very simple form :

κ(ω) = κ 0 (1 + ηω) with κ 0 > 0 and η > 0.
[27] has been considered in several models of litterature.

For the damage evolution, it is usual to assume the normality rule ; the damage rate is then normal to the yield surface and reads then :

     ω = Λ ∂f ∂F ω (F ω , ω) = Λ Λ = 0 if f (F ω , ω) ≤ 0, ḟ (F ω , ω) < 0 Λ > 0 if f (F ω , ω) = 0, ḟ (F ω , ω) = 0 [28]
in which the damage multiplier Λ is determined by the consistency condition. Reporting this evolution law in the rate form of the stress tensor, one gets the incremental formulation of the elastic damage law :

σ = L : ε [29]
with the tangent operator L given by the following expression :

L =        D(ω) if f (F ω , ω) ≤ 0, ḟ (F ω , ω) < 0 D(ω) - ∂F ω ∂ε ⊗ ∂F ω ∂ε h if f (F ω , ω) = 0, ḟ (F ω , ω) = 0 [30]
in [30], tensor D(ω) is given by [49] and h = κ 0 η -∂F ω ∂d .

Different versions of the isotropic damage model

In summary, the complete formulation of the isotropic model requires the choice of D(ω) or equivalently the one of the thermodynamic potential. As the material remains isotropic, D(ω) takes the following general form :

D(ω) = 3K(ω)J + 2µ(ω)K.
[31]

where (K, µ) represent respectively the compressibility and shear moduli for the damaged material. J and K are the two isotropic fourth order projectors de ned by J = 1 3 1⊗1 and K = I-J with 1 the unit second order tensor. I denotes the symmetric fourth order unit tensor.

Differents choices can be made for D(ω). The model presented in subsection 1.2 corresponds to the simplest one and is de ned by [9] which corresponds to the fact that damage affects the compressibility and shear modulus in the same maner. Consequently, Poisson ratio remains constant as already indicated.

A very ef cient way to build macroscopic expression of D(ω) is provided by micromechanics of microcracked materials. In practice, different schemes (dilute scheme, Mori-Tanaka approach etc.)2 can be considered. As an example, [START_REF] Ponte-Castaneda | The effect of spatial distribution of effective behavior of composite materials and cracked media[END_REF] established the following bound for randomly oriented pennyshaped microcracks with a spherical spatial distribution :

K(ω) = K 0 1 - 48(1 -ν 2 0 )ω 27(1 -2ν 0 ) + 16ω(1 + ν 0 ) 2 [32] µ(ω) = µ 0 1 - 480(1 -ν 0 )(5 -ν 0 )ω 675(2 -ν 0 ) + 64ω(4 -5ν 0 )(5 -ν 0 ) [33]
where ν 0 is the Poisson ratio of the undamaged grains.

Macroscopic modelling of anisotropic damage

Representation theorems and expression of the free energy

Although the choice of the most appropriate tensorial representation for the damage variable is still debated [START_REF] Lubarda | Damage tensors and the crack density distribution[END_REF]), a simple continuous representation of the damage can be adequately obtained by using a second-order symmetric tensor D as introduced in subsection 2.1.

Since only induced anisotropy is considered, it is assumed that in the undamaged state the material is isotropic elastic linear. As for the isotropic damage modelling, the next step is to choose an expression of the free energy related to the anisotropic elastic damage behavior. A classical result of the representation theorems (invariants theory) is that any scalar valued anisotropic function of vectors and second order tensors, can be expressed as an isotropic function of the original arguments and of the additional tensorial arguments inducing the anisotropy. Thus, the elastic energy W of the damaged body can be represented by an isotropic function of ε and the damage tensor D. The most general representation of a polynomial isotropic scalar function of the symmetric tensors (ε, D) is a linear combination of the following invariants (see for e.g. [START_REF] Boehler | Application of tensors functions in solids mechanics[END_REF]) :

tr(ε), tr(ε 2 ), tr(ε 3 ), tr(D), tr(D 2 ), tr(D 3 ) tr(ε.D), tr(ε 2 .D), tr(ε.D 2 ), tr(ε 2 .D 2 ) [34]
Since in the undamaged state (D = 0) or at xed D the material is linear elastic, the energy W must be quadratic in ε. Moreover, if the amplitude of damage is moderate, it can be assumed that W is af ne with D. Thus, the representation of W takes the form :

W = 1 2 [λ(tr(ε)) 2 + 2µtr(ε 2 )]+ 1 2 [ηtr(D)(tr(ε)) 2 + γtr(D)tr(ε 2 )] + αtr(ε)tr(ε.D) + χtr(ε 2 .D) [35]
α, β, η, and γ are constants which characterize the damage effects on the material behavior. It must be noted that this form of potential has been used by [START_REF] Murakami | Constitutive and damage evolution equations of elastic-brittle materials based on irreversible thermodynamics[END_REF] for modelling anisotropic damage of concrete ; similar expression without terms in tr(D) has been also considered by [START_REF] Dragon | A model for anisotropic damage by mesocrack growth ; unilateral effect[END_REF].

Formulation of the anisotropic damage model

The macroscopic stress tensor σ, obtained by partial derivation of W reads :

σ = ∂W ∂ε = D(D) : ε [36]
with

D(D) = D 0 + α(1 ⊗ D + D ⊗ 1) + χ(D ⊗ 1 + 1 ⊗ D) + η tr(D)(1 ⊗ 1) + γ tr(D)(1 ⊗ 1) [37]
where D 0 = 3k 0 J + 2µ 0 K is still the elastic stiffness of the undamaged material.

The second state law allows to introduce the thermodynamic force associated to the second order tensor D :

F D = - ∂W ∂D [38]
The derivation of the damage evolution can be done through the same procedure as for the isotropic damage (subsection 2.3). By analogy with [26] a damage surface, expressed in terms of D and D, can be choosen in the form :

f (F D , D) = F D -0 + a 1 tr(D)). [39]
where a 0 and a 1 are material constants : a 0 de nes the initial damage threshold while a 1 describe the manner in which the surface evolves with damage.

The evolution of D is still assumed to follow the normality rule :

Ḋ = 0, if f < 0 or f = 0 and ḟ < 0 Λ ∂f ∂F D , if f = 0 and ḟ = 0 [40]
One obtains the damage multiplier Λ (positive scalar) :

Λ = tr(F D • Ḟ D ) a 1 tr(F D ) [41]

Isotropic elastic damage behavior with accounts of unilateral effects

In this section, the aim is to extend the isotropic model described in section 2 by incorporating unilateral effects due to microcracks closure. It is well known that the mechanical response of a microcracked medium strongly depends on the opening and closure status of the existing defects in the material. In particular, it is commonly observed that the application of a tensile loading followed by a compression in the same direction leads in the compression regime to a total or partial recovery of the Young modulus. From the modelling point of view, the consideration of the unilateral effects in Continuum Damage Mechanics (CDM) still constitutes a challenging task. The main dif culty lies in the necessity to predict both continuous response of the material and partial or total recovery of the elastic constants during the microcracks closure process. During this microcracks closure, the elastic stiffness must also remain symmetric. A critical review of all these dif culties can be found in [START_REF] Chaboche | Damage induced anisotropy : on the dif culties associated with the active/passive unilateral condition[END_REF]. Several other aspects allowing to point out some spurious energy dissipation are discussed in [START_REF] Carol | Spurious energy dissipation/generation in stiffness recovery models for elastic degradation and damage[END_REF]. Moreover, [START_REF] Cormery | A critical review of some damage models with unilateral effect[END_REF] discussed for existing models the existence of the thermodynamics potential when unilateral effects are accounted.

We propose here a mathematically rigorous and physically-motivated model based on micromechanical considerations3 . Again, we assume here that the microcracks are randomly oriented, so that the damage state is isotropic and can be represented by a positive scalar variable ω which may correspond to the microcracks density parameter (see [START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF]). In order to account for the unilateral effects, a separated description is proposed when the microcracks are opened or closed by D o (ω) and D cl (ω) respectively. The simplest partition of the strain space into two half-spaces is provided by a hyperplan Γ which depends on ε. It is assumed that the opening and closure state does not depend on the microcracks density parameter but only on the strain state. Therefore, the microcracks are opened if Γ(ε) > 0 and closed if Γ(ε) ≤ 0. The isotropic stiffness tensor of the damaged material takes then the following form :

D(ω) = D o (ω) = 3K o (ω)J + 2µ o (ω)K if Γ(ε) > 0 D cl (ω) = 3K cl (ω)J + 2µ cl (ω)K if Γ(ε) ≤ 0 [42] (K o , µ o
) and (K cl , µ cl ) moduli represent respectively the compressibility and shear moduli for opened and closed microcracks.

Continuity conditions and microcracks closure criterion

Since the thermodynamic potential has to be continuously differentiable, (K o , µ o ) and (K cl , µ cl ) must ful ll some conditions. Indeed, [START_REF] Curnier | Conewise linear elastic materials[END_REF] (see also [START_REF] Welemane | Une modélisation des matériaux micro ssurés. Application aux roches et aux bétons[END_REF]) have demonstrated that the elastic energy function W is C 1 continuous if and only if :

[D(ω)] =D o (ω) -D cl (ω) = s(ω) ∂Γ ∂ε (ε) ⊗ ∂Γ ∂ε (ε), ∀ε|Γ(ε) = 0 [43]
where s is a continuous scalar valued function depending on ω. We emphasize that D o is the stiffness tensor of the damaged materials when cracks are open, while D cl corresponds to the stiffness when all cracks are closed.

Using [42], it is readily seen that the jump in the stiffness can be expressed as :

[D(ω)] = 3[K o (ω) -K cl (ω)]J + 2[µ o (ω) -µ cl (ω)]K [44]
which must be singular and in fact of rank one to satisfy the condition [43] [START_REF] Curnier | Conewise linear elastic materials[END_REF]. This condition is ful lled if all determinants of second order obtained from a Voigt representation of [D(ω)] are cancelled, that is if :

∀ω, µ o (ω) = µ cl (ω) [45]
The continuity of W is then obtained if the shear moduli does not depend on the microcracks state : µ(ω) = µ o (ω) = µ cl (ω).

It can be noted that at this step no mathematical conditions are imposed on the compressibility moduli.

Moduli recovery conditions

Inspired by some micromechanical considerations on the progressive closure effects of microcracks, we assume as well that the compressibility modulus corresponding to all closed microcracks takes the initial value K 0 of the sound material. Therefore, the modulus restitution condition reads :

K cl (ω) = K 0 [46]
and we will denote

K o (ω) = K(ω).
At this level of development, it remains to precise the expression of the function Γ which de nes the microcracks closure criterion. In view of [44] and [45], the tensor [D(ω)] is written in the form :

[D(ω)] = 3[K o (ω) -K cl (ω)]J = [K o (ω) -K cl (ω)]1 ⊗ 1 [47]
Comparing [43] and [47], it follows that the closure/opening criterion reads :

Γ(ε) = tr(ε) = 0 [48]
It must be emphasized that the simplicity of the microcrack opening-closure transition criterion is due to the isotropy of the damage. In the case of a anisotropic damage, the methodology followed in this work would lead to a more complex criterion which depends in particular on the microcracks orientation.

Final expression of the thermodynamic potential

Finally, the thermodynamic potential takes the following form :

W (ε, ω) = 1 2 ε : D(ω) : ε, D(ω) = 3K(ω)J + 2µ(ω)K if tr(ε) > 0 3K 0 J + 2µ(ω)K if tr(ε) ≤ 0 [49]
In order to provide a physical basis to the damage model, we consider as before for K(ω) and µ(ω) micromechanical results established by [START_REF] Ponte-Castaneda | The effect of spatial distribution of effective behavior of composite materials and cracked media[END_REF] for microcracked media [32].

For the derivation of the incremental formulation we adopt the same loading function and follow the same procedure as in section 2. The main difference is that we have here two different expression of the thermodynamics potential according to the fact that damage is activated or not. It follows the rate form of the stress tensor for the elastic unilateral damage law :

σ = L : ε [50]
with, e being the deviatoric part of ε, the tangent operator L given by : L = 3k 1 J + 2k 2 K -2k 3 (1 ⊗ e + e ⊗ 1) -2k 4 e ⊗ e [51]

k 2 = µ(ω), k 4 = 2(µ ′ (ω)) 2 h , k 1 = K(ω) -H 1 and k 3 = K(ω) -H 3 if tr(ε) > 0 k 1 = K 0 and k 3 = 0 if tr(ε) ≤ 0 [52]
with H 1 = (K ′ (ω)trε) 2 h and H 3 = µ ′ (ω)K ′ (ω) h . The quantity h is expressed as : h = H 0 η+K ′′ (ω)tr(ε) 1+µ ′′ (ω)e where K ′′ (ω) and µ ′′ (ω) are the second derivative of K(ω) and µ(ω) respectively whereas K ′ (ω) and µ ′ (ω) are the rst derivative of K(ω) and µ(ω) respectively. Different versions of the isotropic damage models can be obtained by considering, as previously, expressions of K(ω) and µ(ω) provided by micromechanics (see again [START_REF] Dormieux | Micromechanics of damage propagation in !uid-saturated cracked media[END_REF]). Obviously, recalling that we refer here to a purely macroscopic approach, any suitable function can be considered.

Conclusions

This chapter is devoted to a general presentation of the basic principles of Continuum Damage Mechanics and the associated thermodynamics framework. It may also provide a general view of the classical procedure used for isotropic damage. In order to take into account the damage-induced anisotropy, representation theorems appear as a useful tool ; however, it must be pointed out that the approach presented here is limited to moderate damage and there is need to go beyond this hypothesis. A suitable modelling of the unilateral isotropic damage is also proposed by imposing some continuity conditions as well as the condition of elastic moduli recovery at damage deactivation. This approach for the unilateral damage modelling has been succesfully developed by [START_REF] Welemane | Une modélisation des matériaux micro ssurés. Application aux roches et aux bétons[END_REF] (see also [START_REF] Welemane | An alternative 3D model for damage induced anisotropy and unilateral effect in microcracked materials[END_REF]) in the context of anisotropic damage. Obviously, there other important topics which are not adressed in this introductory paper : Poroelastic damage [START_REF] Dormieux | Microporomechanics[END_REF], coupling between plasticity and damage, effect of initial anisotropy, friction on closed cracks faces, etc.).

Figure 2 .

 2 Figure 2. Elastic response of a damaged material during loading-unloading

ÒFigure 4 .Figure 5 .

 45 Figure 4. R.e.v. of the cracked material

For the methodology of these methods, the reader can refer to paper[START_REF] Dormieux | Micromechanics of damage propagation in !uid-saturated cracked media[END_REF], this issue.

As we will see, some micromechanical results can be incorporated in order to enhance the purely macroscopic modelling.