Introduction

Damage in a material is evidenced by the reduction of the capacity to resist stresses as loading requirements increase. In order to simulate this phenomenon, different approaches have been proposed, which can be summarized in two large groups. The first corresponds to fracture models, which concentrate the crack initiation and growth process as a discontinuity. In the second, the distributed crack models distribute the damage effects on a specific volume, that is, it uses continuous variables related with the density of these defects to describe the materials degradation.

In this perspective, the Continuum Damage Mechanics (CDM) has been used by several authors (^{Lemaitre 1985}); (^{Oliver et al. 1990}); (^{Rodríguez 2012}); (^{Amirpour 2017}) to simulate the degradation process of the elastic properties of the material. This tool provides constitutive equations and damage evolution equations within the framework of the thermodynamics of irreversible process, the theory of internal state variables, and relevant physical considerations. In this manner, the CDM is a phenomenological model that considers damage indirectly through an internal variable.

Different models can define the internal variable. Scalar damage models are the simplest, because they consider the degradation of the stiffness through a single scalar variable that equally affects all the components of the elastic stiffness tensor, while maintaining the isotropy conditions (^{Kachanov 1958}); (^{Lubliner et al. 1989}); (^{Faria et al. 1998}); (^{Grassl et al. 2013}); (^{Juarez-Luna et al. 2014}). When the material’s directional properties are considered, the damage variable should be defined in the tensorial field, giving rise to the anisotropic models (^{Lee et al. 1997}); (^{Carol et al. 2001}); (^{Pelá et al. 2013}); (^{Wang et al. 2014}). At the same time, in the case of reversible loadings, the direction of the force in the frictional materials can activate the crack by reducing the net section (tensile stress case) or remaining inactive (compressive stress case). The consideration of this differentiated behavior generates a second classification, which regards this phenomenon as unilateral damage models (^{Chaboche 1993}); (^{Martín 2001}); (^{Maimì et al. 2007}); (^{Mazars et al. 2015}). They establish two damage variables, one for tensile damage and another for compressive damage, which evolve independently, and they can be either scalar or tensorial variables.

The need to implement an elastic model with damage, whose definition is a simple one, leads to establish this work’s objective as the formulation, implementation and validation of a scalar damage model applied to three-dimensional solids, based on the Finite Element Method. The simplicity of the proposed model is based on the existence of a single internal scalar variable, whose definition is established from four parameters of the material, two associated to the elastic properties (Elasticity Modulus, Poisson’s Ratio) and two related to the damage properties (Maximum Stress, Fracture Energy). Thus, it just requires the stress-strain test to obtain the necessary parameters.

The implementation of the model is made with the ABAQUS finite element software (^{Simulia D.S. 2010}) through the UMAT subroutine. The validation results from the comparison of the results obtained in the experimental tests of (^{Gopalaratnam and Shah 1985}) and the numerical modelling of (^{Paredes et al. 2006}) for an element subjected to increasing axial tension. Additionally, the model for analyzing the behavior of a grooved beam subjected to simple flexural stresses is used.

Scalar damage model

The scalar damage models are very simple and versatile models (^{Oliver et al. 1990}), since the non-linear behavior of the material is defined from an internal scalar variable called damage variable, d. This variable, introduced by (^{Kachanov 1958}) measures the loss of secant stiffness of the material and it is defined between 0, for the virgin material, and 1, for the totally damaged material.

In this model, the Cauchy stress tensor, expressed in indicial notation, is given by:

1

Where σ_ij is the stress tensor in the real space,

Thus, the factor (1 - d) is a reduction factor that produces the transformation from the non-damaged equivalent space to the real damaged space.

Model general formulation

The Helmholtz free energy density for a scalar degradation process is:

2

Particularly for the proposed case, the free energy of the virgin material is obtained as a quadratic function in Ɛ_ij, according to:

3

For isothermal thermomechanical processes, the thermal problem can be separated from the mechanical problem, by applying the Clausius-Plank inequality (^{Oller 2001}):

4

The application of the Coleman-Noll procedure (^{Gurtin et al. 2010}) guarantees a non-negative energy dissipation,

5

6

7

The constitutive equation shows that, for the proposed scalar damage model, all the components of the effective stress tensor are equally deteriorated; this means that the Poisson’s ratio remains constant (isotropic behavior).

The thermodynamic force associated to the damage variable is the free energy of the non-damaged material,

Consequently, to define the damage model it is necessary to establish: (1) an equivalent stress, (2) a damage criterion and (3) the evolution law of the damage internal variable.

3.1 Equivalent Uniaxial Stress

Based on the effective stress tensors, the equivalent uniaxial stress is defined as a function of the second invariant of the stress deviator tensor, J'₂, as follows:

8

Thus, a complex stress state, as the triaxial, can be reduced to a comparable scalar.

3.2 Damage Criterion

The following functions are proposed for the definition of the damage threshold function and the potential damage function, respectively:

9

10

Where

is a scalar function of the effective stress tensor and c(r, d) is a scalar function that depends on the damage threshold, r, and the damage variable, d. The damage threshold is defined in the closed interval between 0 and the initial damage threshold, r⁰, which is a property of the material.

The damage criterion defines a limit surface that distinguishes between a state of elastic behavior and another one showing the degradation process of the material’s properties.

3.3 Evolution Law of the Damage Variable

Analogously to the incremental plasticity theory (^{Lubliner et al. 1989}), an evolution law is used for the damage variable defined by:

11

Where

is a non-negative scalar called damage consistency parameter, which defines the loading, unloading and reloading conditions through the Kuhn-Tucker conditions:

12

These conditions correspond to problems with unilateral restrictions. When , due to the condition (c), the evolution of the consistency factor,

, will be null. Therefore, it is deduced that the temporary variation of the damage has to be null, d =0 and, thus, the material suffers no degradation and behaves elastically. Otherwise, when

>0, there is a damage evolution, d =

and due to the condition (c),

will be null, which means that the material is in the damage threshold.

The magnitude of the damage consistency parameter is defined based on the damage consistency condition:

13

Considering the temporary variation of the potential damage function:

14

Additionally, the temporary variation of function G can be expressed as:

15

Considering (11),

, and taking into account (14) and (15):

16

It is possible to obtain the evolution law of the damage variable as:

17

Finally, based on equation (7), the temporary variation of the damage variable is:

18

Energy dissipation

Damage models assume that the degradation process in a material point starts when the equivalent stress reaches the damage initiation threshold value, when the damage is activated and produces the degradation of the material’s stiffness. During a uniaxial strain test, as the elongation is increased, stress does too, proportionally, until it reaches the maximum capacity of the material. This behavior follows an elastic evolution so that, in an unloading process, the material recovers its original condition. In case the elongation increases beyond the elastic threshold, the cracking process will be initiated, which will be evidenced by a stress decrease as the strain increases, which is a phenomenon called “softening”. The way in which the softening will evolve may vary and it will depend on the type of material. In turn, the behavior of the material point is associated to the amount of energy that the material is capable of dissipating once the damage is initiated. This energy, called fracture energy (Gf), is an intensive property and strictly positive of the material. The fracture energy and the softening curve are specific parameters of each material and its value is determined through laboratory tests (^{Graffe et al. 2010}).

Linear Softening Case

The evolution law can be specifically defined for the linear softening case. Therefore, the explicit function is established, which defines the evolution of the damage internal variable.

19

Where is a constant that depends on the uniaxial elasticity modulus of the virgin material, E⁰, the damage initiation threshold, r⁰, and the specific energy of the continuous model, gf. This last parameter is obtained as the ratio between the fracture energy, Gf, and a characteristic length, lc, depending on the type and size of the element used.

On the other hand, the threshold damage function can be established in terms of the increase of the equivalent stress.

20

Where is the linear softening modulus that determines the contraction of the elastic domain and depends on the material’s parameters and the solid zone whose behavior is dissipative.

Exponential Softening Case

The explicit evolution law is established, which defines the damage internal variable:

21

Where is a constant depending on the same elastic properties as in the linear softening case.

On the other hand, the threshold damage function is defined based on the equivalent stress as follows:

22

Where is the exponential softening modulus.

Figure 1 shows the explicit evolution laws corresponding to the damage variable and the threshold damage function, for the cases of linear and exponential softening mentioned above.

Figure 1 Explicit evolution laws of the damage variable and damage threshold for linear and exponential softening 

Numerical implementation of the proposed model

The proposed model was implemented in the ABAQUS Platform (^{Simulia D.S. 2010}). Therefore, a calculation algorithm was developed, based on the UMAT subroutine (User-defined Mechanical Material Behavior), which allows defining the mechanical constitutive behavior of the materials that are not incorporated in the software’s internal library.

By each time increase, and for each integration point, the software calls the UMAT subroutine and provides the strain tensor at the beginning of the incremental step and the expected strain increase. The UMAT code should allow determining the tangent constitutive tensor, the real tension tensor and the updating of the variables d and r, at the end of each time incremental step. Therefore, it is necessary to declare the initial Jacobian matrix and the laws that allow determining the evolution of the internal variables.

The implemented algorithm solves the non-linear equations by an implicit process (Backward Euler Method). This method allows determining the size of the damage increase while imposing the condition that the residual force between the equivalent stress and the damage threshold of the previous iteration, is null. By comparing two consecutive iterations, it is possible to determine if the convergence can be achieved within a reasonable number of iterations. If the convergence is not likely to happen, the algorithm adjusts the load increase. Although the implicit integration requires a longer time to solve the equations, it is a stronger and more stable method. Figure 2 presents the flowchart of the calculation algorithm developed.

Figure 2 Flowchart of the integration process of the constitutive model 

Validation

The implemented model is validated based on the comparison of the obtained numerical results and the experimental results taken from the literature. Subsequently, a grid sensitivity analysis is carried out.

5.1 Monotonic Axial Tensile Stress

The experimental results published by (^{Gopalaratnam et al. 1985}) are used in simple concrete specimens for validating the model under direct tensile stresses. Additionally, the results obtained in the numerical model with exponential softening, developed by (^{Paredes et al. 2016}), are incorporated.

The geometry of the case study has a cross section of 19.0 mm x 76.2 mm and a length of 82.6 mm. Table 1 summarizes the properties of the material. The model has a total of 312 linear hexahedral elements of 8-nodes, with a total of 504 nodes. In order to obtain a localized failure, a weakening was introduced in the elements located in the anchoring zone.

Table 1 Properties of the material. 

Figure 3 shows the results obtained by the numerical model proposed for both the linear softening and the exponential softening cases, and it makes a comparison with the experimental and numerical results mentioned above.

Figure 3 Stress-Strain curves obtained from the axial tensile stress models with linear and exponential softening, compared with the reference results. 

The fracture energies summarized in Table 2 are determined from the analysis of the areas below the curves. These results allow establishing the dispersion of each model in relation to the reference experimental value. It can be observed that the model with linear softening presents a better correspondence in the magnitude of the dissipated energy; however, the model with exponential softening has a better approximation to the real dissipative process.

Table 2 Fracture energy dissipated under axial tensile stress 

5.2 Grid Sensitivity Analysis

A grid sensitivity analysis is carried out for the model with linear softening. The study proposed three grid densities constituted by 90, 312 and 748 elements, respectively (Figure 4). Figure 5 shows the Stress-Strain curves for each case of analysis. The grid convergence is evidenced, which ensures independence from the size of the element.

Figure 4 Grid density. Grid composed of 90, 312 and 748 elements, respectively. 

Figure 5 Stress-Strain curve for the grid sensitivity analysis 

Application example: grooved beam flexural strength

With the aim of studying the behavior of the model proposed for three-dimensional solid elements subjected to complex stress state, an application example is carried out, which consists in a beam with grooves along the center, supported on its ends and subjected to bending due to a displacement imposed on the center of the element. Its dimensions are 2.00 m long, 0.20 m high and 0.02 m thick, while the groove measures 0.04 m x 0.10 m x 0.02 m, respectively. Given the symmetry characteristics of the element, only one half is modelled and, therefore, a grid formed by 560 linear hexahedral elements of 8 nodes is obtained, adding up to a total of 951 nodes with a grid densification in the area close to the groove and in the groove itself (see Figure 6). The beam is subjected to a vertical displacement in the center with a maximum value of 0.0015 m.

Table 3 summarize the properties of the material. An exponential-type of damage variable evolution law is adopted.

Figure 6 Left: Front view of the geometry and grid of the modelled section. Right: Perspective view of the grid densification in the grooving area. 

Table 3 Properties of the material. 

Due to the instability of the beam’s behavior as the cracking of the element evolves, the Riks method is adopted in ABAQUS/Standard.

Figure 7 shows the model’s vertical stress-vertical strain curve. It evidences that the structure’s deformation process stops when it reaches a displacement of 1.43 x 10-3 m, corresponding to 96% of the final imposed displacement. This is because the nodes forming the grooving have reached such a damage level that the structure is left without enough flexural stiffness to allow displacement increases, thereby producing convergence problems. Figure 8 shows the deformed structure for the ultimate loading state.

Figure 7 Vertical stress vs. vertical strain curve of the model with exponential softening subjected to simple flexural loading 

Figure 9 and Figure 10 show the degradation process of the grooving section as the imposed vertical displacement evolves. The damage has been calculated as the average of the damage recorded in the nodes of the thickness located equidistant from the neutral axis. The damage process is initiated in the nodes subjected to tensile stress for a vertical displacement of 0.23 mm. This process ends when the nodes located outside of the neutral axis record an average damage of 0.98, both in the tensile and compressive stress zones.

Figure 9 Damage analysis in the grooving section as the applied vertical displacement increases 

Figure 11 shows the damage evolution in terms of the displacement increase for the nodes located in the tensile stress zone. Tthe nodes that are farther away from the neutral axis present a similar behavior with a fast damage evolution, reaching a value of 0.90 for a displacement of 0.50 mm. The nodes located at 25.0 mm and 12.5 mm cause that damage level for the vertical displacements of 0.60 mm and 0.79 mm, respectively. According to Figure 9, it can be concluded that there is a similar behavior in the equidistant nodes in the compressive stress zone. This evidences that for an application of 55% of the ultimate displacement, the section has reduced its stiffness by 90%.

Figure 11 Evolution curve of the damage internal variable as the vertical displacement increases for the nodes located in the grooving section. 

On the other hand, Figure 12 shows the specific stress-strain curve in the direction of axis 1 of the integration points located in the maximum compressive and tensile stress zones of the critical section. A similar behavior is evidenced in both ways of the stress, as expected.

Figure 12 Stress-Strain curve in the direction of axis 1 of the integration points located in the beam’s critical section, on the upper and lower face. 

Conclusions

Scalar damage models adequately represent the behavior of isotropic elastic materials with softening, based on a damage variable. These models are very attractive, due to their great simplicity. This paper presented the formulation, implementation and validation of a scalar damage model applied to three-dimensional solid elements.

It regards a model that depends on a damage variable with the same failure surface, both for stress and compression. These simplifications allow defining the model based on four properties of the material determined in a classic stress-strain test. The application field of the model extends to isotropic elastic materials whose limit tensile stress is equal to the compressive stress, assuming the same softening curve. Since the model presents a single internal variable, the degradation of the constitutive tensor as a result of a strong tensile force is kept constant in case the loading direction is inversed.

In order to distinguish between an elastic behavior state and a state where there is a degradation process of the material’s properties, the model does a transformation from the tension tensor to an equivalent scalar and it compares it with a scalar called damage threshold. In this way, the model has an adequate response to complex loading demand states, such as the triaxial ones. The evolution law of the damage variable can be established from the material’s elastic properties, the damage threshold tension and the fracture energy.

The constitutive equations ruling the problem were implemented in the commercial software ABAQUS, through a user-defined UMAT subroutine. The development code allows determining the tangent constitutive tensor, the real tension tensor and the updating of the state variables at the end of each time incremental step. For the updating of the damage internal variable, an implicit-type of calculation algorithm was developed, based on the Backward Euler method.

The comparison between the numerical results and those obtained from the bibliography account for a suitable correlation among them. The model with exponential softening presents a behavior more adjusted to the one observed experimentally than the model with linear evolution. However, the dispersion of the dissipated energy of the former is 5.32%, against a dispersion of 1.77% of the latter. The grid sensitivity analysis shows that it is possible to achieve satisfactory results with a reduced number of elements. On the other hand, the application example allows arriving to the conclusion that it is a powerful tool aimed at the numerical modelling of isotropic materials with stiffness degradation, and challenges the authors to develop models that consider the directed nature of damage in solid elements under triaxial load demands.