Continuum damage modelling in geomechanics

Continuum damage mechanics is a framework for describing the variations of the elastic properties of a material due to microstructural degradations. This chapter presents the application of this theory to the modelling of concrete. Several constitutive relations are devised, including incremental, explicit, and non local damage models. A general framework for damage induced anisotropy is also presented. Coupled damage and plasticity modelling is discussed. Finally, the issue of the experimental determination of the internal length in non local models is tackled.


Introduction
This chapter is concerned with the presentation of several continuum damage models and their implementation in non linear finite element analyses. Continuum damage means that the mechanical effects of progressive microcracking, void nucleation and growth are represented by a set of state variables which act on the elastic and/or plastic behaviour of the material at the macroscopic level. First, we will deal with the scalar damage model and present different versions of this type of constitutive relations. After having recalled how the local integration of the constitutive relation can be performed in the same spirit as for plasticity-based models, we will deal with more specific applications: cases where the evolution of damage is specified in an integrated form, where damage induced inelastic strains are introduced, wheye crack closure effects are modelled, and finally where damage is coupled to plasticity. The following section deals with damage induced anisotropy. A general framework for capturing directionality of damage is recalled. Coupling with plasticity and crack closure effects are also discussed. Non local damage is developed in the fourth section. We present integral and gradient models. In particular, a general scheme for the implementation of the gradient model, implemented into the finite element code "Code_ Aster", is described. The fifth section deals with several computational issues related to damage such as the implementation of the model within secant or tangent algorithms. The chapter concludes on the problem of the calibration of the model parameters in damage models including the internal length. Inverse analysis of size effect tests on geometrically similar specimens is used for the experimental determination of these parameters.

Isotropic Damage Model
When damage is assumed to be isotropic, it is considered that it produces a degradation of the elastic stiffness of the material through a variation of the Young's modulus: (1) where v 0 , E 0 and ~J are the Poisson's ratio and Young's modulus of the undamaged isotropic material and the Kronecker symbol respectively. The elastic (i.e. free) energy per unit mass of material is: (2) where C 0 is the stiffness tensor of the undamaged material. This energy is assumed to be the state potential. The damage energy release rate is defined as the variable associated to the damage state variable in the state potential: ( with the rate of dissipated energy: . oplf/ .
For an isotropic damage model, this equation reduces to: The damage energy release rate is always positive and thus the rate of damage must be positive in order to comply with the second principle of thermodynamics. The rules governing the evolution of damage can be defined following the same principles as those in elasto-plasticity. In fact, the general formalism is that of generalised standard materials where the loading function is assumed to be a pseudo potential of dissipation (Lemaitre, 1992). The evolution of damage requires the definition of a loading function: where Y 0 is a parameter which defines the threshold of damage, and Z is the hardeningsoftening controlling variable. The evolution law is prescribed according to the normality rule: . ·if d==A-OY . .
with the Kuhn-Tucker conditions A 2:: 0, f ~ 0, and Af == 0. The evolutionary equations are completed by the definition of a hardening function: where H can be regarded as the equivalent of a hardening-softening modulus in generalised plasticity. In most applications, this modulus is not constant. The major difference with plasticity is that the loading function and the evolution equation for damage are provided in an explicit form because these quantities are function of the total strain.

Damage Model with Crack Closure
Crack closure effects are of importance when the material is subjected to alternated loads, for instance reinforced concrete structures subjected to earthquakes. The fundamental property of the scalar damage model is that the time derivative of damage is constrained to be positive or zero (second principle of thermodynamics). During load reversals, however, micro cracks are closing progressively and the tangent stiffness of the material should increase while damage is constant (Fig. 1 ).

Stiffness of the damage material
Stiffness recovery due to crack closure Material damaged in tension (path OA) and unloaded up to point B Figure 1. Schematic response of a material subjected to uniaxial compression after being damaged in tension.
Within isotropic damage modelling, one solution is to introduce two damage scalars (d 1 ,dc), instead of one, in order to separate the mechanical effect of micro cracking as a function of the type of loading. La Borderie ( 1991) developed the corresponding constitutive relations and applied it to concrete. The elastic energy is: where (a)+ is the positive part of the stress, i.e. the stress tensor expressed in its principal directions with the positive principal stresses only, and (a)-=a-(a-)+.
Inelastic strains due to the growth of damage are also taken into account in this model with a function /(a) introduced in order to represent the progressive vanishing of inelastic strains due to damage in tension when the material is subjected to compression: where ac is the crack closure stress (material parameter). The stress-strain relation reads: Since the model has two damage variables, the evolutionary equations include two yield functions. These functions are constructed on the same principles as the function used for the one scalar model, except that there are two energy release rates associated to the two damage variables respectively: The two loading functions are: (13) and an associated model can be constructed with two independent hardening functions similar to Eq. (8).
Note that this constitutive relation cannot be integrated explicitly. The stress decomposition into positive and negative parts is not known a priori. Therefore, the conditions of evolution of damage are not known for each integration step. Same as in classical plasticity, the integration of the constitutive relations is performed with an iterative scheme: the equations of the system are the state laws and the evolution laws recast as follows in order to obtain a well-conditioned algebraic system: Here, F;, F;; are calculated from the integrated equation of evolution of damage instead of an incremental one. It does not change the general picture of the integration process, however. The system is written in a matrix form: where matrix [J] can be computed and programmed using formal calculus. A full Newton algorithm is implemented in order to solve this non-linear system. Note that there are two sorts of non-linearities: the first one is due to the evolution of damage (i.e. whether damage grows or not), the second is due to crack closure effects. For constant values of the damage variables, the material is orthotropic in the coordinate system of the principal stresses.

Integrated Damage Model
As we saw above, the evolution of damage is very often related to the state of strain. Moreover, it is defined in an explicit, integrated way, which is easier to handle. The damage loading function in Eq. (6) can be redefined as: where E: is a positive equivalent measure of strain and K is a threshold value. The equation f == 0 represents a loading surface in strain space. For the uniaxial tensile case, the equivalent uniaxial strain in Eq. (16) is straightforward. It is the axial strain. However, for general states of stress, damage evolution should be related to some scalar quantity, function of the state of strain. There are, to this regard, several proposals. For example, an appropriate definition for metals is rooted in the elastic stored energy (Peerlings et al. 1998): which is depicted in Figure 2. It is equivalent to the integrated version of the constitutive relation presented in section 2.1. For concrete, Mazars (1984) proposed the following form: where &i are the principal strains. A third possibility, which is also mentioned by Peerlings, is the modified von Mises definition: I 1 and 1 2 are the first invariant of the strain tensor and the second invariant of the deviatoric strain tensor respectively. Only the modified von Mises criterion leads to a new material parameter, namely the factor k. The parameter k is the ratio between uniaxial compressive and uniaxial tensile strengths. These criteria are plotted in figure 2 (with k == 1 0). The loading surfaces  are closed contours around the origin. The dashed lines represent the constant uniaxial compression and uniaxial tension stress paths. The evolution of damage has the same form as in the previous sections: The function h( K) is specific, depending on different models in the literature. For tension only, exponential softening can be used: where K 0 , a, 17 are model parameters.
In order to capture the differences of mechanical responses of the material in tension and in compression, Mazars proposed to split the damage variable into two parts and used the equivalent strain defined in Eq. (18): where dt and de are the damage variables in tension and compression, respectively. They are combined with the weight coefficients at and ac defined as functions of the principal values of the strains 8~ and 8ij , due to positive and negative stresses (see Mazars, 1984).
In uniaxial tension at= 1 and ae =0. In uniaxial compression ae = 1 and at =0. Hence, dt and de can be obtained separately from uniaxial tests. The purpose of exponent f3 is to reduce the effect of damage on the response of the material under shear compared to tension (Pijaudier-Cabot et al. 1991 ). The evolution of damage is provided in an integrated form, as a function of the variable K (Mazars, 1984):  Be =1500, /3=1.

Damage coupled to plasticity
Elastic damage models or elastic plastic laws are not totally sufficient to correctly capture the constitutive behavior of concrete. In some cases, the calculation of the damage variable is a key point. For instance, it can become essential when coupled effects are considered (coupling between damage and permeability, damage and porosity ... ). Picandet et al. (2001) observed in experiments that there is a relationship between damage and the gas permeability of the material. Damage, measured according to the unloading slope during cyclic compressive loading is related to the variation of permeability, independently from the type of concrete tested (normal concrete, high strength concrete, fiber-reinforced concrete). Thus, the capability of the constitutive model to capture the unloading behavior is essential if a proper evaluation of the permeability is required. An elastic damage model is not appropriate as irreversible strains cannot be captured: a zero stress corresponds to a zero strain and the value of the damage is thus overestimated ( figure 4a). An elastic plastic relation is not adapted (even with softening, see for example Grassl et al, 2002) as the unloading curve always follows the same elastic slope ( figure 4b ).
An appealing alternative consists in combining these two approaches into an elastic-plasticdamage law. The softening behavior and the decrease in the elastic modulus are so reproduced by the damage part while the plasticity effect accounts for the irreversible strains. With this formulation, experimental unloading can be correctly simulated ( The plastic model is governed by the following set of equations: with the effective stress defined as where c, c and ~ are respectively the total, elastic and plastic strains, ' A the plastic multiplier and kh the hardening parameter ( 0 ~ kh ~ 1 ). m and h are the flow vectors defined by : fis the yield surface and t; a function of the stress (Etse and Willam, 1994).
with f1 and 1 2 the effective stress invariant and e the Lode angle function of the second and third stress invariant and ranging from [ -n/6 ; n/6]. The evolution of the plastic multiplier is finally given by: This local problem is solved with an iterative procedure associated to a closest point projection algorithm (see Perez-Foguet et al, 2000). Figure 5 shows the evolution of the yield surface with an increasing hardening parameter kh for simple compression p == K 1 .[:/;,~ == K 2 1:, where K 1 and K 2 are two positive constants). As the hardening parameter has a limited value of 1, once it has reached this critical level, hardening is not allowed any more and the yield surface becomes a failure one (constant effective stress). Figure 6 gives the axial stressstrain curve for a loading in tension. As the development of damage is predominant during simple tension tests, the damage and plastic-damage models are almost equivalent. Cyclic compression is a second elementary test. Experimental results are taken from (Sinha et al, 1964). Figure 7 illustrates the numerical response for simple damage law (without any plasticity). With this type of relation, a zero stress corresponds to a zero strain.
-nJ a.. · · · $--:· : . : : : 1 "~~---· · · · · · : The numerical response of the elastic plastic damage model is given in figure 8. Damage simulates the global softening behaviour of concrete while the plastic part reproduces quantitatively the evolution of the irreversible strains. Experimental and numerical unloading slopes are thus similar, contrary to the simple damage formulation response. If this difference could seem negligible, it is in fact essential if a correct value of the damage needs to be captured. The elastic damage model overestimates d whereas the full plastic damage constitutive law provides more acceptable results.  Figure 9 and 10 illustrates the differences between the two approaches tn term of volumetric behaviors.
-51 or' -2.5 10 7 : : The introduction of plasticity associated with the development of damage plays thus a key role in the numerical simulation of a cyclic compression test. Moreover, the volumetric response, that was totally misevaluated by the elastic damage model, is correctly simulated with the full formulation. In figure 11, numerical results are compared with experiment for different levels of confinement pressures in triaxial experiments (P = 0, 1.5, 4.5, 9, 30 and 60 MPa). Simulations and experiment give similar results. The peak position is quantitatively well reproduced (except for 1.5 MP a ).
The global evolution is also correct: the maximum of the axial stress increases with the pressure and the softening part is less and less significant. When the initial hydrostatic pressure takes higher values, the damage part of the model plays a minor role and plasticity effect becomes predominant.

Damage Induced Anisotropy
Microcracking is usually geometrically oriented as a result of the loading history on the material. In tension, microcracks are perpendicular to the tensile stress direction, in compression microcracks open parallel to the compressive stress direction. Although a scalar damage model, which does not account for directionality of damage, might be a sufficient approximation in usual applications, i.e. when tensile failure is expected with a quasi-radial loading path, damage induced anisotropy is required for more complex loading histories. The influence of crack closure is needed in the case of alternated loads: microcracks may close and the effect of damage on the material stiffness disappears. Finally, plastic strains need to be incorporated in the formulation when the material unloads in compression. The following section describes a constitutive relation based on elasto-plastic damage which addresses these issues (Fichant et al. 1999

Principle
The model is based on the approximation of the relationship between the overall stress (latter simply denoted as stress) and the effective stress in the material defined here by the equation: Depending on the interpolation of the damage variable d (n), several forms of damage induced anisotropy can be obtained. Here, it is defined by three scalars in three mutually orthogonal directions. It is the simplest approximation that yields anisotropy of the damaged stiffness of the material. The material is orthotropic with a possibility of rotation of the principal axes of orthotropy.
The stiffness degradation occurs mainly for tensile loads. Hence, the evolution of damage will be indexed on tensile strains. The evolution of damage is controlled by a loading surface f, which is similar to Eq. the new damage surface is the sum of two ellipsoidal surfaces: the one corresponding to the initial damage surface, and the ellipsoid corresponding to the incremental growth of damage.

Coupling with plasticity
Same as for the isotropic model, we decompose the strain increment in an elastic and plastic one. The evolution of the plastic strain is controlled by a yield function which is expressed in term of the effective stress in the undamaged material. We have implemented in this anisotropic damage model the yield function due to Nadai (1950). It is the combination of two Drucker-Prager functions F} and F2 with the same hardening evolution: where q and r are model parameters, w 0 defines the initial reversible domain in the stress space, and p is the effective plastic strain.

Crack closure effects
A decomposition of the stress tensor into a positive and negative part is introduced again: a =(a) + +(a)-, where < a> +, and < a > _ are the positive and negative parts of the stress tensor. The relationship between the stress and the effective stress defined in Eq. (31) of the model is modified: de (n) is a new damage surface which describes the influence of damage on the response of the material in compression. Since this new variable refers to the same physical state of degradation as in tension, de (n) is directly deduced from d(n ). It is defined by the same interpolation as d(n) and along each principal direction i, we have the relation: ( l a d. 1-8 ..
a is a model parameter. The constitutive relations contain 6 parameters in addition to the Young's modulus of the material and the Poisson's ratio. Their determination benefits from the fact that in tension, plasticity is negligible. Once the evolution of damage in tension has been fitted, the remaining parameters are fitted from a compression test. Figure 12 shows a typical uniaxial compression-tension response of the model corresponding to a concrete with a tensile strength of 3 MP a and a compressive strength of 40 MP a.

Non Local Damage
We turn now attention to structural computations with the isotropic damage model. It is well known that when strain softening is encountered (see. Fig. 3 for instance), some difficulties arts e. Physically, after microcracking has taken place in the first part of the inelastic behavior of concrete, a macroscopic mechanism develops in the form of a damage band. The strains tend to localize in a specific area of finite dimension (macrocrack) but the model predicts this dimension to be zero, with a vanishing energy dissipation at failure. The stress response of a material point cannot be described locally but has also to consider its neighbourhood in order to take into account microcracks interaction (Askes, 2000). Mathematically (Benallal et a/, 1993, Peerlings et a/, 1996, when a material point enters the softening regime, the type of the governing equations changes from elliptic to hyperbolic. The mathematical description becomes ill-posed since the initial and boundary conditions that are meaningful for the elliptic problem are usually not relevant for the hyperbolic part. Numerically, the local formulation of the problem yields a response which is dependant on the smoothness and the orientation of the mesh. Let us consider for example a bar with a variable section under tensile loading whose boundary conditions are given in figure 13. In the damage distribution in figure 14 (a), damage localises only in the most loaded element. If we consider the extreme case of an infinitely fine mesh, this geometrical localisation in a single element yields a fracture without energy dissipation. The local formulation thus exhibits a mesh dependency problem and a regularisation technique has to be added in order to achieve a mesh independent result and damage localisation over several finite elements (Fig. 14b ). It is now established that non locality, in a gradient or integral format, is mandatory for a proper, consistent, modelling of fracture (see e.g. de Borst et al. 1993). It avoids the difficulties encountered upon material softening and strain localisation. Within a single approach, it encompasses both crack initiation (for which continuum models are very well fitted) and crack propagation (for which discrete fracture approaches have been developed).
We shall describe in the following integral and gradient non local damage models.

Integral Model
Consider for instance the scalar damage model in which evolution of damage is controlled by the equivalent strain i:~ introduced by Mazars. The :erinciple of nonlocal continuum models with local strains is to replace the equivalent strain & with its average (Pijaudier-Cabot and Bazant, 1987): where n is the volume of the structure, V,.(x) is the representative volume at point x, and lf/(S) is the weight function, for instance: le is the internal length of the non local continuum. 8 replaces the equivalent strain in the evolution of damage. In particular, the loading function becomes f(e,K) = &-K. As we will see in section 3.5, it should be noticed that this model is easy to implement in the context of explicit, total strain models. Its extension to plasticity and to implicit incremental relations is awkward. The local tangent stiffness operator relating incremental strains to incremental stresses becomes non symmetric, and more importantly its bandwidth can be very large due to non local interactions. This is one of the reasons why gradient damage models have become popular over the past few years.

Gradient damage model
A simple method to transform the above non local model to a gradient model is to expand the equivalent strain into Taylor series truncated for instance to the second order: Substitution in Eq. (39) and integration with respect to variable s yields: where c is a parameter which depends on the type of weight function in Eq. (39). Its dimension is m 2 and it can be regarded as the square of an internal length. Substitution of the new expression of the non local equivalent strain in the non local damage model presented above yields a gradient damage model. Computationally, this model is still delicate to implement because it requires higher continuity in the interpolation of the displacement field.
This difficulty can be solved if an implicit format of the gradient damage model is used. Eq. ( 42) is replaced with Here, the definition of the non local equivalent strain is implicit. It is the solution of a Fredholm equation. As shown by Peerlings et al. ( 1996), the implicit form is in fact an exact representation of the integral relation devised by Pijaudier-Cabot and Bazant ( 1987), provided a specific exponential weight function is used.

Extension of the Gradient Approach
In the anisotropic damage models described in section 3.3, the evolution of damage is directional. The evolution of damage is also directional in the microplane-based models and in the smeared crack models. The extension of the gradient damage model to anisotropy, or to any strain-based damage model, can be performed in a very systematic way (Godard, 2001(Godard, , 2003. The Fredholm equation is written for each strain component Enn : Therefore, a non local (gradient type) strain tensor is computed within each iteration. It is from this tensor that an equivalent strain can be computed, or that directional damage growth can be devised depending on the type of damage model that is implemented. This extended form of a gradient model has been implemented in the general purpose finite element code "Code Aster" at EDF.

Computer implementation of damage models
Since the major effect of damage is a secant stiffness reduction when microcracking progresses, a non-linear computational procedure based on the secant material stiffness seems quite appropriate. This type of algorithm is also widely used when an integral non local damage model is implemented. The major reason is that the tangent stiffness operator (in the case of loading) is far from being straightforward because of non locality.  The integral relation due to the non local term is discretised according to the finite element mesh used for the analysis and an usual quadrature rule is employed for its evaluation. In finite element calculations, the weight function is chopped off: the weights that are less than 0.001 are set to zero. The actual volume of integration does not span over the entire volume of 1he solid and the calculation of the integrals requires less computer time and memory as the number of neighbouring integration points is reduced. It should be kept in mind, however, that the secant stiffness algorithm presents serious difficulties, mainly there is no proof of convergence of such an algorithm. In large scale finite element computations, a standard Newton-Raphson algorithm is highly preferable, as convergence is more robust. Such an algorithm relies on the derivation of the consistent tangent stiffness of the material response.
The implementation of the scalar damage model, coupled to the gradient approach is different. The equilibrium equations and the Fredholm equations are solved as a coupled problem. The non local strain tensor can be discretised according to the same finite element grid as the displacements. For a better robustness, a standard Newton-Raphson algorithm is used with a consistent tangent operator that ensures quadratic convergence, if properly implemented.
In order to illustrate results with the integral and gradient damage models, we consider the classical case of a three point bending beam containing a notch. Two different meshes are used (a COARse 890 element and a MEDium 1938 element). Figure 15 shows the beam geometry. Figure 16 presents the comparison between the non local integral and regularised approaches for the medium mesh (se is equal to the square root of c). Choosing the appropriate non local parameters, the two techniques are almost equivalent. Figure 1 7 presents the damage distribution in the half beam for the gradient model. One can see that it is mesh independent.

Model calibration
The calibration of the model parameters in the non local (integral or gradient) damage model is facing the difficulty of getting at the same time the parameters involved in the evolution of damage K 0 , A 1 , B 1 , and the internal length le . The first set is a point-wise piece of information while the influence of the internal length appears in a boundary value problem with non homogeneous strains (by definition, & = & ). It follows that it is not possible to calibrate the model parameters in the damage evolution law, without considering the computation of a complete boundary value problem up to failure. Hence the model parameter calibration has to rely on an inverse analysis technique.
The model parameters are optimised so that computational results are as close as possible to the experimental results. The optimisation process used here is based on a Levenberg-Marquardt algorithm, which performs the minimisation of the following functionnal 3(P): where Rexpi is the experimental response for test i ( i E [1,3 ]), and Ri (P) is the numerical response corresponding to the vector of model parameters P . As we will see next, the calibration cannot be performed on a single set of experimental data, e.g. three point bend experiments. One possibility is to perform the calibration from several experiments, and for instance size effect test data.

Size effect tests
Three point bend experiments have been carried out by Le Bellego et al. (2000) on notched geometrically similar specimens of various height D = 80, 160, and 320 mm, of length L = 4D, and of thickness b = 40 mm kept constant. The length-to-height ratio is LID = 4, the span-toheight ratio is //D = 3, and the notch length is D/1 0 (Fig. 18). These are mortar specimens with a water-cement ratio of 0.4 and a cement-sand ratio 0.46 (all by weight). The maximal sand grain size is da = 3mm. Experimental and computational data will be also interpreted with the help of Bazant's size effect law, in one of its simplest form (Bazant and Planas, 1998). The nominal strength is obtained with the formula: where F is the peak load. The size effect law reads: D 0 is a characteristic size, ;;' is the tensile strength of the material, and B is a geometryrelated parameter. Figure 19 shows the size effect law in a log-log plot. It can be seen that the nominal value of the material strength decreases as the size of the beam increases. This is typical of a damage process in which the process zone is size independent (controlled by the internal length) compared to the structure. For very small sizes, the process zone covers all the specimen and a material strength criterion is expected, while for infinitely large specimens the process zone size is negligible and Linear Elastic Fracture Mechanics should apply.
The calibration of the model parameters from inverse analysis relies first on an accurate finite element representation of each size of specimen. Experience shows that the accuracy of the finite element discretisation in the fracture process zone depends on the ratio of the size of the finite element to the internal length. Therefore, geometrically similar finite element meshes should not be used for the calibration. In order to avoid some possible mesh bias, the finite element size in the process zone should be kept constant and small enough compared to the internal length (at least one third of the internal length). Obviously, the internal length is not known a priori and this requirement is checked at the end of the calibration. The optimisation technique has to start from an initial set of model parameters obtained from a manual (trial and error) calibration. This manual process depends on the exact form of the constitutive relation. Therefore, it is difficult to devise a general technique. Still it can be carried out in a rather systematic way as follow: first the model parameters are fitted such that the load deflection curve of one size of specimen is correctly described (usually the medium size). As an initial guess, the internal length lies between 3d a and 5 d a, where d a is the maximum aggregate size of concrete according to Bazant and Pijaudier-Cabot (1989). The initial damage threshold and the parameters in the equation controlling the damage growth are varied only (in fact, the Young's modulus of the material could be obtained as well from the initial stiffness of the beam). Then, the internal length is adjusted, the other parameters being kept constant, so that the numerical responses for the three specimen sizes are as close as possible to the experimental results. This process requires many calculations, each trial iteration involves three finite element computations. Figure 20 shows such a fit obtained manually. Because there are experimental errors and because the above calibration process is based on a variation of the internal length which is disconnected to the variation of the other model parameters, there are several sets of model parameters which may provide an acceptable fit of the tests data. By acceptable, we mean that the numerical responses are as much as possible within the extreme curves due to experimental dispersion (or within an arbitrary distance from them). The difference between these sets ought to be magnified. For this purpose, the size effect method is quite useful (see Le Bellego et al. 2003    For each set of model parameters, the constants in the size effect law have been fitted and the plot on Fig. 21 has been constructed. It is clear on this figure that none of them are consistent with the experimental data. While the quality of the fit is quite hard to decide upon on loaddeflexion curves, it is much easier to examine it on the size effect plot. In order to exemplify the lack of objectivity of the model calibration (i.e. the lack of uniqueness of the result) when a single size of beam is used as the experimental data set, Fig.  22 shows the results of the automatic calibration on the middle size specimen. Three fixed different values of the internal length have been considered. For each one, the optimisation has been performed yielding different values of the damage parameters. The computed responses are the same, in good agreement with the experimental data. Bending of notched beams is quite often used for parameter identification because it is relatively easy to perform compared to tensile tests that are very sensitive to boundary conditions (Carmeliet, 1999). This figure shows that the model cannot be calibrated from inverse analysis of a single load deflexion curve. Figure 23 shows the result of the optimisation from the three sizes of specimens simultaneously. Compared to the manual fit in Fig. 20, the present set of model parameter provides computational results that are closer to the experiments. We tried to change the initial set of model parameter in order to investigate its stability. Overall, the final sets did not differ very much. The largest difference was found on the internal length, which lied between 40 mm and 50 mm.

Closure
Continuous damage mechanics is a theory which aims at describing the mechanical effect of cracking and void growth in an elastic material. Isotropic and anisotropic damage models have been presented in this chapter. There is a rather generic technique which permits a rational derivation of damage models, given a level of anisotropy which is defined a priori. In order to be more representative of the concrete mechanical response, damage must be coupled to plasticity. A simple, and classical, technique is to assume that the plastic part of the model is written in term of the effective stress. The same applies to the modelling of rock masses.
Strain-softening, which is typical of concrete and rocks, induces spurious strain localisation. Integral or gradient damage model are solutions to circumvent this problem. A central issue is the calibration of the model parameters, and in particular the experimental determination of the internal length. This can be achieved from size effect test data with the help of inverse analysis of structural responses. It is important, however, to underline that a proper calibration cannot be achieved from the fit of a single experiment. Size effect experiments on geometrically similar specimens, or bending and pure tension tests are required.