Published in Finite Elements in Analysis and Design Vol. 112, pp. 2639, 2016
doi: 10.1016/j.finel.2015.12.011
Progressive fracture in quasibrittle materials is often treated via strain softening models in continuum damage mechanics. Such constitutive relations favour spurious strain localization and illposedness of boundary value problems. The introduction of nonlocal damage models together with a characteristic length parameter controlling the size of the fracture process zone is known to regularize the problem. In order to account for the nonlocality of these models, it is crucial to work with fine spatial discretizations at the damage progress zone. In this paper we present a nonlocal damage model in combination with a meshadaptive finite element technique that can help automatize the analysis of progressive fracture problems in an efficient manner. Classical twodimensional examples are given to illustrate the presented approach.
Keywords: Quasibrittle materials, continuum damage mechanics, nonlocal damage models, finite element method, meshadaptivity
Continuum damage mechanics is a branch of continuum mechanics that describes the progressive loss of material integrity due to the propagation and coalescence of microcracks, microvoids, and similar defects. These changes in the microstructure lead to an irreversible material degradation, characterized by a loss of stiffness that can be observed on the macroscale.
The term “continuum damage mechanics” was first used by [14], but the concept of damage was introduced by Kachanov in 1958 in the context of creep rupture [16]. In that work Kachanov introduced the concept of effective stress, and by using continuum damage he solved problems related to creep in metals. [31] gave the problem physical meaning by suggesting that the reduction of the sectional area was measured by means of the damage parameter. The thermodynamic formalism involved in the irreversible process of damage was developed by [18]. Other important contributions on damage mechanics include: [22, 33, 26, 24, 25, 8, 9], to name but a few.
The behaviour of brittle or quasibrittle materials such as concrete, rocks, mortar or other geomaterials is particularly difficult to predict. In those cases failure is preceded by a gradual development of a nonlinear fracture process zone and a localization of strain. Realistic failure analysis of such quasibrittle structures requires the consideration of progressive damage due to microcracking, modelled by a constitutive law with strain softening. This typically results in highly nonlinear structural responses and so efficient nonlinear solvers based on arclength control are needed for the numerical simulations [13].
If the damage parameter depends only on the strain state at the point under consideration, numerical simulations exhibit a pathological mesh dependence and the energy consumed by the fracture process tends to zero as the mesh is refined. This is the typical behaviour of the socalled local damage models, which are not able to properly describe both the thickness of localization and the distance between damaged zones. They suffer from mesh sensitivity (for size and alignment) and produce unreliable results. Strains concentrate in one element wide zones and the computed forcedisplacement curves are meshdependent. The reason behind these misbehaviours is that the differential equations of motion change their type (from elliptic to hyperbolic in static problems) and the boundary value problem becomes illposed [2].
Classical constitutive models require an extension in the form of a characteristic length to properly model the thickness of localized zones. Such extension can be done by micropolar [23], strain gradient [34], viscous [19] and nonlocal terms [15]. In our model we have worked with the latter approach using a weighted spatial averaging of the internal variables. In this manner the kinematic and equilibrium equations remain standard, and the notions of stress and strain keep their usual meaning.
The first nonlocal models of this type, proposed in the 1960s, aimed at improving the description of elastic wave dispersions in crystals. Nonlocal elasticity was further developed by [12] who later extended it to nonlocal elastoplasticity [11]. Subsequently, it was found that certain nonlocal formulations could act as efficient localization limiters with a regularizing effect on problems with strain localization [30].
Nonlocal models lead to smooth solutions with a continuous variation of strain. However, to resolve narrow bands of highly localized strains using the finite element method it is necessary to use sufficiently fine computational grids. Fortunately, the mesh must be fine only in the damage progression zone, while the remaining part of the structure can be reasonably well represented by a coarser mesh. In general, the localization pattern is not known in advance, and it is actually tedious to suitably construct refined meshes by hand. Thereby, the efficiency of the analysis can be greatly increased by means of an adaptive mesh refinement technique, which automates the whole process [5,29].
In the present work we present a robust nonlocal isotropic damage model for quasibrittle materials that works in a small deformation regime, along with an adaptivemesh finite element technique that permits adapting the spatial discretization in an optimal manner.
The paper is organized as follows. First, the basic concepts on continuum damage mechanics are introduced. Details are given on the basic components of the isotropic damage theory, and on the equivalent strain forms and damage evolution laws that have been implemented in this work.
Next, we review the regularization technique that has been used to overcome the problems associated to strain localization. The fundamental aspects of the integraltype nonlocal damage model derived are presented, pointing out the most relevant aspects of its numerical implementation. The method for estimating the error of the numerical solution and the meshadaptive scheme are explained in some detail.
Finally, two examples are presented showing that the combination of the nonlocal damage model and the meshadaptive technique is a robust method to model the failure of quasibrittle materials.
The simplest damage model for multiaxial stress states is the isotropic damage model with a simple scalar variable. This model is based on the assumption that the stiffness degradation is isotropic, i.e., the stiffness moduli corresponding to different directions decrease proportionally, independently of the direction of loading. Since an isotropic elastic material is characterized by two independent elastic constrains, a general isotropic damage model should deal with two damage variables. The model with a single variable makes use of the additional assumption that the relative reduction of all the stiffness coefficients is the same, in other words, that the Poisson's ratio is not affected by damage. The stressstrain law is postulated as:

(1) 
where is the total stress tensor, is the total strain tensor, is the effective stress tensor, is the elastic constitutive tensor of the intact material, and is the scalar damage variable.
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. Thereby, we can define the damage variable at a generic section of a material as:

(2) 
where and are respectively the total and the effective area of the section, and is the damaged part of the area. An undamaged material is characterized by . Due to propagation and coalescence of microdefects, the damage variable grows and at the late stages of degradation process it approaches asymptotically the limit value , corresponding to a complete damaged material with effective area reduced to zero.
In order to properly determine the evolution of the damage variable regardless of the loading case we must introduce a loading function specifying the elastic domain and the states at which damage grows. The loading function depends on the strain tensor , and on a variable that controls the evolution of the elastic domain. A typical definition for function is

(3) 
where is the equivalent strain, i.e., a scalar measure of the strain level, and represents a scalar measure of the largest strain level ever reached in the previous deformation history of the material up to its current state, i.e.

(4) 
The above expression implies that , where is the damage threshold, a material parameter that represents the value of equivalent strain at which damage starts. In this formula, is not necessarily the physical time (it can be any monotonically increasing parameter controlling the loading process).
We also postulate the loadingunloading conditions in the KuhnTucker form:

(5) 
The first condition indicates that can never be smaller than , while the second one means that cannot decrease. Finally, according to the third condition, can grow only if the current values of and are equal.
The damage evolution law is defined as:

(6) 
which holds not only during monotonic loading but also during unloading and reloading.
There are various damage governing laws that can be effectively used to model damage growth in quasibrittle materials. In this work we adopt the exponential law proposed in [21], which separates the damage in compression and tension as:

(7) 
with:

(8) 

(9) 
where parameters and are associated to residual strength, and parameters and control the slope of the softening branch at the peak of the stressstrain curve. The coefficient in (7) ranges from 0 to 1 and takes into account the character of the stress state. It is evaluated from:

(10) 
where are the principal strains due to positive effective stresses, i.e., the principal values of .
In the second example of this work, however, we use the simplified version of the previous exponential law as:

(11) 
Concerning the definition of the loading function in (3), we must notice that the expression of the equivalent strain plays a role similar to the yield function in plasticity, because it directly affects the shape of the elastic domain.
Numerous forms of equivalent strain can be found in the literature. In this work, we have used two of them: the socalled [20]:

(12) 
and the modified version of [10]:

(13) 
where is a model parameter that sets the ratio between the uniaxial compressive strength and the uniaxial tensile strength, is the Poisson's ratio, is the first invariant of the strain tensor, and is the second invariant of the deviatoric strain tensor.
(a) Mazars.  (b) Modified von Mises. 
Figure 1: Qualitative damage surfaces in the 2D principal stress space. 
As one can see in Figure 1, the introduced forms of equivalent strain lead to nonsymmetric damage surfaces in which the yield value in compression can be several times the value in tension. This is essential to account for the different strength in tension and compression of geomaterials such as concrete and rocks.
We summarize below the basic ingredients of the isotropic damage model:
Although isotropic damage models are quite simple, they are often used to model the failure of concrete and other quasibrittle materials that show important strain localization. If the damage parameter depends only on the strain state at the point under consideration, the numerical simulations exhibit a pathological mesh dependence, and physically unrealistic results are obtained.
The introduction of a characteristic length into the constitutive model, and the formulation of a nonlocal strainsoftening model, have been shown to prevent the spurious localization of strainsoftening damage, to regularize the boundary value problem, and to ensure numerical convergence to physically meaningful solutions [4].
Integraltype nonlocal models abandon the classical assumption of locality and admit that stress at a certain point depends, not only on the state variables at that point, but also on the distribution of the state variables over the whole body, or over a finite neighbourhood of the point under consideration.
In a general manner, the nonlocal integral approach consists in replacing a certain variable by its nonlocal counterpart, obtained by weighted averaging over a spatial neighbourhood of each point under consideration.
Let be some local field in a domain , the corresponding nonlocal field is defined as:

(14) 
where is the chosen nonlocal weighting function.
In an infinite, isotropic and homogeneous medium, the weighting function depends only on the distance between the source point , and the receiver point . Thereby, we usually write , where is usually chosen as a nonnegative function monotonically decreasing for .
One possible choice for is the Gauss distribution function:

(15) 
where the characteristic length is a material parameter reflecting the internal length of the nonlocal continuum.
If a bounded support is preferred, one can also truncate the previous function as follows:

(16) 
where the interaction radius is a parameter usually related to the characteristic length . In the present work, we have considered .
In essence, the interaction radius represents the smallest distance between points and at which the interaction weight vanishes (for weighting functions with a bounded support) or becomes negligible (for weighting functions with an unbounded support). It represents a very important parameter because it controls the size of the softening region.
The interval, circle, or sphere of radius , centered at , is called the domain of influence of point (Figure 2).
Figure 2: Averaging zone near the boundary of a solid. 
In the application to softening materials, it is often required that the nonlocal operator do not alter a uniform field (consistency of order 0) which means that the weighting function must accomplish the normalizing condition:

(17) 
In order to satisfy (17), the weighting function is rescaled as:

(18) 
A suitable nonlocal damage formulation that restores wellposedness of the boundary value problem is obtained if damage is computed from the nonlocal equivalent strain [3].
In the loading function (3), the local equivalent strain is replaced by its weighted spatial average:

(19) 
The internal state variable is then the largest previously reached value of the nonlocal equivalent strain:

(20) 
It is important to note that the damage variable is evaluated from the nonlocal equivalent strain , whereas the strains used in (1) to compute the stresses are considered as local. This way, during the elastic range, when the damage variable remains equal to zero, the stressstrain relation is fully local.
Numerical implementation of the nonlocal damage model based on averaging of the equivalent strain is relatively simple. The evaluation of the stresses remains explicit, and no internal iteration is needed. All that one needs is to implement the algorithm of weighted spatial averaging and, before damage is evaluated, replace the local equivalent strain by its nonlocal counterpart.
The values of the nonlocal equivalent strain must be traced at the individual Gauss integration points of the finite element mesh, because these are the points at which stresses are computed.
Thereby, the averaging integral in (19) is evaluated numerically. By considering the weighting function defined in (18) we can write:

(21) 
where is a coefficient containing the product of the determinant of the Jacobian and the integration weights of Gauss point , and is the weight of nonlocal interaction between points and , defined as:

(22) 
In the previous two equations, subscript represents the receiver point under consideration, whereas indexes and correspond to source points. Since the chosen weighting function has bounded support (16), vanishes if the distance between points and is larger than the interaction radius . Therefore, the sums in (21) and (22) do not need to be taken over all Gauss points, but only over those that are located inside the domain of influence of point .
In this regard, note that at the damage process zone one must always use an element size smaller than the interaction radius, in order to account for the nonlocal interaction. Otherwise the damage model would become local.
Each Gauss point must have a nonlocal interaction table that gives access to its neighbours. This table must be constructed at the beginning of the problem from a search of nonlocal neighbours. In this work, the search of neighbours is performed by means of a gridbased algorithm. A general rectangular grid is defined in the entire domain and all the integration points are positioned in the cells. This way, the neighbour search that must be performed for each Gauss point is restricted to a limited number of cells, i.e., the ones that fall inside the domain of influence of the considered point (Figure 3).
Figure 3: Gridbased nonlocal search. 
The stress evaluation procedure, repeatedly called during the incrementaliterative strategy, makes use of the nonlocal interaction tables when the nonlocal equivalent strain is computed. To obtain we first compute the local equivalent strains at all Gauss points, and then we calculate the nonlocal counterpart using (21).
If one aims to obtain an iterative solver with quadratic convergence when working with a nonlocal damage model, it is necessary to construct the tangent stiffness matrix in a consistent manner.
Let us start expressing the vector of internal forces as a numerical integration over the Gauss points. In order to make the explanation more understandable, from now on we will be using Voigt notation, i.e.

(23) 
Using the stressstrain law (1) and the classic straindisplacement relation , we expand (23) as follows:

(24) 
The tangent stiffness matrix is obtained by differentiating the internal forces with respect to the vector of nodal displacements . Since the damage variable depends on the nodal displacements through the equivalent strain, we will compute first the derivative of the damage variable with respect to the displacement vector . Taking into account that (6), depends on (20), and depends on through the interpolated strains, we use the chain rule to write the derivative of the damage variable for an integration point as follows:

(25) 
where is the derivative of the damage evolution law with respect to the internal state variable , and is the loadingunloading factor that is 0 in an elastic loading or in an unloading regime, and 1 in a loading regime with growing damage. Thus,

(26) 
Using expression (21), we can differentiate the nonlocal equivalent strain of a Gauss point with respect to the nodal displacements as

(27) 
The derivative of the equivalent strain with respect to the strain vector is a row matrix that will depend on the chosen form of equivalent strain.
Let us now rewrite the expression of the internal forces in (24) as follows:

(28) 
At this point, we can already differentiate (28) and substitute (25) to get a first expression for the tangent stiffness matrix, as

(29) 
Note that the first term in the right hand side (r.h.s.) of (29) is the secant stiffness matrix, that coincides with the tangent matrix in an elastic loading or in an unloading regime (). The second term in the r.h.s. is the nonlocal part of the tangent stiffness matrix. Substituting (27) into (29) yields:

(30) 
Defining for convenience the column matrix , the row matrix , and the coefficient , Eq. (30) can be rewritten like:

(31) 
The double index of the sum, caused by the nonlocal interaction, implies that the term on the right part of Eq. (31) cannot be assembled from element contributions only. Essentially, each pair of Gauss points and contributes to the global stiffness matrix with a block of the same size as that of the classical element stiffness matrix. The difference is that the assembling routine differs from the usual one because in this case one needs to take into account the elements of both points and (Figure 4). Consequently, the global stiffness matrix is always nonsymmetric and its bandwidth increases due to the nonlocal interaction.
Figure 4: Nonlocal assembly process. 
In order to avoid the additional nonzero entries that the nonlocal interaction introduces into the global stiffness matrix, one could neglect the nonlocal terms by using , where is the Kronecker delta. This way, Eq. (31) reduces to

(32) 
where the sum is performed over one index only. Note that the resulting local tangent matrix is no longer consistent, and quadratic convergence is lost.
Probably the most important issue caused by nonlocality is the evolutionary character of the profile of the stiffness matrix. For the simulation of quasibrittle materials like concrete, the consistent stiffness matrix remains local through the elastic branch, and so the initial distribution of nonzero entries is the same as in the local case. However, when the damage threshold is exceeded and the damage zone starts propagating, new nonzero entries appear due to the nonlocal interaction between Gauss points belonging to different elements, and the profile of the stiffness matrix must be dynamically adapted. The number of additional nonzero entries will depend on each particular case, but if a finer mesh is used in the expected softening zones, this number can be relatively high.
The nonlocal averaging certainly increases the computational cost with respect to the corresponding local model. However, the nonlocal model completely removes the pathological sensitivity to the mesh size and partially alleviates the meshinduced directional bias. Consequently, this extra computational effort is indeed worthwhile.
In this work we have used the error estimation and adaptive mesh refinement strategy presented by Oñate and coworkers [5, 6, 2729].
The general algorithm of nonlinear adaptive FEM analysis implemented in this work is described as follows. After reaching the equilibrium state and updating the solution, an error estimation is performed in order to evaluate the error distribution over the mesh. Then, a remeshing criterion uses the information about error distribution and determines the required mesh density. From this analysis, we obtain a new spatial discretization using a mesh generator interface.
In a truly adaptive approach, after generating a new discretization, the data structures corresponding to the newly generated mesh are created, and the transfer of displacements and internal variables from the old mesh to the new one is performed. After the mapping, the internal variables are used together with the strain computed from the mapped displacements to update the internal state of each integration point on the new mesh (to achieve local consistency). Once the transfer has finished, the old discretization is deleted and the mapped configuration is brought into global equilibrium through iteration. Afterwards, the solution continues with the next incrementaliterative step.
Another possibility is to restart the analysis from the initial state after the new discretization is generated. This approach does not require the transfer of the current state from the old discretization to the new one, but from the computational point of view is less efficient than the truly adaptive approach, especially if the remeshing is done frequently.
We present below the main stages of the implemented adaptive procedure:
The error estimator is based on the stress evaluation. We define the error as the difference between an exact value of the effective stresses and an approximate one :

(33) 
Note that we compute the error with the effective stresses and not with the total stresses because the latter tend to zero as damage grows.
A reasonably good value of the exact stresses is obtained from their extrapolation from the Gauss points to the nodes and a posterior smoothing [27]. The approximate value of stresses is simply the default evaluated stresses of the FEM solution. Therefore, the estimated error can be defined as the difference between the smoothed effective stresses , and the calculated effective stresses :

(34) 
We will work with integral measures of the error. The energy norm of the error over an element is defined as:

(35) 
The square of the global error can be obtained from the sum of the squares of all the elemental errors:

(36) 
The smaller is the distance between the nodes of the mesh, the smaller will be the difference between the smoothed and nonsmoothed stresses. Thereby, the presented energy norm tends to zero as the size of the element diminishes, i.e., where is the element size and is the order of the polynomial defining the shape functions.
Determining whether a mesh must be refined or not requires to previously define some quality conditions based on the estimated error.
First, the energy norm of the global error should be smaller than a certain percentage of the deformation energy:

(37) 
where is the permissible percentage of global error, prescribed a priori, and the deformation energy is obtained from:

(38) 
with

(39) 
To determine whether condition (37) is fulfilled, it is convenient to work with a global error parameter defined as:

(40) 
Thereby, when the global error condition is perfectly fulfilled. For the mesh should be refined, and for the average mesh size could be larger.
Taking into account that , the new element size can be obtained with:

(41) 
If one aims at obtaining a selective refinement method, apart from (37), another condition concerning the error of each element must be simultaneously imposed. In this paper, we have chosen a remeshing criterion based on a global error equidistribution [28].
This criterion distributes the global error uniformly among all the elements of the mesh, and so the elemental error must accomplish the following condition:

(42) 
where is the number of elements of the mesh.
Like for the global error, we can work with a local error parameter:

(43) 
with the same meaning that in the global case: indicates an optimal element size, whereas and imply that the element is too large and too small, respectively.
This local parameter allows us to define the new element size that accomplishes (42):

(44) 
where is the space dimension of the problem.
In the end, the remeshing criterion results from the combination of the global error condition (37) and the local one (42). In consequence, the final refinement parameter of the element can be defined as:

(45) 
And the new element size is obtained with:

(46) 
where

(47) 
The mesh generator interface is the preprocessing and postprocessing GiD software used to run the simulations [1]. Therefore, not only GiD meshes the geometry at the beginning of the numerical solution, but it is also GiD which allows to obtain the new spatial discretization every time we need to adapt the mesh.
Thereby, after the error estimator and the remeshing criterion are applied, a "background mesh" file is generated with the information of the new element sizes. Then, GiD allows us to load this background mesh file along with the original mesh file so as to generate a new mesh according to the refinement parameter of each old element.
Before one can continue with the solution of a problem, there is a last stage that must be performed in a truly adaptive approach: the mapping of primary unknowns and internal variables.
In essence, if one aims at continuing the analysis from the current state, instead of restarting it from scratch after every mesh refinement, it is necessary to apply some transfer algorithms for the displacements and internal history variables (in the present case, the state variable governing the damage evolution).
Mapping of the primary unknowns (nodal displacements) is achieved by using the shape function projections. To do so, we must first place each new node inside an old element, by means of a gridbased search, and then we interpolate the displacements of the old nodes to the new ones (Figure 5a).
Mapping of the internal state variables is done through a weighted spatial averaging, similar to the one used for the computation of the nonlocal equivalent strain. The difference is that, in this case, the source points are the integration points of the old mesh, and the receiver points are the integration points of the new mesh (Figure 5b). Once again, another gridbased search is performed in order to determine the Gauss points of the old mesh that fall inside the interaction radius of each new Gauss point.
(a) Nodal displacements mapping. 
(b) Internal state variables mapping. 
Figure 5: Mapping of variables. 
In this section we will present two examples of application of the nonlocal damage model and the mesh adaptive technique proposed. The examples are the threepoint bending test and the fourpoint shear test. For each case we will first assess the robustness of the nonlocal damage model when working with different spatial discretizations, and then we will solve the problem again applying the meshadaptive technique, in order to point out the strengths and limitations of the code in its current state.
Regarding the solving strategy, we will be following the equilibrium path of the problem for a series of steps through a global arclength method in an incrementaliterative process. Self weight will not be considered.
This test is performed with a notched beam subjected to threepoint bending (TPB). The beam has a square cross section of 40 320 mm, a span of 1280 mm. The notch is 3 mm thick and extends over one tenth of the beam depth (Figure 6).
Figure 6: TPB test. Problem statement. Distances in mm 
Plane stress conditions have been assumed. The geometry was meshed by 4noded quadrilaterals with 2 2 integration points.
As it has been stated before, the nonlocal damage model was defined with the Mazars model, regarding the equivalent strain (12) and the damage evolution law (7). The material parameters were obtained from [17] and are summarized in Table 1.
Parameter  Value 
Young modulus ()  38500 MPa 
Poisson's ratio ()  0.24 
Damage threshold ()  
Parameter in compression ()  1.25 
Parameter in compression ()  1000 
Parameter in tension ()  0.95 
Parameter in tension ()  9000 
Characteristic length ()  40 mm 
In order to assess the robustness of the nonlocal damage model, we have also solved the problem with a local damage model based on the damage evolution law of [24] and a modified definition of the equivalent strain presented by Simo and Ju in [33]. The parameters for this model are shown in Table 2.
Parameter  Value 
Young modulus ()  38500 MPa 
Poisson's ratio ()  0.24 
Compressive strength ()  45 MPa 
Tensile strength ()  3.8 MPa 
Fracture energy ()  100 
Limit fracture length ()  5 mm 
The problem was solved for different spatial discretizations. In this case we used three unstructured meshes of 4noded quadrilaterals with a minimum size of 15 mm, 7 mm, and 3 mm, respectively (Figure 7).
(a) Mesh 1: 2024 4noded elements.  
(b) Mesh 2: 2679 4noded elements.  
(c) Mesh 3: 6543 4noded elements.  
Figure 7: TPB test. Spatial discretizations. 
Figure 8 shows the relation between the applied load and the vertical deflection of the beam for each mesh. As one can see from the discontinued equilibrium curves, we had serious difficulties in tracing the response of a full test. The reason behind the convergence problems could be the use of a too global arclength method. Indeed, to account for the localized nature of quasibrittle failure, a more specific control parameter, like the Crack Mouth Opening Displacement (CMOD), could help improving the convergence near snapback zones.
That aside, if we look at the curves in Figure 8a we can see that, although there is no relevant difference between the response obtained with meshes 1 and 2, the peak load actually decreases with the finest mesh 3, and so the total dissipated energy. On the other hand, the response diagram in Figure 8b shows practically the same peak load for the three meshes.
(a) Local damage model. 
(b) Nonlocal damage model. 
Figure 8: TPB test. Forcevertical deflection diagrams. 
Figure 9 shows clearly the different behaviour of the local and nonlocal approaches. While in the local model the thickness of the damage pattern is of the order of the element size, in the nonlocal case the damage pattern is controlled by the interaction radius and so the response remains virtually unaltered regardless of the mesh.
Furthermore, in order to check whether the implemented model can properly reproduce the behaviour of the real test, we have compared the results obtained in the nonlocal model with experiments performed in [17].
Figure 10 shows the forcevertical deflection diagram of the experimental solution and the computed solution with mesh 3.
Figure 10: TPB test. Nonlocal model validation. 
As we can see, the peak load is properly captured, but the post peak branch of the numerical solution falls faster than in the experimental solution, and even shows a certain amount of snapback behaviour. Looking at the elastic branch of the responses, it seems that the stiffness degradation starts before in the computed solution and, in consequence, the peak is slightly displaced to the right. Therefore, we can say that the behaviour of the numerical model seems more brittle than that observed in the experiment. This difference could be due to the approach used in the laboratory test to control the postpeak behaviour.
The same problem was solved again, using the meshadaptive technique described along with the nonlocal damage model.
As mentioned in Section 5.2, we use a remeshing criterion based on the global error equidistribution. The objective of the example was to obtain efficient spatial discretizations for the changing damage states, limiting the global error to a 10% of the deformation energy . The initial mesh is an unstructured mesh of 4noded quadrilaterals with an average size of 10 mm.
In this case three refinements have been performed. In order to understand the meshadaptive process, it is interesting to see the evolution of damage for the different meshes (Figure 11).
A quantitative measure of the suitability of the mesh is given by the global refinement parameter . A value close to 1 means that the global error is similar to the goal percentage of error . The value of 1.33 in Figure 11a denotes that the mesh must be refined in those areas of the structure with higher stress variations, i.e., the zone near the notch, the point of application of the load, and the supports. On the other hand, since this procedure also optimizes the mesh enlarging the elements on those zones with smoother stresses, the number of elements can be reduced from 4453 to 2977.
In the subsequent meshes (Figures 11b, 11c and 11d) the global refinement parameters are smaller than 1, showing that the meshadaptive technique properly keeps the error below the prescribed threshold.
We can say that the main advantage of the adaptive mesh technique is that it automatically optimizes the mesh for the problem we are solving, so that the number of elements can be drastically reduced, and so the computational cost of the problem. However, this procedure still does not distinguish compressive stresses from tensile stresses when estimating the error. In materials like concrete, this can lead to excessively small elements in compressed zones, and so it could be convenient to give more weight to the tensile stresses.
The test will be performed with a singleedge notched beam subjected to fourpoint shear (FPS). The analysed beam has a square cross section of 100 200 mm, a span of 840 mm, and the notch is 10 mm thick and 40 mm depth (Figure 12).
Figure 12: FPS test. Problem statement. Distances in mm 
Plane stress conditions have been assumed. The geometry was meshed using standard linear 3noded triangles with one integration point.
The nonlocal damage model is defined with the equivalent strain form of the modified von Mises model (13), and the simplified exponential damage evolution law (11). The material parameters for the nonlocal approach have been obtained from [32] and are summarized in Table 3.
Parameter  Value 
Young modulus ()  28000 MPa 
Poisson's ratio ()  0.1 
Damage threshold ()  
Compressive to tensile strength ratio ()  10 
Parameter ()  0.8 
Parameter ()  9000 
Characteristic length ()  10 mm 
Like in the first example, we have compared the nonlocal damage model with a local damage model (Table 4).
Parameter  Value 
Young modulus ()  28000 MPa 
Poisson's ratio ()  0.1 
Compressive strength ()  35 MPa 
Tensile strength ()  3.2 MPa 
Fracture energy ()  140 
Limit fracture length ()  5 mm 
In this case, we used three unstructured meshes of 3noded triangles with a minimum size of 8 mm, 5 mm, and 3 mm, respectively (Figure 13).
(a) Mesh 1: 2216 3noded elements. 
(b) Mesh 2: 3502 3noded elements. 
(c) Mesh 3: 7183 3noded elements. 
Figure 13: FPS test. Spatial discretizations. 
In Figure 14 we represent the relation between the applied load and the Crack Mouth Sliding Displacement (CMSD) for each mesh. The curves in Figure 14a show that in this case the peak and the dissipated energy also decrease as the mesh is refined. However, this reduction is not so clear as in the previous example. On the other hand, Figure 14b shows virtually the same peak load for the three meshes. Also, although we can notice some oscillations in the postpeak regions of meshes 1 and 2, we can say that the residual force at the right part of the graph seems to be pretty similar for all cases.
(a) Local damage model. 
(b) Nonlocal damage model. 
Figure 14: FPS test. ForceCrack Mouth Sliding Displacement diagrams. 
Figure 15 shows the damage progression for both models stressing the more localized nature of the local damage model. Regarding the nonlocal model, if we compare this example with the threepoint bending test, we can clearly see that the damage progression zone here is restricted to a narrower region than before. The reason is that now the characteristic length defining the interaction radius is quite smaller than before.
Figure 16 shows the relation between the applied load and the Crack Mouth Sliding Displacement of a reference experimental solution [7] and the computed solution with the finest mesh.
Figure 16: FPS test. Nonlocal model validation. 
We can see that the obtained response is very similar to the experimental solution. What should be noticed, though, is that the failure obtained in the computed solution is slightly more brittle as compared to the experimental one. This is possibly consequence of an imprecise tracing of the equilibrium path, caused by the global arclength strategy used.
This example is performed with the same permissible percentage of global error as in the previous one , but the initial mesh is an unstructured mesh of triangles with an average size of 7 mm.
In this case the number of refinements is four. The evolution of damage for the different spatial discretizations is represented in Figure 17.
Figure 17a shows no damage before the first refinement. It can be deduced then that the subsequent spatial discretization is mainly caused by the high compressive stresses at the loading points and at the supports.
However, from Figures 17c, 17d and 17e it is clear that the adaptive procedure properly captures the damage process zone and adapts the mesh accordingly.
Another aspect to consider here is directly related to the abrupt failure of the computed solution seen in Figure 16. Indeed, when damage grows abruptly the global refinement parameter increases up to a value of 2.20 (Figure 17c). Consequently, the adaptive technique refines excessively all the affected zone, leading to a considerably large number of elements and a much smaller global refinement parameter (Figure 17d). Nonetheless, once damage progress stabilises, the adaptive procedure readjusts the spatial discretization and the resulting mesh is actually efficient (Figure 17e).
Finally, we can notice that the damage pattern in Figure 17e is a bit wider than the one shown in Figure 15f. The reason could be that some extra damage diffusion is introduced due to the errors appearing during the transfer of internal variables between meshes. Therefore, in the near future it would be important to assess the effect of these errors in the traced equilibrium path, and improve the transfer algorithms if necessary.
The implemented nonlocal damage model has shown a robust behaviour for changing spatial discretizations, avoiding pathological mesh sensitivity and meshinduced directional bias. It has also proved to properly capture the experimental peak load of the response diagrams in the two examples analysed: the threepoint bending test and the fourpoint shear test. However, the numerical solution has shown a more brittle postpeak branch as compared to the experimental results. An important aspect to note about this integraltype nonlocal approach is that the averaging performed to account for the nonlocal interaction, considerably increases the computational cost of the solution, as compared to the classical local approach, and also modifies the traditional assembly process of the global tangent stiffness matrix.
That aside, it must be stated that the global arclength strategy used in the simulations has shown an irregular success when tracing the nonlinear solution, and so it is probably not the best approach for the kind of problems we have been solving. Indeed, the implementation of an advanced arclength method with a specific control parameter (CMOD or CMSD) that accounts for the localized nature of quasibrittle materials could help improving the accuracy of the converged solution.
The implemented meshadaptive technique allows one to capture the process damage zone and adjust the spatial discretization accordingly. This procedure generates large elements in zones with smooth stress distribution, and small elements in zones with strong variations of stresses. Hence, efficient mesh distributions can be obtained. Nevertheless, the mesh adaptive technique is still a recent implementation and there are some features that will be improved in subsequent work. For instance, the procedure does not distinguish compressive stresses from tensile stresses when estimating the error, which can result in excessively fine meshes at compressed zones that are not actually transcendent for the damage progression of materials like concrete. Furthermore, errors in the transfer operations of internal variables introduce some extra damage diffusion that should be analysed in order to assess its influence in the traced equilibrium path.
The authors are thankful to Mr. Miquel Portabella for his collaboration in the implementation of the meshadaptive technique.
This research was partially funded by the Advanced Grant project SAFECON of the European Research Council and project NUMEXAS of the FP7 of the European Commission.
[1] GiD the personal pre and post processor. http://www.gidhome.com/, April 2015.
[2] Z. P. Bazant, T. B. Belytschko, and T.P. Chang. Continuum theory for strainsoftening. Journal of Engineering Mechanics, 110(12):1666 1692, 1984.
[3] Z. P. Bazant and M. Jirásek. Nonlocal integral formulations of plasticity and damage: survey of progress. Journal of Engineering Mechanics, 128(11):11191149, 2002.
[4] Z. P. Bazant and G. PijaudierCabot. Nonlocal continuum damage, localization instability and convergence. Journal of applied mechanics, 55(2):287293, 1988
[5] G. Bugeda and E. Oñate. Adaptive mesh refinement techniques for aerodynamic problems. In Proceedings of the International Conference on Métodos Numéricos en Ingeniería y Ciencias Aplicadas. Universidad de Concepción, Chile. 1620 November 1992, pages 513522. Published by CIMNE, Barcelona, 1992.
[6] G. Bugeda and E. Oñate. A methodology for adaptive mesh refinement in optimum shape design problems. Computing Systems in Engineering, 5(1):91102, 1994.
[7] A. Carpinteri, S. Valente, G. Ferrara, and G. Melchiorrl. Is mode ii fracture energy a real material property? Computers & structures, 48(3):397413, 1993.
[8] M. Cervera and M. Chiumenti. Mesh objective tensile cracking via a local continuum damage model and a crack tracking technique. Computer Methods in Applied Mechanics and Engineering, 196(1):304320, 2006.
[9] M. Cervera, M. Chiumenti, and C. A. de Saracibar. Shear band localization via local j 2 continuum damage mechanics. Computer Methods in Applied Mechanics and Engineering, 193(9):849880, 2004.
[10] J. De Vree, W. Brekelmans, and M. Van Gils. Comparison of nonlocal approaches in continuum damage mechanics. Computers & Structures, 55(4):581588, 1995.
[11] A. C. Eringen. On nonlocal plasticity. International Journal of Engineering Science, 19(12):14611474, 1981.
[12] A. C. Eringen and D. Edelen. On nonlocal elasticity. International Journal of Engineering Science, 10(3):233248, 1972
[13] M. Fafard and B. Massicotte. Geometrical interpretation of the arclength method. Computers & structures, 46(4):603615, 1993
[14] J. Hult. Creep in continua and structures. In Topics in applied continuum mechanics, 137155, 1974.
[15] M. Jirásek and S. Rolshoven. Comparison of integraltype nonlocal plasticity models for strainsoftening materials. International Journal of Engineering Science, 41(13):15531602, 2003.
[16] L. Kachanov. On the time of fracture under conditions of creep. Izv. AN SSSR, Otd. Tekh. Nauk, (8):2635, 1958.
[17] C. Le Bellego and É. normale supérieure de Cachan. Couplages chimiemécanique dans les structures en béton attaquées par l'eau: Etude expérimentale et analyse numérique. LMTENS Cachan, 2001.
[18] J. Lemaitre and J. Chaboche. Mécanique des matériaux solides, 1985. Dunod, Paris.
[19] B. Loret and J. H. Prevost. Dynamic strain localization in elasto(visco) plastic solids, part 1. general formulation and onedimensional examples. Computer Methods in Applied Mechanics and Engineering, 83(3):247273, 1990
[20] J. Mazars. Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure. PhD thesis, 1984.
[21] J. Mazars. A description of microand macroscale damage of concrete structures. Engineering Fracture Mechanics, 25(5):729737, 1986.
[22] J. Mazars and G. PijaudierCabot. From damage to fracture mechanics and conversely: a combined approach. International Journal of Solids and Structures, 33(20):33273342, 1996.
[23] H. Mühlhaus. Shearband analysis in antigranulocytes material by cosserattheory. Ingenieur Archiv, 56(5):389399, 1986.
[24] J. Oliver, M. Cervera, S. Oller, and J. Lubliner. Isotropic damage models and smeared crack analysis of concrete.In Second International Conference on Computer Aided Analysis And Design of Concrete Structures, volume 2, 945958, 1990
[25] J. Oliver, M. Cervera, S. Oller, and J. Lubliner. A simple damage model for concrete, including long term effects. In Second International Conference on Computer Aided Analysis And Design of Concrete Structures, volume 2, 945958, 1990
[26] S. Oller, E. Oñate, J. Oliver, and J. Lubliner. Finite element nonlinear analysis of concrete structures using a ‘’plasticdamage model’’. Engineering Fracture Mechanics, 35(1):219231, 1990.
[27] E. Oñate. Structural Analysis with the Finite Element Method. Linear Statics: Volume 1: Basis and Solids, volume 1. Springer Science & Business Media, 2009.
[28] E. Oñate and G. Bugeda. A study of mesh optimality criteria in adaptive finite element analysis. Engineering Computations, 10(4):307321, 1993.
[29] E. Oñate and J. Castro. Adaptive mesh refinement techniques for structural problems. In The finite element method in the 1990's, 133145. Springer, 1991.
[30] G. PijaudierCabot and Z. P. Bazant. Nonlocal damage theory. Journal of Engineering Mechanics, 113(10):15121533, 1987
[31] Y. Rabotnov. Damage from creep. Zhurn. Prikl. Mekh. Tekhn. Phys, 2:113123, 1963
[32] A. Rodríguez Ferran, I. Arbós, A. Huerta, et al. Adaptive analysis based on error estimation for nonlocal damage models. 2001.
[33] J. Simo and J. Ju. Strainand stressbased continuum damage modelsi. formulation. International journal of solids and structures, 23(7):821840, 1987.
[34] H. Zbib and E. Aifantis. A gradientdependent ow theory of plasticity: application to metal and soil instabilities. Applied Mechanics Reviews, 42(11S):S295S304, 1989.
Published on 11/07/19
DOI: 10.1016/j.finel.2015.12.011
Licence: CC BYNCSA license