# Compaction and Tensile Damage in Concrete: Constitutive Modelling and Application to Dynamics

Nicolas Burlion, Fabrice Gatuingt, Gilles Pijaudier-Cabot, Laurent Daudeville

## ► To cite this version:

Nicolas Burlion, Fabrice Gatuingt, Gilles Pijaudier-Cabot, Laurent Daudeville. Compaction and Tensile Damage in Concrete: Constitutive Modelling and Application to Dynamics. Computer Methods in Applied Mechanics and Engineering, 2000, 183 (3-4), pp.291-308. 10.1016/S0045-7825(99)00223-6. hal-01005864

HAL Id: hal-01005864

https://hal.science/hal-01005864

Submitted on 24 Feb 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Compaction and tensile damage in concrete: constitutive
modelling and application to dynamics

Nicolas Burlion, Fabrice Gatuingt, Gilles Pijaudier-Cabot, Laurent Daudeville LMT-

Cachan, ENS Cachan/CNRS/Université P. et M. Curie, 61 avenue du Président Wilson, F-94235 Cachan Cedex, France

The objective of this study is to develop a model for concrete with an emphasis on tension and compaction. Compaction of concrete is physically a collapse of the material voids. It produces plastic strains in the material and, at the same time, an increase of the bulk modulus. The model is based on mechanics of porous materials, damage and plasticity. The computational implementation has been carried out in the Lagrangian finite element code DYNA3D. In order to show the influence of compaction, simulations of a split Hopkinson test performed on confined concrete and on a concrete rod submitted to an impact have been carried out. The examples demonstrate the importance of compaction during an impact, which has a tendency of strengthening the concrete structure.

# 1. Introduction

In a concrete structure subjected to a shock or to an impact, for instance a concrete slab, the material is subjected to various states of stresses which yield different failure modes. Near the impactor, severe hydrostatic compression is observed. This state of stress produces irreversible compaction of the material. Farther from the impact location, the confinement stresses decrease and the material experiences compression with a moderate triaxial state of stress. Finally, compressive wave reflection may occur and can result in a tensile wave which will interact with compressive waves and produce spalling, i.e. tensile cracking induced by wave interaction. Computational analysis of concrete and reinforced concrete elements subjected to this type of loading history requires the implementation of a constitutive relation capable of capturing the major features of the material response under such loads: tensile cracking, compression failure and the effect of confinement on the ultimate stress, and finally compaction which induces an increase of the tangent and elastic stiffness of the material, an increase of the wave speed and consequently substantial modifications of wave interactions which may affect tensile cracking and failure. Fig. 1 shows experimental results from [6] on which the hydrostatic stress has been plotted versus the volumetric strain for two different loading histories: triaxial hydrostatic compression and oedometric compression inspired from the testing procedure devised by Bazant et al. [3] (i.e. with zero, or quasi-vanishing, transverse strains). Looking at the material response and from a phenomenological point of view, constitutive models should incorporate irreversible plastic strains, an increase of stiffness due to compaction of concrete, and material damage due to progressive cracking in tension which does not occur in the experimental results shown in Fig. 1. Furthermore, these experimental results show that it is not possible to separate the deviatoric and hydrostatic responses of the material, as opposed to the assumption commonly used in most material
<figure><img src="image_2.png" /><figcaption>Fig. 1. Experimental hydrostatic and oedometric test results (after Burlion, 1997).</figcaption></figure>

models for impact problems. Would the hydrostatic and deviatoric responses be uncoupled, the two curves in Fig. 1 would be exactly the same.

Except maybe for the microplane model which is derived from different modelling principles [1,2]), phenomenological damage or plasticity models are usually capable to describe correctly the behaviour of concrete if one of the principle strains is positive and more generally in biaxial situations (see e.g. [21,10]). Most of these models are not capable, however, to capture the effects of severe compaction because of two difficulties:

1. Plastic strain growth under an hydrostatic state of stress and material stiffening (loading and unloading stiffness). This is due to the plastic/damage loading functions which should be closed on the hydrostatic axis, e.g. with a cap (see e.g. [20]).

2. Theoretical difficulties which prevent a decrease of damage, and thus an increase of elastic stiffness. From the viewpoint of thermodynamics of irreversible processes, damage cannot decrease in a closed system, unless another dissipative mechanism occurs at the same time which results in a positive global dissipation rate (which would be negative otherwise).

In this contribution, we present a model which has been developed, keeping in mind these two difficulties. The originality of the model compared to the approaches based on cap models in plasticity is the description of the variations of the material elastic stiffness. The degradation of the elastic moduli is described in the model by two damage scalars: a tension damage variable and a compression damage variable. Tension damage is controlled by the positive elastic strains. In order to capture the evolution of plastic strain we use a modified Gurson's yield function with associated flow rules [15,23]. The evolution of the volume fraction of voids entering in the Gurson's yield function is directly related to compression damage. When it decreases, it produces an increase of the material stiffness and therefore a decrease of material damage. Hence, the evolution of compression damage and plastic strain are entirely coupled. Section 2 presents this constitutive model and its calibration from experimental results. In Section 3, the implementation of this model in the finite element code DYNA3D is proposed. Finally, computations on test cases which exhibit the major properties and characteristics of the model are presented in Section 4. We do
not intend to deal here with thorough comparisons with experimental results, but rather to exhibit the influence of its major characteristics on the response of small structures in transient dynamics.

## 2. Constitutive model

In the following, we will use the small strain assumption and the standard decomposition of strain increments into elastic and plastic parts:

$$d\epsilon_{ij}(x, t) = d\epsilon_{ij}^{e}(x, t) + d\epsilon_{ij}^{p}(x, t), \qquad (1)$$

where $\epsilon_{ij}(x, t)$ are the total strain components, $\epsilon_{ij}^{e}(x, t)$ are the elastic strains, and $\epsilon_{ij}^{p}(x, t)$ are the irreversible strain components. We are going to present first the elastic-damage part of the model. Plasticity will be added later and coupled to the foregoing damage model.

### 2.1. Damage model

We focus in this section on the relationship between the elastic strain tensor and the stress tensor. There are two mechanisms which induce a variation of the elastic moduli of the material:

* The first one is microcracking. In tension, where we expect that the plastic strains are very small compared to those in compression, the material nonlinear response is due almost solely to microcrack growth. As we will see in Section 2.2, tension damage will be considered to be controlled by positive strains which produce microcrack opening and microcrack growth.

* In compression, and more specifically for states of stress where the confinement prevents microcrack opening, the other damage mechanism is the crushing of the cement or mortar matrix in concrete. It induces a variation of the volume fraction of voids which will later be considered to be driven by the irreversible plastic strains.

An important issue is now the type of damage variables to be used is whether the directionality of damage (understood here as microcracking and variation of void volume fraction) is important or not. Damage due to the variation of the volume fraction of voids, later on denoted as compression damage, can easily be considered to be an isotropic phenomenon. From a physical point of view, it is very clear, however, that microcracking in tension is geometrically oriented. Microcracks develop in planes which are perpendicular to the applied stress. One approach would be then to consider that tension damage growth induces in this case elastic orthotropy. Such models have been devised in the literature, among others, by Sidoroff [26], Chaboche et al. [7], Dragon and Mroz [9], Berthaud et al. [4], Valanis [28]. Recent results by Fichant et al. [11], however, have shown that in situations where material failure is essentially controlled by one-dimensional extensions, scalar damage models provide numerical predictions at the structural level which are very similar to orthotropic damage models. In fact, damage induced anisotropy was found to be important when the material is subjected to multiaxial extensions which induce damage growth, and also when the strain history applied to the material is severely non-radial.

In view of the complexity of anisotropic damage models compared to the isotropic (scalar) damage model, both from the viewpoint of calibration and of numerical implementation, we will consider that tension damage $d$ and compression damage $\delta$ are scalars, keeping in mind the limits of the simplifying assumption of isotropy and, in some instances, its merits.

We assume now that both damage mechanisms act on the elastic moduli independently. It follows that the global variation of damage $\dot{D}$, which measures the overall variation of the elastic moduli of the material, can be written as:

$$\dot{D} = \dot{d} + \dot{\delta}. \qquad (2)$$

It should be underlined here that compression and tension damages are rooted in different micromechanical mechanisms. It is their consequences on the elastic response of the material which is folded into a single variable $D$, which is not an internal variable in the thermodynamic sense [19]. Tensile damage and compression damage are independent internal state variables which refer each to a specific degradation
mechanism of the material at the micromechanical level. Their evolution is such that $D \in [0, 1]$ in order to have an overall stiffness which varies between two bounds: the stiffness of an uncracked, unvoided material and a vanishing stiffness.

This type of superposition of two damage variables was introduced in a similar form by Mazars [21], with a slightly different meaning since it was a weighted summation of damage due to uniaxial tension and damage due to uniaxial compression where the weights were non-dimensional functions of the state of strain.

The stress–strain relations are:

$$ \sigma_{ij} = (1 - D)[\lambda\epsilon_{kk}^e \delta_{ij} + 2\mu\epsilon_{ij}^e], \qquad (3) $$

where $\sigma_{ij}$ are the stress components, and $(\lambda, \mu)$ are the Lamé coefficients. The underlying assumption in Eq. (3) is that the bulk and shear (unloading) moduli of the material degrade at the same rate. Other combinations of tensile and compressive damage could be envisioned (see e.g. [17]).

Another issue which is encountered in the description of the response of concrete in tension is strain softening and strain and damage localisation (see e.g. de Borst et al. [5]). We do not intend to deal with this problem in this contribution, but the gradient approach to damage developed by Frémond and Nedjar [12], Peerlings et al. [24], could be implemented in order to circumvent this problem. Since we have two damage variables in the model, the gradient approach should be applied to (i) tensile damage (see e.g. [25]); and (ii) compressive damage. Leblond et al. [18] have developed a non-local approach for ductile damage which can be transformed into a gradient model and extended to the compression damage model presented in the following sections.

## 2.2. Tensile damage growth

The loading function defined in Eq. (4) and evolution equations (5) are used:

$$ f(\tilde{\epsilon}, \kappa) = \tilde{\epsilon} - \kappa, \qquad (4) $$

if $f(\tilde{\epsilon}, \kappa) = 0$ and $\frac{\partial f(\tilde{\epsilon}, \kappa)}{\partial \epsilon_{ij}^e} \epsilon_{ij}^e > 0$ then $\begin{cases} d = g_1(\kappa), \\ d\tilde{\epsilon} = d\kappa, \end{cases}$

else $d\kappa = 0$ \qquad (5)

with

$$ g_1(\kappa) = 1 - \frac{\kappa_0(1 - A_t)}{\kappa} - \frac{A_t}{\exp[B_t(\kappa - \kappa_0)]}, $$

where $\kappa$ is the damage history variable, $A_t$ and $B_t$ are experimentally determined model parameters, and $\tilde{\epsilon}$ is the effective tensile strain defined by Mazars [21]:

$$ \tilde{\epsilon} = \sqrt{\sum_{i=1}^{3} (\langle \epsilon_i^e \rangle_{+})^2}, \qquad (6) $$

where $\langle \epsilon_i^e \rangle_{+}$ is the positive part of the principal reversible strain $\epsilon_i^e$. These equations were introduced by Mazars [21] and can be also found for instance in [22]. The damage loading function defines the domain of reversible behaviour. When $f(\tilde{\epsilon}, \kappa) < 0$ damage does not grow. Eq. (5) is such that the damage history variable $\kappa$ is, at each material point, the maximum between the largest value of the effective tensile strain encountered during the loading history and the threshold value $\kappa_0$. Initially, before any load is applied on the virgin material, $\kappa = \kappa_0$.

## 2.3. Compression damage and plasticity

In order to capture material stiffening observed in Fig. 1, which is due to a collapse of the porous structure of concrete (mortar), compression damage is coupled to plasticity. More precisely, the incremental
evolution of the compression damage variable will be related to the first invariant of the irreversible plastic strain.

A modified version of the Gurson's plastic yield function is chosen [15,23]:

$$F_{NT}(\sigma_{ij}, \sigma_M, f^*) = \frac{3J_2}{\sigma_M^2} + 2q_1 f^* \cosh\left(q_2 \frac{I_1}{2\sigma_M}\right) - (1 + (q_3 f^*)^2) = 0, \quad (8)$$

where $I_1$ is the first stress tensor invariant, $I_1 = \sigma_{kk}$, $J_2$ the second invariant of the deviatoric stress tensor defined by:

$$J_2 = \frac{1}{2}(s_{ij}s_{ij}) \qquad s_{ij} = \sigma_{ij} - \frac{1}{3}\sigma_{kk}\delta_{ij}, \qquad (9)$$

where $\sigma_M$ is the equivalent yield stress in the matrix, $q_1$, $q_2$, and $q_3$ are model parameters. In Gurson's yield function $f^*$ is the material porosity (volume fraction of voids). $f^*$ increases with void development in tension, and decreases with void closure in triaxial compression. The material porosity will be later on related to compression damage. Note also that the deviatoric material response is coupled to the hydrostatic response, in agreement with the observations in Fig. 1. The evolution of modified Gurson's yield surface with the decrease of porosity is presented in Fig. 2. Due to symmetry, the yield surface is only plotted on a quarter of the $(I_1/\sigma_M, 3J_2/\sigma_M)$ space. The yield surface grows due to the closure of the material porosity (or volume fraction of voids) ($f^* \to 0$).

It is important to point out at this stage that the Gurson's yield surface is not expected to control the material response in tension. The yield limits in tension and compression are the same according to the Gurson's yield surface. Plastic strains will not develop in tension because the level of (positive) stress required in order to reach plastic yielding is way above the limit of elastic reversible behaviour defined by the damage loading function. Hence, in tensile loadings, the material will experience damage growth due to microcracking and not plastic strain growth. Finally, it must be observed that when the volume fraction of voids is equal to zero the loading function reduces to the classical Von Mises expression. It is however quite difficult to decide whether the Von Mises expression is adequate or not because a vanishing volume fraction of void corresponds to the limit case of an entirely compacted material which, in fact, would suffer severe microstructural reorganisation due to the application of a very high hydrostatic compression (several GPa). Furthermore, such states of stress induce phase changes both in mortar and in the aggregates which are not

<figure><img src="image_3.png" /><figcaption>Fig. 2. Modified Gurson yield function: evolution with the decrease of porosity ($q_1 = q_2 = q_3 = 1$)</figcaption></figure>
accounted for in the present model whose range of application in the hydrostatic compression regime is above −1 GPa.

The normality rule gives the standard expressions for the irreversible strains:

$$d e_{ij}^{p} = d \lambda \frac{\partial F_{NT}}{\partial \sigma_{ij}} \qquad (10)$$

with the Kuhn–Tucker relations

$$F_{NT} \le 0, d\lambda \ge 0, \text{ and } F_{NT} d\lambda = 0, \qquad (11)$$

where $d\lambda$ is the plastic multiplier. The work hardening equation for a material with voids was given by Needleman and Tvergaard [23]:

$$\sigma_{ij} d e_{ij}^{p} = (1 - f^{*}) \sigma_{M} d e_{M}^{p}, \qquad (12)$$

where $e_M^p$ is the irreversible equivalent strain in the matrix (i.e. the material without voids) associated to the matrix equivalent stress $\sigma_M$. The relation between the matrix equivalent strain and the matrix equivalent stress is defined by

$$\epsilon_M = \begin{cases} \frac{\sigma_M}{E} & \text{if } \sigma_M \le \sigma_y \\ \frac{\sigma_y}{E} \left( \frac{\sigma_M}{\sigma_y} \right)^n & \text{if } \sigma_M > \sigma_y, \end{cases} \qquad (13)$$

where $\sigma_y$ is the elastic strength of the matrix, and $n$ the hardening exponent.

The decrease of the void volume fraction $f^*$ is controlled by the plastic flow. Same as in the Gurson's model, we assume that microvoid evolution is controlled by the irreversible volumetric strain [6]:

$$d f^* = k(1 - f^*) f^* d e_{kk}^p, \qquad (14)$$

where $k$ is a model parameter which controls the rate at which the porosity changes as a consequence of the plastic flow. For standard concrete, the initial value of $f^*$ is roughly equal to 0.25, corresponding to the expected initial porosity of the material.

We need now to relate the variation of volume fraction of voids to the evolution of the elastic moduli. From micromechanics, one can use existing relations available in the literature. An example was provided by Colantonio and Stainier [8] who used the same plastic yield function and Mori–Tanaka theory for a two phase material in order to relate the variation of porosity to the variation of the elastic constants. The difficulty lies, however, in the combination of compression and tension damages which would be more complex than the one assumed in Eq. (2). In a first approximation we will set:

$$f^* = \delta. \qquad (15)$$

According to this equation, the variation of porosity is coupled to the variation of the Young's modulus of the material and yield elastic stiffening/softening effects which are not present in standard Gurson's type models.

## 2.4. Model response and calibration

The determination of the model parameters benefits from the fact that in uniaxial tension damage due to microcracking occurs without plastic strain growth and conversely, in compression dominated regime (with a large enough confinement in order to avoid positive reversible strains) plastic strain growth is observed with the inherent variation of compression damage. In order to calibrate the model, one needs three types of tests: a tensile test is required in order to calibrate damage growth due to tensile strain. Two compression tests are required in order to calibrate the compression damage–plasticity model. The reason is that in order to calibrate the evolution of the porosity with the plastic strain, one needs an experiment where the compression hydrostatic stress is large enough. An oedometric test and a hydrostatic compression test can be used for the model calibration (see [6]). Note that the calibration of the model in the compressive regime
where microcracking does not occur bears the same difficulties as the calibration of usual Gurson's models for ductile damage, except for the compression damage growth, which requires loading/unloading cycles in the experiments.

Fig. 3 shows the fits with the experimental data obtained by Burlion [6] on hydrostatic compression experiments (a) and uniaxial confined compression test (b). The values of model parameters obtained for a mortar are given in Table 1. The plot in Fig. 3(b) shows the axial stress vs. axial strain response of the material obtained experimentally. Further details on the experiments can be found in [6]. Note that this set of two experimental data is exactly the same (in a slightly different coordinate system for Fig. 3(b) as the one in Fig. 1). Hence the model is capable of reproducing the coupling between the hydrostatic and deviatoric responses underlined in the introduction.

<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Comparison between the experimental results and the model simulation for mortar under hydrostatic compression
(b) Comparison between the experimental results and the model simulation for mortar under oedometric compression.</figcaption></figure>
<table><caption>Table 1
Model parameters for the comparisons and simulations in Figs. 3–6</caption><thead><tr><td>Parameters</td><td>Units</td><td>Value</td></tr></thead><tbody><tr><td>q<sub>1</sub></td><td>–</td><td>0.5</td></tr><tr><td>q<sub>2</sub></td><td>–</td><td>0.7</td></tr><tr><td>q<sub>3</sub></td><td>–</td><td>0.8</td></tr><tr><td>σ<sub>y</sub></td><td>MPa</td><td>25</td></tr><tr><td>n</td><td>–</td><td>10</td></tr><tr><td>k</td><td>–</td><td>1</td></tr><tr><td>d<sub>0</sub></td><td>–</td><td>0.3</td></tr><tr><td>f<sub>0</sub><sup>*</sup></td><td>–</td><td>0.3</td></tr><tr><td>E</td><td>MPa</td><td>26 000</td></tr><tr><td>ν</td><td>–</td><td>0.2</td></tr><tr><td>A<sub>t</sub></td><td>–</td><td>0.8</td></tr><tr><td>B<sub>t</sub></td><td>–</td><td>20 000</td></tr><tr><td>κ</td><td>–</td><td>1 · 10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Numerical simulation: stress–strain evolution in uniaxial and hydrostatic tension.</figcaption></figure>

Fig. 4 shows numerical simulations in uniaxial tension and in hydrostatic tension. The model parameters are similar to those used by Mazars [21] since the present constitutive provides a similar response. It is very clear from loading-unloading cycles that plastic strain do not grow and that the non linear response of the material is entirely controlled by damage due to microcracking. Fig. 5 shows the variation of the Young's modulus of the material versus the total strain in hydrostatic compression followed by hydrostatic tension. We can observe that the main characteristics involved in the compaction of the material are captured as expected. The Young's modulus increases due to the triaxial compression, and decreases when the material is subjected to hydrostatic tension. Fig. 6 shows the stress–strain relation obtained with the model in uniaxial compression. We can see that the model response in uniaxial compression is not very well described by the model. This particular point should be the objective of future developments. Nevertheless, the global form of a classical uniaxial compression curve is acceptable in the computations as the main characteristics are preserved (compressive strength, strain at peak stress, softening regime, ...).
<figure><img src="image_6.png" /><figcaption>Fig. 5. Numerical simulation: modulus evolution in hydrostatic compression and hydrostatic tension.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Numerical simulation: stress–strain evolution in uniaxial compression.</figcaption></figure>

3. Time integration of the constitutive model

THE main objective of the constitutive relations presented in the previous section is the description of the response of concrete in transient dynamics. The constitutive relations have been implemented in the commercial finite element code DYNA3D. It is a vectorised explicit three-dimensional finite element code for analysing the large deformation dynamic response of inelastic solids. The equations of motion are integrated with the explicit central difference method.
At this point, it is necessary to consider the complexity of the constitutive relations developed in Section 2. The model has two loading functions defined in two different spaces. The damage loading function is defined in the one-dimensional space of the effective tensile strain and the plasticity loading function is defined in the stress space. Because the stresses depend on the elastic strains and elastic stiffness, the two loading functions are entirely coupled. Plastic damage models are not that complex to implement usually. In most cases [16,27,11], it is possible to separate the effects of damage and plasticity because the plastic yield function is expressed as a function of the effective stress which is independent of damage. Standard algorithms can be implemented for the integration of the plastic part of the model and the damage part is computed afterwards. Generally, damage is an explicit function of the elastic strain as defined in Eqs. (4)–(6).

We have chosen here to implement an explicit, Euler forward integration scheme because time integration needs to be extremely fast in 3D transient dynamics explicit codes where time increments are very small and where the number of degree of freedom is very high. Clearly, this integration scheme is less accurate than an implicit one since the tensile damage and yield conditions are not satisfied at the end of each time step. Implicit return mapping algorithms would be more appropriate in quasistatic computations but not very efficient as far as computing time is concerned in the present context.

We are going in a first step to derive the incremental stress–strain relations under the assumption that damage and plastic strain grow. The consistency condition on the plastic yield surface yields:

$$dF_{\text{NT}}(\sigma_{ij}, \sigma_M, f^*) = \frac{\partial F_{\text{NT}}}{\partial \sigma_{ij}} d\sigma_{ij} + \frac{\partial F_{\text{NT}}}{\partial \sigma_M} d\sigma_M + \frac{\partial F_{\text{NT}}}{\partial f^*} df^* = 0, \quad (16)$$

where

$$\begin{align}
\frac{\partial F_{\text{NT}}}{\partial \sigma_{ij}} &= \frac{3}{\sigma_M^2} s_{ij} + \frac{q_1 q_2 f^*}{\sigma_M} \left[ \sinh \left( \frac{q_2 I_1}{2 \sigma_M} \right) \right] \delta_{ij}, \nonumber \\
\frac{\partial F_{\text{NT}}}{\partial \sigma_M} &= -2 \frac{J_2}{\sigma_M^3} - \frac{q_1 q_2 f^* I_1}{\sigma_M^2} \left[ \sinh \left( \frac{q_2 I_1}{2 \sigma_M} \right) \right], \tag{17} \\
\frac{\partial F_{\text{NT}}}{\partial f^*} &= 2q_1 \left[ \cosh \left( \frac{q_2 I_1}{2 \sigma_M} \right) \right] - 2q_3 f^*. \nonumber
\end{align}$$

The porosity evolution given by Eq. (14) is rewritten with Eq. (10) as:

$$d f^* = d \lambda k (1 - f^*) f^* \frac{\partial F_{\text{NT}}}{\partial \sigma_{ij}} \delta_{ij}. \quad (18)$$

The relation between the equivalent matrix strain and the equivalent matrix stress is deduced from Eq. (13):

$$d\sigma_M = E_t d\varepsilon_M \quad \text{with} \quad \frac{1}{E_t} = \frac{\partial[\sigma_y/E(\sigma_M/\sigma_y)^n]}{\partial\sigma_M} = \frac{n}{E} \left(\frac{\sigma_M}{\sigma_y}\right)^{n-1}. \quad (19)$$

Substitution of this equation in Eq. (12) provides an expression of the increment of equivalent matrix stress as a function of the hardening modulus $E^*$, the initial porosity, and the plastic multiplier:

$$d\sigma_M = d\lambda \frac{\sigma_{ij}(\partial F_{\text{NT}}/\partial\sigma_{ij})E^*}{(1-f^*)\sigma_M} \quad \text{with} \quad \frac{1}{E^*} = \left(\frac{1}{E_t} - \frac{1}{E}\right). \quad (20)$$

Substitution of the previous equations (18) and (20) in the consistency condition (Eq. (16)) provides an equation from which the plastic multiplier is the single unknown:

$$\left[ \frac{3}{\sigma_M^2} s_{ij} + \frac{q_1 q_2 f^*}{\sigma_M} \left[ \sinh \left( \frac{q_2 I_1}{2 \sigma_M} \right) \right] \delta_{ij} \right] d\sigma_{ij} - \left[ 2 \frac{J_2}{\sigma_M^3} + \frac{q_1 q_2 f^* I_1}{\sigma_M^2} \left[ \sinh \left( \frac{q_2 I_1}{2 \sigma_M} \right) \right] \right] d\lambda \frac{\sigma_{ij}(\partial F_{\text{NT}}/\partial \sigma_{ij})E^*}{(1-f^*)\sigma_M} \\ + \left[ 2q_1 \left[ \cosh \left( \frac{q_2 I_1}{2 \sigma_M} \right) \right] - 2q_3 f^* \right] d\lambda k (1-f^*) f^* \frac{\partial F_{\text{NT}}}{\partial \sigma_{ij}} \delta_{ij} = 0. \quad (21)$$
This equation can be recast in the following abbreviated form:

$$d\lambda = fct(d\sigma_{ij}) \qquad (22)$$

We may now consider the incremental stress strain relation where the plastic multiplier is assumed to be determined from Eq. (22):

$$d\sigma_{ij} = (1-d)C_{ijkl}\left(d\varepsilon_{ij} - d\lambda\frac{\partial F_{NT}}{\partial\sigma_{ij}}\right) - ddC_{ijkl}\varepsilon_{kl}^{e}, \qquad (23)$$

where $C_{ijkl}$ is the elastic tensor. Assuming that tensile damage, plastic strains and compression damage are evolving at the same time, the incremental stress-strain relation becomes:

$$d\sigma_{ij} = (1 - D)C_{ijkl}\left(d\varepsilon_{ij} - d\lambda\frac{\partial F_{NT}}{\partial\sigma_{ij}}\right) - \left[\frac{\partial g_1(\tilde{e})}{\partial\tilde{e}}\frac{\partial\tilde{e}}{\partial\varepsilon_{ij}^e}\left(d\varepsilon_{ij} - d\lambda\frac{\partial F_{NT}}{\partial\sigma_{ij}}\right) + d\lambda\frac{\partial g_2(f^*)}{\partial f^*}k(1 - f^*)f^*\frac{\partial F_{NT}}{\partial\sigma_{ij}}\delta_{ij}\right] \times C_{ijkl}\varepsilon_{kl}^e. \qquad (24)$$

Substitution of Eq. (22) in this expression provides a 6 × 6 system of equations from which the explicit final form of the stress increment (6 components) are computed as functions of the total strain increment and quantities related to plasticity and damage which are known at the beginning of the increment.

Within a predictor-corrector scheme, an elastic predictor is computed first. The effective tensile strain is computed and the damage loading function is tested.

* If the effective strain Eq. (4) lies within the damage loading surface, we assume that damage due to tensile strains is equal to zero. Then, the plasticity yield function is tested according to the predictor and an explicit plastic correction is computed if it is necessary. This correction is obtained from Eqs. (22) and (23) assuming that damage is constant.

* If the effective strain lies outside the damage loading function, damage will be considered to evolve in Eq. (23). The plasticity yield function is tested again and if the trial stress lies outside the yield function, Eq. (24) is used to compute the new incremental stress.

Damage is always tested first because if there is no plasticity, the tensile strains will be in most cases the largest possible. Nevertheless, the damage loading function is checked when the elastic strains have been determined according to the plasticity model. Corrections are made if damage should evolve, the full equation (24) is used.

Because it is explicit, the time integration of the constitutive relation may fail if the time increments in the finite element calculation are too large. Figs. 7 and 8 show two examples on which the influence of the strain increment size has been tested. The corresponding model parameters are indicated in Table 2. In Fig. 7, the loading history is such that in the principal direction $\varepsilon_1 : \varepsilon_2 : \varepsilon_3 = -1 : 0.2 : 0.2$. In the elastic regime, it corresponds to uniaxial compression. Outside the elastic regime, a transverse stress is induced by the transverse plastic strains. Fig. 7(a) shows the longitudinal stress $\sigma_1$ as a function of the longitudinal strain $\varepsilon_1$ and Fig. 7(b) shows the transverse stress $\sigma_2$ as a function of the longitudinal strain. The second test case is hydrostatic compression where $\varepsilon_1 : \varepsilon_2 : \varepsilon_3 = -1 : -1 : -1$ (Fig. 8). In both situations, the time integration converges as the strain increments size is decreased and the error and its accumulation does not seem to be extremely important. The largest error is observed on Fig. 7(b) where the maximum transverse strain is quite overestimated for large strain increments. When yielding is initiated, there is no search for the contact stress on the yield surface. This is most probably, with the damage-plasticity interaction which changes, the reason for the errors observed on Fig. 7. It is interesting to point out that in structural computations with DYNA3D, the time increment are generally very small due to the element size which should be small in regions where material nonlinearities occur. Consequently, the strain increments remain relatively small, typically lower than the smallest ones on Figs. 7 and 8. Therefore, one can consider that the time integration is accurate enough, or at least is a compromise between accuracy and computational speed. This scheme can be enhanced by subdividing the strain increment if it is too large. An automatic procedure for this purpose would need to be developed.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. (a) Longitudinal stress vs. longitudinal strain and (b) Transverse stress vs. longitudinal strain for several strain increments.</figcaption></figure>

## 4. Computational examples

Two types of computations are going to be presented in the following. The first one is a simulation of the split Hopkinson test aimed at demonstrating the influence of the variation of porosity and inherent stiffening on wave propagation. The second test case is an impact on a concrete rod.

In the split Hopkinson test, the input bar is impacted with a striker at an initial velocity which is an experimental parameter. In the input bar, a stress wave is developed. This wave arrives on the specimen and becomes the specimen loading. In the output bar, new waves are developed (Fig. 9). The experimental data are obtained by measuring the strains in the input and output bars. With these strains, we can obtain the experimental velocities and forces (along the cylinder axis) applied on the two faces of the concrete specimen in contact with the bars. In order to show the influence of confinement, a special specimen has been
<figure><img src="image_10.png" /><figcaption>Fig. 8. Pressure vs. volumetric strain for different strain increments.</figcaption></figure>

<table><caption>Table 2
Model parameters for the finite element computations</caption><thead><tr><td>Parameter</td><td>Unit</td><td>Value</td></tr></thead><tbody><tr><td>q<sub>1</sub></td><td>–</td><td>1.5</td></tr><tr><td>q<sub>2</sub></td><td>–</td><td>0.9</td></tr><tr><td>q<sub>3</sub></td><td>–</td><td>1.5</td></tr><tr><td>σ<sub>y</sub></td><td>MPa</td><td>70</td></tr><tr><td>n</td><td>–</td><td>10</td></tr><tr><td>d<sub>0</sub></td><td>–</td><td>0.3</td></tr><tr><td>f<sub>0</sub><sup>*</sup></td><td>–</td><td>0.3</td></tr><tr><td>E</td><td>MPa</td><td>40 000</td></tr><tr><td>ν</td><td>–</td><td>0.2</td></tr><tr><td>k</td><td>–</td><td>105</td></tr><tr><td>A<sub>t</sub></td><td>–</td><td>0.8</td></tr><tr><td>B<sub>t</sub></td><td>–</td><td>20 000</td></tr><tr><td>κ</td><td>–</td><td>1 · 10<sup>-4</sup></td></tr><tr><td>E steel</td><td>MPa</td><td>200 000</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Split Hopkinson test.</figcaption></figure>
designed, made of a concrete cylinder embedded in a metal jacket. The metal jacket controls the radial
deformation of the specimen and therefore applies a confinement stress on concrete which avoids splitting.
The friction between the concrete sample and the metal jacket has been neglected because, experimentally,
the contact surface is coated with Teflon.

Such experiments are currently being performed [14]. Numerical simulations of this split Hopkinson test with different evolutions of the porosity are presented here. Fig. 10 shows the various hydrostatic stress versus volumetric variation curves. These different curves where obtained by varying the model parameter *k* in Eq. (14). The other parameters are constant and correspond to those on Table 2. In these simulations the input and output bars are not represented in the finite element model. Only the concrete specimen is described with special boundary conditions which account for wave transmission and reflection. Axial velocities at the boundaries are prescribed by mean of rigid surfaces. These boundary conditions have in fact little influence on the computed forces [13]. The external radial displacements of the specimen (steel jacket) are free. The input and output experimental velocities are prescribed to the confined concrete specimen and the input and output forces are computed. Fig. 11 shows that the computed forces increase as the material becomes more and more prone to compaction. This is expected because for a given volumetric strain, the hydrostatic stress increases with the compaction (Fig. 10). One important consequence of this result is that in impact problems, because of the confinement near the impactor, the material should be able to carry out higher stresses. The impactor speed should be very much decreased because the material can oppose higher stresses.

The second example is an impact problem on a rod shown on Fig. 12. The concrete rod is 20 cm long, with a square cross section of 1 cm×1 cm. It is discretised with 500 constant size elements (of cross section 1 cm×1 cm). Radial displacements are blocked and longitudinal displacements are allowed only. The steel striker of length 0.5 cm impacts the rod with a speed of 100 m/s. The model parameters are those in Table 2. Fig. 13 shows the strain wave propagation in the case where the rod is elastic. Times $t_1-t_5$ were chosen in order to observe distinctly the stress wave shape when it propagates in the rod. The wave reaches the free boundary of the rod just before the instant $t_4$ and becomes a tensile wave as it is reflected. Fig. 14(a) and (b) show results with the plastic-damage model. Times $t'_1-t'_5$ and $t_1-t_6$ are different and were chosen in order to observe easily the stress wave evolution during propagation.

<figure><img src="image_12.png" /><figcaption>Fig. 10. (a) Pressure vs. volumetric strain for different values of K.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Input (a) and output (b) forces for different porosity evolutions.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Finite Element Mesh.</figcaption></figure>

Compaction does not occur in the computation shown in Fig. 14(a) ($k = 0$), i.e. the volume fraction of voids remains constant. The amplitude of the compressive wave is decreased as it propagates. At time $t_4$ plastic strains have disappeared. The plateau observed corresponds to the yield limit. The signal duration increases too because the elastic wave propagates faster than the plastic one. The plastic wave corresponding
<figure><img src="image_16.png" /><figcaption>Fig. 13. Stress waves propagation for an elastic model.</figcaption></figure>

to the peak on the signal (largest negative stress) vanishes because of the unloading wave which is part of the signal and propagates faster than the plastic one and also because energy is consumed during the plastic strain growth. Upon reflection of the wave on the free boundary, a tensile wave is generated and will produce spalling of the bar if the model incorporates tensile damage. Same as in classical modelling of spalling, the size of the spall is controlled by the wave interaction.

Fig. 14(b) shows the same computation where compaction is included. Note that, in this example, the state of strain is similar to that of the oedometric test shown in Fig. 1 where compaction effects are already very important. It is important to remark that the amplitude of the stress is much lower. Because of the compaction, the unloading wave propagates faster and cancels the plastic wave sooner than in the situation where compaction is omitted. The plastic zone in the bar is smaller and it is expected that the structure will be less damaged due to the impact if compaction is described.

**5. Conclusions**

A damage plasticity model has been presented. The constitutive relations possess several characteristics:

* Two different mechanisms of degradation of the elastic moduli are considered: tension damage due to microcracking and compression damage due to microvoid growth or collapse. These two types of material damage are described by scalars and their influence on the elastic stiffness is superimposed, with the assumption that the damaged material remains isotropic.

* Tensile damage growth is controlled by positive elastic strains. A modified Gurson yield function is used in order to capture the decrease of the porosity in concrete due to plastic strains observed on experimental results. The porosity of the material is directly related to compression damage. The model describes compaction phenomena quite well. Nevertheless, some improvements are expected for uniaxial compression.

Identification of the parameters of the model requires 3 types of experiments: the first one is a tension test needed to calibrate the model in tension; the second and the third ones are two different confined compression tests. They are required in order to find the parameters for plasticity model and thus the evolution of compression damage.

The model has been implemented in the explicit finite element code DYNA3D. An explicit time integration of the constitutive relation was used for the sake of simplicity and of computational efficiency in the
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Stress wave propagation for the compaction model: (a) without variation of the volume fraction of voids, (b) with variations of the volume fraction of voids.</figcaption></figure>

context on transient dynamics 3D computations. The examples demonstrate the importance of compaction on the energy dissipated during the impact. It is expected that compaction will greatly modify the response of concrete structures. In an impact problem it should induce a decrease of the size and depth of the plastic zone under the impactor. In perforation problems, it should contribute at lowering the output speed of the projectile.

## Acknowledgements

Financial supports from TDA Armements SAS under a project directed by Dr. Th. Bouet and from the French Research Network GEO are gratefully acknowledged (Program coordinator, Dr. P. Bailly).
References

[1] Z.P. Bazant, P.C. Prat, Microplane model for brittle plastic materials. I. Theory, J. Eng. Mech., ASCE 114 (10) (1988) 1672-1688.

[2] Z.P. Bazant, P.C. Prat, Miroplane model for brittle plastic materials. II. Verification, J. Eng. Mech., ASCE 114 (10) (1988) 1689-1702.

[3] Z.P. Bazant, F.C. Bishop, T.P. Chang, Confined compression tests of cement paste and concrete up to ksi, ACI J.33 (1986) 553-560.

[4] Y. Berthaud, C. Laborderie, S. Ramtani, Damage modelling and crack closure effects, in: J.W. Ju et al. (Eds.), Damage Mechanics in Engrg. Materials, ASME AMD vol. 109 1990, pp. 263-276.

[5] R. de Borst, L.J. Sluys, H.B. Muhlhaus, J. Pamin, Fundamental issues in finite element analyses of localization of deformation, Engrg. Comput. 10 (1993) 99-121.

[6] N. Burlion, Compaction des bêtons: éléments de modélisation et caractèrisation expérimentale, thèse de Doctorat de l'ENS de Cachan, Cachan, France, 1997.

[7] J.L. Chaboche, P.M. Lesne, J.F. Maire, Phenomenological Damage Mechanics of Brittle Materials with Description of the Unilateral Damage Effect, in: Z.P. Bazant et al. (Eds.), Fracture and Damage in QuasibrITTLE Structures, E&FN Spon Pubs, London, 1994, pp. 75-84.

[8] L. Colantonio, L. Stainier, Numerical Integration of Visoplastic Constitutive Equations for Porous Materials, in: J.-A. Desideri, P. Le Tallec, E. Onate, J. Periaux, E. Stein, (Eds.), Numerical Methods in Engineering'96, Wiley Pubs., 1996, pp. 28-34.

[9] A. Dragon, Z. Mroz, A continuum model for plastic-brittle behaviour of rock and concrete, Int. J. Engrg. Sci. 17 (1979) 121-137.

[10] P.H. Feenstra, Computational aspects of Biaxial Sress in Plain and Reinforced Concrete, Ph.D. Dissertation, Delft University of Technology, The Netherlands, 1993.

[11] S. Fichant, Ch. LaBorderie, G. Pijaudier-Cabot, Isotropic and Anisotropic Descriptions of Damage in Concrete Structures, Int. J. of Mechanics of Cohesive-Frictional Materials 5 (1998) 339-359.

[12] M. Frémond, B. Nedjar, Endommagement et principe des puissances virtuelles, Comptes Rendus Acad. Sci. série II (1993) 857-864.

[13] F. Gatuingt, L. Daudeville, Numerical simulation of split Hopkinson tests, First European LS-DYNA3D users conference, Stratford England, 1997.

[14] F. Gatuingt, G. Gary, N. Burlion, L. Daudeville, Comportement des ouvrages en dynamique rapide : campagne expérimentale et modélisation, Proc. Colloque GEO'97, Aussois, France, 1997.

[15] A.L. Gurson, Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media, Engrg. Materials and Technology 99 (1977) 2-15.

[16] J.W. Ju, On energy-based coupled elasto-plastic damage theories: constitutive modeling and computational aspects, Int. J. of Solids Structures 25 (7) (1989) 803-833.

[17] P. Ladeveze, Sur une théorie de l'endommagement anisotrope, Int. report N 34, Laboratoire de Mécanique et Technologie, Cachan, France; presented at the CNRS Int. Colloquium on failure criteria of structured media, Villars de Lans, June, France; published in English in the proceedings Failure Criteria of Structured Media, J.P. Boehler (Ed.), Balkema Pubs, 1993, 1983, pp. 355-364.

[18] J.B. Leblond, G. Perrin, J. Devaux, Bifurcation effects in dustile metals with nonlocal damage, J. of Applied Mech. ASME 61 (1994) 236-242.

[19] J. Lemaitre, A course on Damage Mechanics, Springer Pubs, Berlin, 1992.

[20] V.A. Lubarda, S. Mastilovic, J. Knap, Brittle-Ductile transition in porous rocks by cap model, J. of Engrg. Mech. ASCE 122 (7) (1996) 633-642.

[21] J. Mazars, Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure, Thèse de Doctorat ès Sciences, Université Paris 6, 1984.

[22] J. Mazars, G. Pijaudier-Cabot, Continuum damage theory- application to concrete, J. Engrg. Mech. ASCE 115 (1989) 345-365.

[23] A. Needleman, V. Tvergaard, An analysis of ductile rupture in notched bars, J. Mech. Phys. Solids 32 (1984) 461-490.

[24] R.H.J. Peerlings, R. De Borst, W.A. Brekelmans, J.H.P. De Vree, Gradient enhanced damage for quasi-brittle materials, Int. J. Num. Meth. Engrg. 39 (1996) 3391-3403.

[25] G. Pijaudier-Cabot, N. Burlion, Damage and localization in elastic materials with voids, Int. J. of Mechanics of Cohesive-Frictional Materials 1 (1996) 129-144.

[26] F. Sidoroff, Description of Anisotropic Damage Application to Elasticity, in: J. Hult, J. Lemaitre (Eds.), Proceedings IUTAM Symposium on Physical Nonlinearites in Structural Mechanics, Springer Pubs, Berlin, 1981, pp. 237-244.

[27] J.C. Simo, J.W. Ju, Stress and strain based continuum damage models-II. Computational aspects, Int. J. of Solids Structures 23 (1987) 841-869.

[28] K.C. Valanis, A global damage theory and the hyperbolicity of the wave problem, J. of Applied Mech. ASME 58 (1991) 311-316.