Calibration of nonlocal damage model from size effect tests

The calibration of nonlocal models which contain an internal length has been among the major issues conditioning the implementation of this kind of failure models. Direct calibration from uniaxial testing, where the state of strain remains homogeneous throughout the specimen, is impossible. The softening law is not directly accessible because the strains cannot remain homogeneous during the entire test. In the absence of local information on the displacement field and on micro cracking in the fracture process zone, the calibration has to rely on inverse analysis. This paper presents such a procedure based on three point bend size effect tests on notched specimens. The complete load deflection curves are used for the identification of the constitutive relations. Manual calibration is discussed first. It is emphasised that calibration on the load deflexion curve from a single experiment is not objective. We show that Bazant's size effect law, which is related to peak loads only, may serve as a helpful guide to reach the closest fit. Then, automatic calibration is described. An optimal set of model parameters can be obtained within a reasonable number of iterations.


Introduction
Nonlocal models are necessary in order to achieve a consistent modelling of failure due to strain softening in concrete (Pijaudier-Cabot and Bazant, 1987;Peerlings et al., 1996), metal alloys (Fleck et al., 1994;Leblond et al., 1994;Tvergaard and Needleman, 1994), composites (Geers et al., 1998), and quasi-brittle heterogeneous materials. Hence, their calibration has to rely on sufficiently robust methods for an efficient implementation and an effective use of such models. Whether they are in a gradient or in an integral form, these constitutive relations incorporate an internal length that must be calibrated from experiments in addition to other, more classical, model parameters. The fact that the models contain this internal length raised many questions in the past on the most appropriate technique for measuring it. The early proposal for concrete by Bazant and Pijaudier-Cabot (1989) suggested to compare the response of two specimens subjected to tension. A first specimen was made of concrete and steel fibers. The fibers, glued on concrete, were acting as reinforcements and induced multiple cracking in concrete. In the second specimen, the fibers were not glued on concrete and a single crack propagated. The comparison of 1 the energy dissipated by the two specimens provided an approximate value of this internal length. A recent reappraisal of this experimental technique showed that the accuracy of the results could be very much influenced by friction and interface cracks between steel and concrete (Boudon-Cussac et al., 1999).
As of today, direct experimental techniques that would provide this internal length have not been devised according to our knowledge. Looking at the consequences of the nonlocal approach on the theoretical description of the fracture process, two possibilities based on inverse approaches of calibration, emerged: -In the course of the failure process (e.g., mode I cracking), a fracture process zone forms ahead of the crack. From the modelling point of view, the nonlocal approach is capable of describing this process zone. Its shape is controlled by the strain softening response of the material and by the internal length. Therefore, an accurate measurement of the displacement field within the process zone provides information, which can be compared to a finite element model with the purpose of model calibration. Geers et al. (1996) devised the proper experimental data interpretation for such comparisons and applied it to fiber reinforced composite materials. Manhken and Kuhl (1999) provided also an identification procedure for the parameters in a gradient damage model. -Another technique relies on the fact that the size of the fracture process zone in the material is independent from the size of the structure, provided it does not interfere with its boundaries. Hence, the response of geometrically similar specimens is not geometrically similar and there is a size effect. A complete explanation of this size effect can be found in the textbook by Bazant and Planas (1998). It is due to energy redistribution from the rest of the structure (whose size changes from one specimen to another) where elastic energy is stored, to the fracture process zone whose size is constant. It is only with failure models that contain an internal length that this size effect can be described numerically. It can be the cohesive crack model, in which the characteristic length is related to the length of the fracture process zone (see Bazant and Planas, 1998), or plasticity or damage based nonlocal models (see, e.g., Mazars et al., 1991;de Borst and Gutierrez, 1999) where the characteristic length is proportional to the width of the fracture process zone approximately (Mazars and Pijaudier-Cabot, 1996). Therefore, fits of size effect tests may provide the model parameters in nonlocal models. In this paper, we will focus on this calibration technique based on size effect tests.
The calibration of the model parameters from inverse structural analysis is often a manual trial and error process that can become very demanding in terms of computational effort. Automatic calibration methods have not been devised, except in contributions by Carmeliet (1999) who used a Markovian optimisation, and by Mahnken and Kuhl (1999) who used least square fitting. According to Carmeliet, the experimental data to be fitted had to be derived from tension and three point bend specimens of different sizes. He showed also that it is possible to optimise the set of experiments in order to achieve the best accuracy of the numerical inverse analysis. He did not use, however, the entire experimental information contained in the size effect tests and considered the peak loads only. Additional tensile tests were necessary in order to complete the model calibration. The procedure devised by Mahnken and Kuhl uses the complete experimental information of a panel in tension, but relies on the knowledge of the displacements at several different points of the structure inside and outside the localisation zone. Although this method seems to perform very well on synthetic data, its implementation on real experimental data is still pending. In the present contribution, we will show that the model parameters entering into the description of the nonlinear response of the material in tension, including the internal length, can be obtained from a fit of one set of size effect tests. The load deflexion curves obtained from three point bend notched specimens of three geometrically similar sizes will serve as experimental input data. The nonlocal damage model will be used as an illustrative example. It is recalled in Section 2. Section 3 presents the experimental results used in this study and the sensitivity analysis of the constitutive relations, carried out with the help of Bazant's size effect law. A manual calibration of the model is performed. This calibration serves as an initial guess in an automatic calibration process discussed in Section 4.

Constitutive model
We are going to use the nonlocal version of the scalar damage model for which all the model parameters are recalled in this section. The influence of micro cracking is introduced via a single scalar damage variable d ranging from 0 to 1. The stress strain relation reads: E and ν are the Young's modulus and the Poisson's ratio of the undamaged material. ε ij and σ ij are the strain and stress components, C ij kl and C 0 ij kl are the damaged and initial (elastic) secant stiffness of the material, and δ ij is the Kronecker symbol. The evolution of damage is based on the amount of extension that the material is experiencing during the mechanical loading. In the nonlocal damage model, it is controlled by the averageε of an equivalent strainε defined as: where · + is the Macauley bracket and ε i are the principal strains. Its average is: where Ω is the volume of the structure, V r (x) is the representative volume at point x, and ψ(x − s) is the weight function: l c is the internal length of the nonlocal continuum. The loading function of damage is: κ is the threshold of damage growth. Initially, κ = κ 0 . In the course of loading κ assumes the maximum value of the equivalent strain ever reached during the loading history: h(κ) is the generic expression of the damage evolution law. In order to capture the difference of mechanical responses of the material in tension and in compression, the damage variable is split into two parts: d = α t d t + α c d c , where d t and d c are the damage variables in tension and compression respectively. They are combined with the weight coefficients α t and α c (see Pijaudier-Cabot et al., 1991). The evolution of damage is provided in an integrated form, as a function of the variable κ: Standard values of the model parameters in the damage have been given by Mazars (1984). Same as in the paper by Carmeliet (1999), we will focus attention on the response of the model in tension only. It means that in the course of all the calculations d c = 0 and that the model parameters considered for the identification are the Young's modulus E, the Poisson's ratio ν, the initial damage threshold κ 0 , constants A t and B t , and the internal length l c . The elastic constants are easily obtained from uniaxial compression tests, and will be considered to be known in the present study although the Young's modulus could be obtained from the fit of the experimental data. The calibration should aim at obtaining four parameters: κ 0 , A t , B t , and the internal length l c . In the case of uniaxial tension it is possible to provide a meaning to some of them, as it helps for manual calibration or initial guesses in an automatic calibration process. In a first approximation, κ 0 is the peak strain. 1 − A t is the ratio between the residual stress for large strains and the peak stress. The residual stress is zero when 1 − A t = 0 and when 1 − A t = 1 the uniaxial response looks like an elastic perfectly plastic one. B t controls the softening slope. The larger B t , the faster the stress decreases in the post peak regime.
Note that the general calibration technique discussed in this paper is not dependent on the type of nonlocal model considered. The uniaxial tension responses of most models rely on an average number of four model parameters. Obviously, the sensitivity of each model parameter on the computational results remains strongly dependent on the exact mathematical form of the constitutive relation. The scalar damage model, which is widely used for quasi-brittle materials, is considered as an illustrative example here, and a procedure, same as the one depicted in the next sections, could be equally applied with another type of constitutive relation.

Manual calibration
The calibration of the model parameters in the nonlocal (integral or gradient) damage model is facing the difficulty of getting at the same time the parameters involved in the evolution of damage κ 0 , A t , B t , and the internal length l c . The first set is a point wise piece of information while the influence of the internal length appears in a boundary value problem with nonhomogeneous strains. The model is devised such that in the case of a homogeneous state of strainε =ε. Theoretically, it could be possible to calibrate the evolution of damage on tests in which the strain distribution remains homogeneous over the specimen. This is in fact impossible because whenever strain softening occurs, the strain and damage distributions happen to become strongly nonhomogeneous due to bifurcation and strain localisation. 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. This is true for integral or gradient models based on plasticity or damage theories. With nonlocal models, analytical solutions can hardly be derived. It is the reason why a robust finite element implementation is a crucial stage in the implementation of nonlocal failure models.
We are going first to recall the experimental data from which the calibration process will be carried out. Manual calibration from these data and sensitivity analyses will be discussed afterwards.

Experimental results
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 L/D = 4, the span-to-height ratio is l/D = 3, and the notch length is D/10 (Fig. 1). 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 d a = 3 mm. Fig. 2 shows the load deflection responses for the three different sizes. Experimental and computational data will be interpreted with the help of Bazant's size effect law, in one of its simplest form (see 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, f t is the tensile strength of the material, and B is a geometry-related parameter. In a log-log plot, a strength of material criterion is represented by a horizontal curve and a LEFM (linear elastic fracture mechanics) criterion is depicted by a line with the slope −1/2. The two lines cut each other at the abscissa D/D 0 = 1. D 0 and Bf t are obtained from a linear regression which provides also the fracture energy G f .
where k 0 is a geometry-related parameter. Fig. 3 shows the fit of the size effect law (Eq. (9)) obtained from the experimental data in Fig. 2. The larger the beam, the lower the nominal strength.

Trial and error method of calibration
The calibration of the model parameters from inverse finite element analysis relies first on an accurate finite element model 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. Fig. 4 shows an example of set of finite element meshes. The following computations have been performed with the finite element code CASTEM 2000, in which a secant nonlinear solver is implemented. The detailed calibration process depends on the exact form of the constitutive relation. Therefore, it is difficult to devise a general technique as it relies on the sensitivity and meaning of the model parameters (as provided in Section 2). Still it can be carried out in a rather systematic way as follow: first the model parameters are fitted such that the load deflexion 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 5d 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 (Eq. (7)) 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 in a second phase. The other parameters are 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. When implementing this method in its first phase, the results are observed to be insensitive on the choice of the internal length. Fig. 5 shows some possible fits for the middle size of beam with different characteristic lengths. In the first set, the internal length is l c = 7 mm, in the second it is l c = 13 mm and l c = 40 mm in the third one. Table 1 shows the model parameters for each fit. Several equivalent fits can be obtained for several (fixed) different internal lengths. As we will see in Section 4, it is even possible to find a set of model parameters for each different internal length that provides the same structural response on one size of beam. Therefore, it is emphasised here that fittings of models that contain an internal length cannot rely on a single experiment with a single load deflexion curve. Several possible sets of model parameters provide the same result and calibration is clearly not objective, i.e., it does not yield a unique solution. Therefore, it is quite useful to have a rather accurate first guess of the value of the internal length (e.g., as provided by Bazant and Pijaudier-Cabot (1989)).  Additional information is required in order to sort out the best fit among the available ones. Experiments where the evolution of displacements is known in the course of loading inside and outside the strain and damage localisation zone might provide the correct result (Mahnken and Kuhl, 1999). Obviously, experiments in which displacement fields are directly measured (see, e.g., Geers et al., 1996) are even better because they provide a map of the displacements over the specimen and their evolution in the course of failure, against which finite element results can be compared. The experimental techniques involved in the accurate measurement of displacements remains however quite delicate to implement. Size effect tests are another possible candidate that is investigated in this paper. Because there are experimental errors and because the above calibration process is based on a variation of the internal length, which is disconnected from the variation of the other model parameters, there are still several sets of model parameters that may provide an acceptable fit of the size effect test 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). Fig. 6 shows three of these fits. Note that the rather large variation of load after peak on the computations of the middle and large size beams in Fig. 6(c) is probably the result of a lack of stability of the computed response under displacement controlled conditions. During the considered displacement step, snap back might be observed and the load deflexion curve could be farther off the experimental data set. The fits in Fig. 6 are very much equivalent to each other: the error on the peak loads (for instance) is given in Table 2 for each size. It is in the range of 10-15% for each set. If these fits were equally good, it would again mean that the mechanical model cannot be calibrated objectively because there are several, quite different sets of parameters which provide identical fits. This is, however, not exact and the size effect law can be a useful tool that magnifies the difference between these three sets. As we will see in the next section, the interpretation of the computations with the help of the size effect law reveals better the differences among the three sets of model parameters, which can hardly be distinguished from the analysis of load deflexion curves.

Sensitivity analysis on the size effect law
Before we analyse the different fits obtained in the previous section with the help of the size effect law, we may investigate the sensitivity of the model parameters on the size effect predictions. Starting from a reference set of model parameters, each model parameter has been varied while the others were kept constant. The load deflexion curves for the three sizes of specimen have been computed and then the parameters in the size effect law were obtained from the linear regression in Eq. (10). Because the nonlocal model captures size effect, each set of computed data fits well to the size effect law.
The variations of the model predictions have been plotted in the normalised size effect coordinates (σ/Bf t versus D/D 0 ). Note that for each set of computations (on the three sizes), a different set of parameters Bf t and D 0 is obtained and used in the normalised size effect plot. The difference between each fit is the location of the corresponding points on the size effect curve. For instance and for a given specimen size D, if D 0 changes from one set of parameter to another, the ratio D/D 0 will change Fig. 6. Different manual fits of the model parameters. too and the corresponding point in the normalised size effect plot will be shifted on the horizontal axis. When D 0 decreases, the set of points shifts to the right which indicates that the response of the beams is more "brittle", i.e., closer to a LEFM criterion. Conversely, when D 0 increases, the set of points shifts to the left of the curve, which means that the response of the beams is more "ductile", i.e., closer to a strength of material criterion.   Fig. 7 shows the influence of A t and B t on the size effect plot. The reference set of model parameters is given in Table 3 (except in Fig. 7(b) where B t = 15 000). It can be observed that A t has a very small effect on the size effect plot. When B t increases, the beam responses become more brittle, the set of data is shifted on the right on the size effect plot. κ 0 and the internal length l c have much more influence on size effect as shown in Fig. 8. The reference set of model parameters is again the one given in Table 3 (except in Fig. 8(a) where A t = 0.9, l c = 7 mm). When the damage threshold is increased, or when the internal length is decreased, the beam responses become more brittle. Fig. 9 shows the influence of the internal length on the parameter D 0 in the size effect law. For small internal lengths, less than 40 mm approximately, D 0 increases linearly with increasing internal length. There is a change of regime when the internal length becomes very large (although D 0 is still proportional to l c ). In fact, the fracture process zone is very large when the internal length is very large. It becomes larger that the beam thickness and the beam depth. It is also much larger than the specimen notch and in this range, the simple size effect law implemented here is not applicable because notch sensitivity is entirely lost. Still, it is interesting to keep in mind that the internal length is proportional to the characteristic size D 0 as it may be used for manual calibration.
We may now try to discriminate the various sets of model parameters obtained in the previous section with the help of the size effect law. Same as in the sensitivity analysis, the results of the computations for each set of model parameters served for the determination of constants D 0 and Bf t . Then, the corresponding data points have been plotted in the normalised coordinate system (Fig. 10). We have plotted upper and lower experimental bounds between which the best fit should be in Fig. 10 too. Note that the distance between these bounds is quite large, which means that (i) the test data are rather dispersed even for the largest size, and (ii) that a better fit should be obtained with a fourth size of specimen as suggested by Bazant and Planas (1998).
It becomes clear on this figure that among the sets of model parameters available, the best one corresponds to an internal length equal to 40 mm. It corresponds to the model parameters provided in Table 3 and the model predictions depicted in Fig. 6(c). The three data points are almost within the experimental bounds for this set while it is not the case for the others. It can be also concluded that the optimum fit has not been obtained since it was not possible to place all the points computed with  Even if the model predictions provide rather acceptable results, compared to the experimental data obtained on the specimen sizes tested, their extrapolations to structures of very different sizes will not be accurate because the size effect law is not well captured. In the present contribution, the size effect law is used as a zoom, focusing on the differences between equivalent sets of model parameters obtained from a manual calibration. We compare here experiments and computations according to an empirical formula within the range of specimen sizes tested. Under these conditions, the physical relevance of the size effect law is not really important. Its sensitivity at distinguishing the different sets of finite element results is the important property. Whether the size effect law applies or not is a different issue. Extrapolation to structures of very different sizes requires that the size effect law remains valid. To this respect, it is necessary to underline here that a simplified version of the size effect law has been implemented which holds approximately for a size range of 1 : 15, and for notched structures (Bazant and Planas, 1998). For a very broad size range, it must be generalised (Bazant, 1999(Bazant, , 2001. The counter part of such a generalisation is that the size effect law contains more parameters, and that additional experiments (on different sizes so as to cover the widest range) need to be performed.

Automatic calibration
Although the optimal set of model parameters has not been found, the best one may be used now as an entry within an optimisation procedure. In the present application, the model parameters are optimised so that the three computed load deflexion curves are as close as possible to the experimental results. The optimisation process uses the responses of the three specimen sizes simultaneously. A Levenberg-Marquardt algorithm is implemented. The algorithm performs the minimisation of the following functional ( P ): where R i exp is the experimental response for the size i (i ∈ [1, 3]), and R i ( P ) is the numerical response for size i corresponding to the vector of model parameters P . Several definitions can be used for the responses. Here, we use a simple one: it is the stress computed according to Eq. (8) from the applied load measured for 100 values of the deflexion equally spaced and covering the experimental data range. Because the experimental responses and numerical ones have rather different values from one specimen size to another, it is necessary to express them in a nondimensional way (or to use a weighted functional). This is the reason why the error between the numerical and experimental responses is divided by the maximal experimental response for each size. The functional is minimised within 5 iterations approximately. We observed that the error reaches a minimum and remains constant afterwards. At this stage, the calculation was stopped.
In addition to the load deflexion curves, the size effect law could have been entered in to the functional. There at least two possibilities: one could add on the right-hand side of Eq. (11) the quadratic expression: where ω is a given weight, P max is the maximum load computed according to the finite element calculation (one for each size), and P sel is the maximum load computed according to the size effect law. In addition to the model parameters entering in the constitutive relations, the parameters entering in the size effect law would be considered as unknown and calibrated at the same time (although they are not independent from the previous ones). This is probably a useful enhancement compared to optimisation process used in this paper where the size effect law is calibrated afterwards, under the obvious condition that the test data are placed within the range of validity of the size effect formula. A second possibility is to compare with experiments. The added term in the functional could be where P exp sel is the experimental maximum load. In this solution, the fit would be constrained to follow the experimental size effect law (with the values of D 0 and Bf t determined from the experiments) as closely as possible. In fact it would increase the relative weight of peak load values in the optimisation. This solution has the disadvantage that the size effect formula cannot be used as an additional check of the quality of the optimal set of model parameter because it is included in the calibration process.   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. 11 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 threshold and of the parameters A t and B t . 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. This figure shows that the model cannot be calibrated from inverse analysis of a single load deflexion curve. Fig. 12 shows the comparison between the average test data for the three sizes of specimens and the computations for the optimal fit of model parameters given in Table 4. Same as in Fig. 6, the extreme experimental values have been plotted in this figure. Compared to the manual fits in Fig. 6, the present set of model parameter provides computational results that are closer to the experiments. The optimisation process is rather fast because it starts from a relatively correct initial set of model parameters. 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.
This result is transposed into the size effect plot in Fig. 13. Fig. 13(a) shows that on the size effect normalised plot (same procedure as in Fig. 10), the model prediction is now within the lower and upper bounds. For the sake of completeness, we have plotted in Fig. 13(b) the comparisons between the average experimental results and the numerical predictions, in a diagram, which is now normalised with respect to the parameters of the size effect law identified from the experimental results only. In this plot, the agreement looks much better because the same values of D 0 and of Bf t are used in order to normalise the experimental and numerical data. The normalised size effect plot, each set of data point using its own optimal factors (D 0 , Fig. 13. Size effect predictions for the optimal fit; (a) each set of data is normalised with different parameters obtained from linear regression; (b) plot normalised according to the same set of parameters obtained from a regression with the average experimental data.
Bf t ), acts like a zoom which enlarges the discrepancies between the experimental data and the model predictions as described in Section 3.

Conclusions
The development of nonlocal continuum models with local strains corresponds to the aim of providing sound and computationally robust failure continuum models. Failure is predicted with a finite -nonzero -energy dissipation. The structural size effect, observed for most quasi-brittle materials, is well captured. Their calibration, however, is facing the usual difficulties involved in inverse analysis. We have shown that such models cannot be calibrated from the inverse analysis of a load displacement curve, as it is usually performed for standard (local) models, e.g., from bending tests. Structural scaling of fracture can be a convenient tool, which helps at getting a good calibration. The model parameters can be calibrated first manually from size effect three point bend tests on notched beams. In this trial and error process, the normalised size effect plot is a helpful guideline as it enlarges the difference between the experimental results and the computations. Furthermore, the internal length is proportional to the characteristic size D 0 in the size effect plot, provided that the fracture process zone does not interfere with the boundaries of the solid, which corresponds to large values of internal lengths. Automatic calibration by the minimisation of a quadratic functional and a Levenberg-Marquardt algorithm provides an optimal set of model parameters. The manual calibration is used as the initial starting point of the optimisation, which converges within a small number of iterations. The optimal set of model parameters is consistent with the size effect experimental predictions and is within the experimental dispersion. It should be noted that according to the size effect plot, the dispersion of the experimental results is very large. A smaller experimental dispersion would be obtained if an additional larger size would be tested. Therefore a closer optimal fit could be achieved with the help of size effect tests on at least four different sizes.
Finally, let us point out that the present calibration technique is efficient because the material considered in this study exhibits size effect. A similar size effect is not necessarily observed for ductile materials and a different set of experimental data ought to be devised in this case.