Introduction
=============

Damage mechanics has often been used to study the mechanical behaviour of concrete structures. Damage models can be used to simulate the quasi-brittle behaviour of cementitious materials, since they are capable of representing cracking under various loading conditions. In this post, we present the implementation, regularization and validation of the classical isotropic Mazars’ damage model for concrete in Comsol Multiphysics. We furthermore explore its limitations and present the implementation of an alternative damage model, also proposed by Mazars and co-workers, more suitable for multiaxial loading conditions.

Background
=============

Simulating concrete load bearing capacity and post-peak behaviour is essential to assess the mechanical behaviour of concrete subject to mechanical and environmental loads. The tool presented in this post allows accurate mechanical behaviour characterization through a simple isotropic model. The classical Mazars’ damage model (Mazars, 1986) has been implemented within the Comsol user interface. The regularization of this method to avoid mesh-dependency of the results is carried out by the implicit gradient damage formulation (Peerlings et al., 1996). In addition, we have also implemented a recently developed damage model with updated features (Mazars et al., 2015) to more accurately describe concrete behaviour under biaxial and triaxial loading conditions. The final goal of this work is to develop concrete mechanics models in the Comsol interface that allow their coupling with other processes, such as chemical degradation for durability issues, or moisture and heat transport.

Mazars’ damage model implementation in the Comsol interface
================================================

Classical models following damage mechanics theory (Kachanov, 1958) are based on the definition of a scalar damage variable *d*, varying between 0 for an intact material and 1 for a fully damaged state. This variable affects directly the stiffness tensor, leading to a nonlinear relation between strain (**ε**) and stress (**σ**, MPa) tensors.

One of the most widespread damage models, perhaps due to its simplicity, is the one proposed by Mazars (1986). This model defines an isotropic scalar damage variable described by different evolution laws whether the material is under tensile or compressive conditions. The damage variable *d*, is calculated as a linear combination of compressive and tensile damage contributions based on the evolution of a history variable that stores the maximum value of strains reached. The main model material parameters are At, Bt,Ac and Bc.

In this line, Comsol developed a functionality to access external material models for structural mechanics in a post by Ed Gonzalez (2015). This functionality is intended to enable the programming of any desired constitutive model as an external material model to be used in Comsol applications. In that post, an example is presented with the implementation of Mazars’ damage model. Although this functionality is quite powerful, one of the main drawbacks is the complexity to access the variables defined in the external material code, which complicates the coupling with other physics.

In this post, we present a different implementation of the Mazars’s damage model, in this case using a generic physics (Domain Ordinary Differential Equation) and a specific solver configuration to allow for storing the history variables. The main feature of the implementation is the definition of two history variables needed in the mechanical damage model. Two additional degrees of freedom (DOFs) are added through this new physics.

A compression test with loading/unloading cycles has been modelled with the two above-mentioned approaches. Figure 1 presents the cyclic displacement imposed at the top face of the concrete cylinder and the evolution of the damage variable.

<figure className="image">
<StrapiVideo width="725" height="395" src="/uploads/Mech_Dam_Cyclic_Loading_cce1792b0a.mp4" type="video/mp4"/>
 <figcaption>Figure 1. Imposed displacement in the top face and resulting damage variable evolution with time for the compression test with loading/unloading cycles.</figcaption>
</figure>

Figure 2 shows a comparison of the two implementation methods that confirms that both methods yield the same results. The main advantages of the new implementation are that: (1) it can be fully coupled with other physics and any other constitutive law, (2) all variables are readily available within the Comsol interface, thus facilitating the post-process of results, and (3) the mechanical damage model is easier to adjust or reformulate in the physics tree, without the need of code compilation.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_2_0aec1ae1a0.png" alt=""/>
 <figcaption> Figure 2. Results of the verification test in terms of stress-strain curves using two different damage model implementations.</figcaption>
</figure>

Regularization by the implicit gradient damage formulation
============================================

The use of damage models for simulation of nonlinear stress-strain response including the softening regime has the disadvantage that the results depend on the size of the finite elements, i.e. they are mesh dependent. To remedy this, several methods have been proposed in the literature to regularize the solution so that the results converge to a single solution upon mesh refinement. One of such methods is the so-called gradient-enhanced formulation, which we have implemented in Comsol following the implicit gradient formulation proposed by Peerlings et al. (1996) and Simone (2007). This model consists of a Helmholtz type differential equation in terms of an equivalent strain, . The formulation includes the use of a characteristic length l (m) to compute the nonlocal variable . Thus, in the regularized model, the history variable κ is calculated with this nonlocal variable instead of the equivalent strain. The regularization of Mazars’s damage model has been implemented in Comsol using the built-in Helmholtz equation physics.

To validate this regularized implementation, a classical 3-point bending beam test has been simulated. This kind of test has been widely studied experimentally and computationally by several authors (Kormeling and Reinhardt, 1983; Jirásek, 2011).

Results are presented with and without regularization in Figure 3. Without regularization, the damage model results in different peak loads depending on the finite element size. The model convergence is also sensitive to the element size, resulting in stability problems in the post-peak regions for finer meshes (similar to Jirásek, 2011). On the other hand, with the implicit gradient formulation, the results are much more objective and little mesh dependence is observed. In this case, the characteristic length is 1 mm. The results correspond very well with the regularized test reported by Jirásek (2011) when comparing load-displacement curves. Therefore, it is concluded that the regularized mechanical damage model implementation works properly.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_3_1f0a09af90.png" alt=""/>
 <figcaption> Figure 3. Load (N) vs. displacement (mm) curve of the unnotched three-point beam test with no regularization method (left), and with the regularized damage model (right). Results are shown for different mesh refinements.</figcaption>
</figure>

Mu model
========

Based on Mazars’ classical damage model (Mazars, 1986) a new formulation named “μ model” (Mazars et al., 2015) has been recently proposed. The improvements of this formulation include among others a more accurate representation of concrete behaviour against biaxial and triaxial compression. This is a useful development when modelling materials under multi-axial loading regimes, as is the case of underground structures. In these cases, high compression stresses are normally acting in the three principal directions of the stress tensor.

It is known that the Mazars’ classical damage model under plane stress underestimates the strength of concrete under biaxial compression (e.g. Jirásek, 2011). We have implemented the μ model in Comsol following a similar approach than for the Mazars’ damage model.

The results of the Mazars’ damage model implemented in Comsol under plane stress conditions for different biaxial loads are compared to classical experimental data from Kupfer et al. (1969) in Figure 4. It is clearly observed that under biaxial compression and plane stress conditions the model predicts a significantly lower concrete strength compared to experimental data. On the other hand, strength increase under biaxial compression is obtained with the μ model, as shown in Figure 4 from the current implementation in Comsol.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_4_9a77569630.png" alt=""/>
 <figcaption> Figure 4. Comparison of the results of several biaxial loading tests from Kupfer et al. (1969) and the model results using the Comsol implementation of (1) the Mazars’ damage model and (2) the μ model.</figcaption>
</figure>

To quantify the improvements of the μ model under confined conditions, a set of simulations representing triaxial compression tests were implemented. The triaxial test simulated here were carried out by Ramtani (1990) and were extracted from Mazars et al. (2015). The experiments are carried out at several confinement levels (S= 0, 20, 35, 50 and 100 MPa). The results of the μ model obtained by Mazars et al. (2015) and those extracted from the Comsol implementation are presented in Figure 5. The model implementation is again validated, as results from both models are perfectly matching. Although the implemented scheme gives accurate results in terms of strength, some convergence problems were found after the peak stress with increasing strains.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_5_b4e74e3f30.png" alt=""/>
 <figcaption> Figure 5. Comparison of triaxial test modelled results from Mazars et al. (2015) in blue and the implemented μ model results in Comsol Multiphysics in red.</figcaption>
</figure>

Conclusions
==========

In this post, we presented an alternative implementation of damage models for concrete mechanical behaviour to the external materials functionality of Comsol. We used a generic physics (Domain Ordinary Differential Equation) and a specific solver configuration to allow for storing the history variables of the constitutive model. The presented models can represent concrete mechanical behaviour under diverse loading paths. Models are implemented in Comsol taking special care to history variables description. Classical mesh-dependency problems are overcome by the regularization of the damage variable through a gradient-enhanced scheme. The models were validated by comparison with results from different experimental tests. The implemented models are able to predict concrete mechanical behaviour under confined conditions and can be readily couple with other mechanical processes, as well as chemical degradation or hydraulic and thermal processes.


Acknowledgements
==============

This work has received funding from the European Union's Horizon 2020 Research and Training Programme of the European Atomic Energy Community (H2020-NFRP-2014/2015) under grant agreement n° 662147 (CEBAMA) and from the Swedish Nuclear Fuel and Waste Management Company (SKB), which are gratefully acknowledged. Fruitful discussions with Henrik von Schenck and Per Mårtensson from SKB are also gratefully acknowledged.

References
==============

**González E, 2015.** Accessing External Material Models for Structural Mechanics. COMSOL Blog December 2015, accessible at: https://www.comsol.com/blogs/accessing-external-material-models-for-structural-mechanics/

**Jirásek M, 2011.** Damage and smeared crack models. In Numerical Modeling of Concrete Cracking. Editors: Günter Hofstetter, Günther Meschke, p. 1-49.

**Kachanov L M, 1958.** Time of the rupture process under creep conditions. Isv. Akad. Nauk. SSR, 8, 26–31.

**Kormeling H A, Reinhardt H W, 1983.** Determination of the fracture energy of normal concrete and epoxy modified concrete. Stevin Laboratory, Delft University of Technology, Report No. 5-83-18, 1983.

**Kupfer H, Hilsdorf H K, Rusch H, 1969.** Behavior of concrete under biaxial stresses. Journal of the American Concrete Institute, 66, 656–666.

**Mazars J, 1986.** A description of micro- and macroscale damage of concrete structures. Engineering Fracture Mechanics, 25(5–6), 729-737.

**Mazars J, Hamon F, Grange S, 2015.** A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings. Materials and Structures, 48, 3779–3793.

**Peerlings R H J, de Borst R, Brekelmans W A M, de Vree J H P, 1996.** Gradient-enhanced damage for quasi-brittle materials. Int. J. for Numerical Methods in Engng., 39, 3391-3403.

**Ramtani S, 1990.**Contribution à la modélisation du comportement multiaxial du béton endommagé avec description du caractère unilateral. PhD thesis. University Paris 6, 1990.

**Simone A, 2007. **Explicit and implicit gradient-enhanced damage models. Revue Européenne de Génie Civil, 11(7-8), 1023-1044.

small_comsol_mechanical_damage_2.png

thumbnail_comsol_mechanical_damage_2.png

comsol_mechanical_damage_2.png

Damage mechanics has often been used to study the mechanical behaviour of concrete structures

Solid-mechanics

Comsol implementation of mechanical damage models for concrete