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 A<sub>t</sub>, B<sub>t</sub>,A<sub>c</sub> and B<sub>c</sub>.

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><b>Figure 1.</b> 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> <b>Figure 2</b>. 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> <b>Figure 3</b>.  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> <b>Figure 4</b>. 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> <b>Figure 5</b>. 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

Comsol implementation of mechanical damage models for concrete

Background
============
Several coupled Hydro-Mechanical-Chemical (HMC) processes are involved during the construction and operation of a heap leach pad. Due to the complexity of the problem, common existing heap leaching models solve the different phenomena in an uncoupled or partially coupled way. Also, often modeling of leaching of copper oxide ores assumes that mineral dissolution/precipitation is described by the diffusion-controlled shrinking core model and the system is composed by a single copper-containing mineral and a single mineral representing the gangue. Although these simplifications help to solve the problem, they could result in poor predictions of the composition of the pregnant solution (PLS) and consequently, an over- or under- estimation of the copper recovery and acid consumption. Further, improving the design of heap leach pads requires performing an integrated analysis of the HMC phenomena. Mechanics stability and mineral recovery efficiency cannot be estimated accurately without considering multiple coupling of the complex chemistry, soil mechanics and hydrodynamics. Here we present a coupled HMC model that accounts for the kinetics of several mineral species and the evolution of a muti-lift heap. Numerical simulations were performed with Comsol-PhreeqC (iCP), an interface that allows solving coupled unsaturated flow, soil mechanics and reactive transport using the Sequential Non Iterative Approach (SNIA).

 

Numerical implementation in iCP
==========================
The HMC governing equations are solved with the numerical interface Comsol-PhreeqC (iCP) (Nardi et al., 2013). iCP is a Java interface that couples the general purpose finite element framework Comsol Multiphysics (www.comsol.com) and the geochemical simulator PHREEQC (Parkhurst and Appelo, 1999; http://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/). In brief, the tool is based on implementing an operator splitting technique for solving the hydro-mechanical problem (a coupled or uncoupled system of partial/ordinary equations and algebraic differential equations) on Comsol and the chemical system of algebraic and differential equations on PHREEQC.

Simulation of a multi-lift heap leach pad
====================================
The present HMC model was applied to simulate the construction and operation of a synthetic copper heap leach pad consisting of five ore layers settled on a terrain with a slope of 2%. Each layer is simulated sequentially, mimicking the sequence of construction and operation of a real heap leach pad. We assumed that the construction of each layer is performed during 10 d, followed by an irrigation period of 90 d at a constant irrigation rate of 7.5 l/h-m2. The stacked ore consists of an artificial mineral mixture made from Chrysocolla (20% w/w), Chalcopyrite (10% w/w), Pyrite (8% w/w), Albite (30% w/w) and K-Feldspar (32% w/w), whose chemical evolution is assumed to be governed by kinetics reactions. Equilibrium precipitation-dissolution reactions of Jarosite-K, Jurbanite, Basaluminite, FeOH3(a) and Schwertmannite were also taken into account.

The figure below shows the evolution of concentrations of selected dissolved species at the output of the drainage system. The evolution of pH shows the arrival of the acidic solution to the liner system in a time of around 60 d (top left). Kinetic dissolution of Chrysocolla and Chalcopyrite are triggered by the infiltration of the acidic solution and, as a result, an increase of Cu concentrations is observed in the outlet flow (top right). The chemical evolution of impurities like Si and Al (bottom) clearly shows that under this low pH conditions kinetic dissolution of alumina-silicates is also favoured. This process is highly dependent on the pH and, over a long period of time, porosity can be significantly modified. As expected, the composition of the PLS responds to the irrigation and drainage rates of the system (shown in grey colour).



<figure className="image">
  <StrapiImage  className="mx-auto"  src="/uploads/copper_heap_leaching_figure_1_024bc9026d.png" alt=""/>
  <figcaption>Evolution of pH (a) and concentrations of Cu (b), Si (c) and Al (d) in the pregnant leach solution.</figcaption>
</figure>


 

Chrysocolla and oxide minerals leach rapidly under heap leach conditions. Chalcopyrite hardly reacts at all in heap leaching. The video below shows that the Chalcopyrite concentration (middle) practically does not change, while the Chrysocolla (top) is almost completely dissolved at the end of the leaching of each layer. The latter has a major impact on the porosity of the ore bed (bottom).

 
<figure className="image">
<iframe width="560" height="315" src="https://www.youtube.com/embed/vYcgOlo9O_o" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
  <figcaption>Evolution of Chrysocolla (top) and Chalcopyrite (mid) concentrations and porosity due to chemical reactions (bottom).</figcaption>
</figure>
 

The copper recovery and acid consumption curves (see the figure below) directly reflect the different dissolution rates of the mineral species considered in the model, specially the Chrysocolla (fast dissolution rate) and the Chalcopyrite (slow dissolution rate). The first layer provides a copper recovery of about 30%, while the copper recovered from the leaching of subsequent layers is about 8%. At the end of the operation (500 d), approximately the 60% of the copper has been recovered. The copper recovering increases due to the leaching of residual mineral from previous layers, depending on the thickness of the loaded layers and the irrigation times. The competition between Chrysocolla and Chalcopyrite is also evidenced by the acid sulphuric consumption. The acid consumption during the leaching of a new layer is progressively lower than the acid consumption of previous leaching stages. This is because the mass of mineral of the new layer is lower and the impact of acid consumption from previous leached layers is minor.

 
<figure className="image">
  <StrapiImage  className="mx-auto"  src="/uploads/copper_heap_leaching_figure_3_20d4aa9811.png" alt=""/>
  <figcaption>Predicted copper recovery (top) and acid consumption (bottom).</figcaption>
</figure>

Changes in porosity and permeability due to mineral dissolution/precipitation and irrigation time can affect significantly the hydromechanical behavior of the heap leach pad. First, the hydrodynamics behavior of the system is conditioned by the time that takes the irrigation solution to saturate the heap as it flows to the drainage system. After a short period of time, a quasi-steady state is achieved for every layer. Secondly, the mechanical effects are traduced in an increase in the vertical stress developed mainly on the rock foundation as the layers are progressively loaded. This vertical stress increase also reaches the top of the third layer, where the vertical displacement reaches a maximum. The evolution of liquid saturation, Darcy velocity field and vertical stress is shown in the next video.

<figure className="image">
<iframe width="560" height="315" src="https://www.youtube.com/embed/iM1oiOpNlOI" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
  <figcaption>Evolution of liquid saturation and Darcy velocity field (top) and vertical stress (bottom) in the leach pad.</figcaption>
</figure>
 

There are zones of the ore where the porosity increases due to the mineral dissolution and other zones where the consolidation of the layers leads to a porosity decrease. In the present application case, chemical predominates over hydromechanical effects. The figure below (top) shows the combined hydromechanical and chemical effects on the global porosity at an observation point (white dot in the previous video). The consolidation of the lifts related to the porosity decrease is illustrated in the figure below (bottom) in terms of the vertical displacement at the cross section A-A´ shown in the previous video. The maximum vertical displacement is obtained at mid height, which seems to be in agreement with field observations.

 
<figure className="image">
  <StrapiImage  className="mx-auto"  src="/uploads/copper_heap_leaching_figure_5_1b37d3a52f.png" alt=""/>
  <figcaption>Evolution of porosity at the centre of the first lift (top) and vertical displacement (bottom) at the cross section A-A’.</figcaption>
</figure>


Conclusions
============
The results show that changes in porosity and permeability due to mineral dissolution/precipitation and irrigation time affect significantly the hydromechanical behavior of a heap leach pad, and vice versa. A strong advantage of the present model is that it can treat multiple chemical species simultaneously, providing solutions that are best for interpreting all the chemical interactions within a leaching system. The stability of the pad and the ore recovery efficiency can be predicted with more accuracy because the model considers the interaction between variable saturated flow, mechanical deformation and chemical reactions during the construction and operation of the heap.

Acknowledgements
============
I would like to express my gratitude and appreciation to Albert Nardi and Mayu Tincopa who laid the foundations for this research. I also want to acknowledge the partners of the iMaGe consortium (www.image-modelling.net) for supporting the iCP development.

References
============
*Nardi, A., Idiart, A., Trinchero, P., de Vries, L. M. and Molinero, J. *(2013). **‘Interface Comsol-PHREEQC (iCP), an efficient numerical framework for the solution of coupled multiphysics and geochemistry’, **Computers & Geosciences, submitted.

*Parkhurst, D. L. and Appelo, C. A. J. *(1999).** ‘User's guide to PHREEQC (version 2)-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations’, **U.S. Geological Survey Water-Resources Investigations, vol Report 99–4259, 312 p.

thumbnail_copper-heap-leaching_overview.png

copper-heap-leaching_overview.png

Heap leach pads constitute a common practice in the mining industry.

Hydro-mechanical-chemical coupled modeling of copper heap leaching with iCP

We are pleased to announce the release of version 2.0 of iCP! iCP 2.0 includes new capabilities such as adding extra dimensions and chemo-mechanical couplings and has been updated for compatibility with the latest releases of Comsol and PhreeqC. 

You can download iCP 2.0 [here](/static/downloads/iCP-v2.0.zip) and you can request license by filling in the [license form](/form/icp-license-request/).

The following list outlines all the new capabilities/improvements added in iCP 2.0: 
- **Extra dimensions** can now be added to your reactive transport models. Extra dimensions extend a standard geometry with additional spatial dimensions and can be used to represent a rock matrix in a fracture-matrix system. The developments carried out in the iCP code allow to include chemical reactions in both domains, the extra dimension representing the matrix and the fracture. The main benefit of this methodology is the possibility to include matrix diffusion in discrete fracture networks (DFNs) of thousands of fractures with complex intersections, see for instance the animation below. It can be challenging or impossible to mesh the fractures and rock matrix with such a complex geometry. In this context, the extra dimensions of COMSOL can be very handy when the matrix has low permeability and transport is diffusion dominated. An example (Example 3) has been added to the examples library to demonstrate how to set up a model with extra dimensions. 

<figure className="image">
      <StrapiImage
        src="/uploads/extradim_dcd8cc3c5a.gif"  className="mx-auto" alt=""
      />
</figure>

- **Mechanical modelling** performed in Comsol can now be coupled with reactive transport with iCP.  These capabilities can be used to model the chemical degradation of concrete and its impact on mechanical properties (see animation below) or to predict the swelling potential of bentonite under different geochemical conditions. Chemo-mechanical coupling is demonstrated with Example 4 in the examples library. 

<figure className="image">
      <StrapiImage
        src="/uploads/mechanics_gif_9d4abcd0ff.gif"  className="mx-auto" alt=""
      />
</figure>

- The amount of exchanger and surface complexes specified in PhreeqC can now be updated from Comsol. This can be used to model variations in number of sites or cation exchange capacity due to mineral dissolution/precipitation. 

- iCP 2.0 is compatible with the latest releases of Comsol: versions 5.4, 5.5 and 5.6. iCP 2.0 by runs Comsol in Standalone mode (without Comsol server) by default, which improves computational performance and memory handling. Note that this option is not available with Comsol 5.5. 

- iCP 2.0 is compatible with Comsol physics that include internal DOFs. This means that physics can now be run with the Compute boundary fluxes option enabled.

- The IPhreeqC library has been updated to the latest version, 3.7.0.

- The default water content used by iCP (when it is not coupled) is now the value specified in the input PhreeqC files. Note that in previous versions this value was overwritten by default by 1 kg of water. 

- Error logging has been improved in iCP 2.0. When an error occurs in the chemical step, iCP now specifies the coordinates at which the error occurred. Also, when an error occurs in the transport step the full error message given by Comsol is displayed.


icp20.png

We are pleased to announce the release of version 2.0 of iCP! iCP 2.0 includes new capabilities such as adding extra dimensions and chemo-mechanical couplings.

Geomechanics

Comsol implementation of mechanical damage models for concrete

Hydro-mechanical-chemical coupled modeling of copper heap leaching with iCP

iCP v2.0 release