# An elastic plastic damage formulation for concrete

Ludovic Jason, Antonio Huerta, Gilles Pijaudier-Cabot, Shahrokh Ghavamian

► To cite this version:

Ludovic Jason, Antonio Huerta, Gilles Pijaudier-Cabot, Shahrokh Ghavamian. An elastic plastic damage formulation for concrete. Computer Methods in Applied Mechanics and Engineering, 2006, 195 (52), pp.7077-7092. 10.1016/j.cma.2005.04.017. hal-00351799

HAL Id: hal-00351799

https://hal.science/hal-00351799v1

Submitted on 3 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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 CC0 1.0 - Universal - International License
An elastic plastic damage formulation for concrete: Application
to elementary tests and comparison with an isotropic damage model

Ludovic Jason<sup>a,b</sup>, Antonio Huerta<sup>c</sup>, Gilles Pijaudier-Cabot<sup>a</sup>, Shahrokh Ghavamian<sup>b</sup>

<sup>a</sup> R&DO, GeM—Institut de Recherche en Génie Civil et Mécanique, Ecole Centrale de Nantes, Université de Nantes—CNRS,
1 rue de la Noé, BP 92101, F44300 Nantes Cedex 3, France

<sup>b</sup> EDF R&D, 1 avenue du Général de Gaulle, 92141 Clamart Cedex, France

<sup>c</sup> Laboratori de Càlcul Numèric, Departament de Matemàtica Aplicada III Universitat Politècnica de Catalunya, Jordi Girona 1-3 E-08034 Barcelona, Spain

Pure elastic damage models or pure elastic plastic constitutive laws are not totally satisfactory to describe the behaviour of concrete. They indeed fail to reproduce the unloading slopes during cyclic loading which define experimentally the value of the damage in the material. When coupled effects are considered, in particular in hydro-mechanical problems, the capability of numerical models to reproduce the unloading behaviour is essential, because an accurate value of the damage, which controls the material permeability, is needed. In the context of very large size calculations that are needed for 3D massive structures heavily reinforced and pre-stressed (such as containment vessels), constitutive relations ought also to be as simple as possible. Here an elastic plastic damage formulation is proposed to circumvent the disadvantages of pure plastic and pure damage approaches. It is based on an isotropic damage model combined with a hardening yield plastic surface in order to reach a compromise as far as simplicity is concerned. Three elementary tests are first considered for validation. A tension test, a cyclic compression test and triaxial tests illustrate the improvements achieved by the coupled law compared to a simple damage model (plastic strains, change of volumetric behaviour, decrease in the elastic slope under hydrostatic pressures). Finally, one structural application is also considered: a concrete column wrapped in a steel tube.

Keywords: Concrete; Damage; Plasticity

# 1. Introduction

The evaluation of the integrity of a concrete structure first supposes the knowledge of its mechanical properties. When continuum mechanics is considered, elastic damage models or elastic plastic constitutive laws are generally the standard approach to describe the behaviour of concrete. In the first case, the mechanical effect of the progressive microcracking and strain softening are represented by a set of internal state variables which act on the elastic behaviour (decrease of the stiffness) at the macroscopic level (see e.g. [42,43,29]). In plasticity models, softening is directly included in the expression of a plastic yield surface by means of a hardening–softening function generally [14,5,30].

For some concrete structures in sensitive environments, such as confinement buildings of nuclear power plants, durability is not only related to mechanical load bearing capacities but also, for instance, to some transfer properties. In containment vessels which is the background of this contribution, the leakage rate upon a growth of internal pressure is a key safety factor. Therefore, the interaction between the mechanical properties and the permeability of concrete must be
<figure><img src="image_2.png" /><figcaption>Fig. 1. Experimental relation between damage and permeability. $k_v$ and $k_e$ are the intrinsic permeabilities of the initial and damaged material, respectively, [36]. OC = ordinary concrete, HPC = high performance concrete, HPFC = high performance fiber reinforced concrete.</figcaption></figure>

captured as accurately as possible. Recent experimental results [36] have shown that the permeability of concrete to gas was mainly controlled by the degradation of the unloading stiffness of the material (Fig. 1). This observation was made in the regime of diffuse damage, prior to localisation and failure. In the localised cracking regime, it is not the amount of damage that is the controlling parameter on the leakage rate, but rather the cracks opening and roughness. For the industrial applications aimed at in this paper, it is the diffuse microcracking regime which is of practical interest since nuclear containment vessels are not expected to contain major through cracks for obvious safety reasons. Design and commissioning standards restrict the stresses applied on the material to the pre-peak regime, with in fact as little as possible incursions in the non-linear regime. A theoretical background to the correlation between damage and permeability in the diffuse regime has been proposed by Chatzigeorgiou et al. [4]. It is shown that the variation of the permeability of a damaged sample (unloaded) is related to the variation of its elastic stiffness.

In this type of problem, it is therefore very important to capture the variations of the elastic (unloading) stiffness of the material upon mechanical loading. Continuum damage is the right theoretical framework for that. But continuum damage models cannot, alone, be implemented. The material undergoes also some irreversible deformations during loading. As sketched in Fig. 2, even if damage or plasticity models taken separately are capable of capturing the same material response upon monotonic loading, both theories do not capture the evolution of unloading stiffness accurately. At a given point on the stress-strain response, neglecting plastic strains according to a pure damage approach would result in an artificial increase of damage (secant unloading slope). Neglecting damage effects in a pure plastic model would result in a permeability that does not evolve, or is not described objectively when it is a function of the applied strains or stresses [4]. By objective, we mean that for the same value of strain (or stress) invariant combinations, different values of permeability are to be expected. Hence, coupled damage-plasticity models are a requisite in hydro-mechanical problems dealing with concrete structures.

There are several possibilities for coupling plasticity and damage effects in a single constitutive relation. Historically, damage has first been coupled to plasticity (see e.g. [25]) in the so-called ductile failure approaches for metal alloys. The underlying assumption was that void nucleation is triggered by plastic strains. Applications to concrete were proposed among others, by Oller et al. [33], Meschke and Lackner [31], Kratzig and Polling [21]. In these models, damage growth is a

<figure><img src="image_3.png" /><figcaption>Fig. 2. Unloading response of elastic damage (a), elastic plastic (b) and elastic plastic damage (c) models.</figcaption></figure>
function of the plastic strains. There is a difficulty, however. In uniaxial tension there is little plasticity and quite a lot of damage while in uniaxial compression, the picture is reversed with little damage and important plastic strains. Furthermore, it can be hardly explained how plastic strain may develop in concrete prior to microcracking. A common assumption is that irreversible strains are due to microcrack sliding and internal friction. Such a process requires the formation of internal surfaces (microcracks) prior to occur and develop.

The second approach, that is more suited to both tension and compression responses, uses the effective stress. The plastic yield function is no longer written in term of the applied stress. Rather, it is a function of the effective stress, that is the stress in the undamaged material in between the microcracks (see e.g. [42] and the discussion by Ju [20]). Among others Simo and Ju [43], Hansen and Schreyer [15], Yazdani and Schreyer [47], Fichant et al. [9], Mazars and Pijaudier-Cabot [29], Salari et al. [39], Lee and Fenves [24], Faria et al. [8], and Jefferson [19] applied this approach to isotropic and anisotropic damage coupled to elasto-plasticity. It has been extended to other sources of damage, for instance to thermal damage by Nechnech et al. [32].

A last possibility is what could be called the strong coupled approach. As opposed to the above where the plastic yield function is written in term of the effective stress, the applied stress appears in the plastic process, which becomes coupled to damage. Luccioni et al. [27] and Armero and Oller [2] provided the thermodynamic consistent grounds and discussed the algorithmic aspects of such a model. It has been also used by Gatuingt and Pijaudier-Cabot [10] for modelling the transient response of concrete.

We are aiming at a constitutive model that will be used for the computations of nuclear vessels. These are massive, heavily reinforced 3D pre-stressed concrete structures on which coupled hydro-mechanical computations should be performed in the future [11]. Here, we will deal only with the mechanical material response, but in the development of the constitutive equations, this final aim ought to be kept in mind. Hence, we look for a sufficiently simple constitutive model and we will choose the isotropic damage approach.

The major criteria that the model should satisfy are:

* a good description of the tensile material response from which cracking is expected. Such a load may occur when the vessel is pressurised,

* a good description of the material response for unconfined and confined compression as these are the major service loads,

* a good description of the compactant-dilatant volumetric response of concrete, e.g. in uniaxial compression.

This last issue is important because the dilatant material response induces transverse positive strains that may control damage growth (according to a strain driven process). Surprisingly, it is not very often investigated in studies on elasto-plastic damage or damage models. In fact very few of the above mentioned papers display the stress-volumetric strain curves and as we will see, the isotropic damage approach does not, alone, capture the volumetric material response correctly.

We have chosen here the effective stress approach because it provides a simple way to separate the damage and plastic processes. Plastic effects, driven by the effective stresses, can be described independently from damage ones and vice versa. One of the main interest is to ease the numerical implementation which is implicit/explicit. The plastic part is implicit, as usually, and the damage part is explicit, same as in classical continuum damage computations. As a consequence, existing robust algorithms for integrating the constitutive relations can be implemented. The calibration of the material parameters is also easier to handle as a consequence of the separation of damage and plasticity processes.

In this contribution, the damage process is (elastic) strain controlled. We will use here the classical isotropic model due to Mazars [28] which incorporates strain-softening. The plastic process shall be described with the help of a yield function inspired from Etse and Willam [7] and modified by Crouch and Tahar [5]. It is a hardening process in the present model. Softening is controlled by damage, while plasticity controls hardening, in compression especially. The model is discussed in Section 2. Section 3 deals with its algorithmic aspects. Elements of validation and comparisons between the damage and damage plasticity approaches are presented in Section 4. The constitutive law is applied to three elementary loading cases: simple tension, cyclic compression and triaxial compression. The axial capacity of a circular concrete-filled tube column, generally designed for seismic purposes, is also investigated. Applications to structures that are representative of a nuclear vessel have been presented in [18]. They are left, with the hydro-mechanical computations, for a future contribution.

## 2. Model formulation

In the following, we will restrict considerations to infinitesimal strains. The total strain tensor is denoted as $\epsilon$, the effective stress is $\sigma'$ and a standard, additive, elastic plastic strain decomposition ($\epsilon^e$ and $\epsilon^p$ respectively) is used. The isotropic damage variable is denoted as $D$ and the applied stress tensor is $\sigma$. The effective stress is classically defined as

$$ \sigma = (1 - D) \cdot \sigma'. \qquad (1) $$
## 2.1. Plasticity

Plasticity is governed by the following classical set of equations in which the effective stress has been substituted to the applied stress:

$$ 
\begin{aligned}
\dot{\boldsymbol{\epsilon}} &= \dot{\boldsymbol{\epsilon}}^c + \dot{\boldsymbol{\epsilon}}^p, \\
\boldsymbol{\sigma}' &= \boldsymbol{E}\boldsymbol{\epsilon}^c, \\
\dot{\boldsymbol{\epsilon}}^p &= \dot{\boldsymbol{\lambda}}\boldsymbol{m}(\boldsymbol{\sigma}', \boldsymbol{k}), \\
\dot{\boldsymbol{k}} &= \dot{\boldsymbol{\lambda}}\boldsymbol{h}(\boldsymbol{\sigma}', \boldsymbol{k}),
\end{aligned}
\qquad (2) $$

where $\boldsymbol{E}$ is the elastic stiffness tensor, $\boldsymbol{m}$ is the flow vector, $\boldsymbol{k}$ is the set of internal variables and $\boldsymbol{h}$ is the plastic modulus. The dot denotes time derivatives. The plastic multiplier $\dot{\boldsymbol{\lambda}}$ is given by the loading-unloading criterion (Kuhn-Tucker form)

$$ F(\boldsymbol{\sigma}', \boldsymbol{\kappa}) \le 0, \quad \dot{\boldsymbol{\lambda}} \ge 0, \quad F(\boldsymbol{\sigma}', \boldsymbol{\kappa}) \dot{\boldsymbol{\lambda}} = 0, \qquad (3) $$

where $F$ is the plastic yield function defined in the effective stress space. Note that we will use standard associated plasticity. It means in particular that

$$ \boldsymbol{m} = \frac{\partial F}{\partial \boldsymbol{\sigma}'} \qquad (4) $$

The yield function has been chosen according to two main objectives: (1) a correct volumetric behaviour is expected in uniaxial compression with a change from a contractant to a dilatant response; (2) an appropriate response in triaxial tests with confinement has to be obtained, with the apparition of plasticity for high levels of hydrostatic pressures especially. This supposes a closed function along the first invariant and eliminates Drucker-Prager equations [6]. The plastic yield surface depends on the three normalised effective stress invariant $(\bar{\rho}, \bar{\xi}, \theta)$ and one internal variable $k_h$

$$ \bar{\xi} = \frac{\operatorname{tr}(\boldsymbol{\sigma}')}{\sqrt{3}r_c}, \qquad \bar{\rho} = \frac{\sqrt{\boldsymbol{s}' : \boldsymbol{s}'}}{r_c}, \qquad \theta = \frac{1}{3} \arcsin \left( -\frac{\sqrt{3}}{2} \frac{\boldsymbol{s}' \cdot \boldsymbol{s}' \cdot \boldsymbol{s}'}{(\boldsymbol{s}' : \boldsymbol{s}')^{3/2}} \right), \qquad (5) $$

where $\boldsymbol{s}'$ is the deviatoric effective stress tensor, and $r_c$ is a parameter of the model. It is defined by means of three functions $\hat{k}$, $\bar{\rho}_c$ and $r$

$$ F = \bar{\rho}^2 - \frac{\hat{k}(\bar{\xi}, k_h)\bar{\rho}_c^2(\bar{\xi})}{r^2(\bar{\xi}, \theta)}. \qquad (6) $$

The hardening function $\hat{k}$, quantifies the growth of the loading surface with respect to the inelastic deformation. It is given by

$$ \hat{k} = k^{2p} \left( 1 - \frac{\bar{\xi}^2}{\bar{\xi}_h^2} \right), \qquad \bar{\xi}_h = \frac{A}{1-k}, \qquad k = k_0 + (1-k_0) \sqrt{k_h(2-k_h)} \qquad (7) $$

with $p$, $A$ and $k_0$ three parameters of the model. The deviatoric invariant $\bar{\rho}_c$ is defined as

$$ \bar{\rho}_c = \left(\frac{1}{6}\right)^\gamma \sqrt{\frac{2}{3}} \left(-n + \sqrt{n^2 - 12\sqrt{3}n\bar{\xi} + 36}\right)^\gamma, \qquad (8) $$

where

$$ n = \frac{3\left(1 - \bar{f}_t^2\right)}{\bar{f}_t + 2\bar{f}_t^\frac{\gamma}{2}}, \qquad \bar{f}_t = \frac{r_t}{r_c} \qquad (9) $$

and $\gamma$, $r_t$ and $r_c$ are three constants. Finally, $r$ is the deviatoric shape function given by [3], and it corresponds to a specific form of the elliptic surface developed in [46]

$$ r = \frac{2d_0}{d_1 - \sqrt{d_1^2 - 4d_0d_2}} \qquad (10) $$
with

$$
\begin{align}
d_0 &= c_1 \cos^2 \theta - c_2 \sin^2 \theta + c_3 \sin \theta \cos \theta, & c_0 &= \frac{(2 - \sqrt{3}B_1)(2B_0 - \sqrt{3}B_1)}{(B_1(1 + B_0) - \sqrt{3}B_0)^2}, \nonumber \\
d_1 &= 2(c_4\sqrt{3} \cos \theta - c_5 \sin \theta), & c_1 &= 3 - c_0(1 + B_0)^2, \nonumber \\
d_2 &= B_0(4 - 3B_0c_0), & c_2 &= 1 + 3c_0(1 - B_0)^2, \nonumber \\
B_0 &= \frac{\bar{\rho}_e}{\bar{\rho}_c}, & c_3 &= 2c_0\sqrt{3}(1 - B_0^2), \tag{11} \\
B_1 &= \frac{\sqrt{3}(1 - \alpha)B_0}{1 + B_0} + \frac{2\alpha B_0}{\sqrt{3}}, & c_4 &= (1 + B_0)(1 - B_0c_0), \nonumber \\
\bar{\rho}_e &= \left(\frac{1}{3}\right)^\gamma \sqrt{\frac{2}{3}} \left(-n + \sqrt{n^2 - 3\sqrt{3}n\bar{\xi} + 9}\right)^\gamma, & c_5 &= (1 - B_0)(1 - 3B_0c_0), \nonumber
\end{align}
$$

where $\alpha$ is a model parameter. Following the definition by Etse and Willam [7], the internal variable is a function of con-
finement too:

$$
\begin{align}
\dot{k}_h &= \frac{\sqrt{\frac{2}{3} \dot{\boldsymbol{\varepsilon}}^P : \dot{\boldsymbol{\varepsilon}}^P}}{\zeta(\boldsymbol{\sigma}')} \quad \text{if } k_h \le 1, \tag{12} \\
\dot{k}_h &= 0 \quad \text{if } k_h = 1, \nonumber
\end{align}
$$

where $\zeta$ depends on the first normalised invariant

$$
\begin{align}
\zeta &= -A_h + \sqrt{A_h^2 - B_h \bar{\xi} + C_h} \quad \text{if } \bar{\xi} \le 0, \tag{13} \\
\zeta &= -A_h + \sqrt{A_h^2 + C_h} \quad \text{if } \bar{\xi} > 0 \nonumber
\end{align}
$$

and $A_h$, $B_h$ and $C_h$ are three model parameters.

Note that Eq. (12) assumes that $k_h$ ranges between 0 and 1. When $k_h = 1$, the yield surface becomes a limit surface with no hardening in a pure elasto-plastic model. Overall, the plastic part of the constitutive relation relies on 10 model parameters. This is quite reasonable compared to some cap models which contain more than parameters (e.g. 15 in the MRS model [35]).

Combining Eqs. (12) and (2) yields the following expression for the plastic modulus $h$:

$$
\begin{align}
h &= \frac{\sqrt{\frac{2}{3} \frac{\partial F}{\partial \boldsymbol{\sigma}'} : \frac{\partial F}{\partial \boldsymbol{\sigma}'}}}{\zeta(\boldsymbol{\sigma}')} && \text{if } k_h \le 1, \tag{14} \\
h &= 0 && \text{if } k_h = 1. \nonumber
\end{align}
$$

Fig. 3(a) illustrates the evolution of the yield surface with the internal variable. Fig. 3(b) shows the plastic surface as a function of the third normalised invariant, the lode angle $\theta$ ranging from $-\frac{\pi}{6}$ (uniaxial tension) to $\frac{\pi}{6}$ (uniaxial compression) for $k_h = 0$. The plastic behaviour is not symmetric: depending on the value of $\theta$, the plastic threshold is not the same. Finally, Fig. 3(c) represents the failure surfaces for two lode angles, once the internal variable has reached the limited value of 1 and hardening is not allowed any more.

2.2. Damage model

The evolution of damage is defined according to the isotropic model developed by Mazars initially in (1984). It is slightly
modified here in order to account for plastic strains. These strains do not appear in the damage growth relationships. An
equivalent strain $\epsilon_{eq}$ is computed from the elastic strain tensor $\boldsymbol{\epsilon}^e$

$$
\epsilon^e = E^{-1} \sigma', \tag{15}
$$

where $E^{-1}$ is the inverse of the elastic tensor. This strain measure (or “equivalent strain”),

$$
\epsilon_{eq} = \sqrt{\sum_{i=1}^{3} (\langle e_i^e \rangle_{+})^2} \qquad (16)
$$
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Evolution of the yield surface with the hardening parameter for a constant lode angle (ε̄/ε₀), (b) evolution of the yield surface with the lode angle for a fixed hardening parameter and (c) failure surfaces for two lode angles.</figcaption></figure>

is used to control the damage evolution, and ⟨ε̄i⟩ represents the ith positive principal value of the elastic strain tensor. The loading surface g is defined as

$$g(\epsilon^e, D) = \tilde{d}(\epsilon^e) - D, \qquad (17)$$

where damage D is also the history variable. It takes the maximum value reached by d̃ during the history of loading

$$D = \text{Max}_{t}(\tilde{d}(t, 0)) \qquad (18)$$

and d̃ is defined by an evolution law which distinguishes the mechanical responses of the material in tension and in compression with the help of two pairs of scalars, one for tension (αt, Dt) and one for compression (αc, Dc). Namely,

$$\begin{align*}
\text{if } \epsilon_{eq} \ge \epsilon_{D_0}, \\
\tilde{d}(\epsilon^e) &= \alpha_t(\epsilon^e)D_t(\epsilon_{eq}) + \alpha_c(\epsilon^e)D_c(\epsilon_{eq}), \\
D_t &= 1 - \frac{\epsilon_{D_0}(1 - A_t)}{\epsilon_{eq}} - \frac{A_t}{\exp[B_t(\epsilon_{eq} - \epsilon_{D_0})]}, \tag{19} \\
\alpha_t &= \sum_{i=1}^{3} \frac{\epsilon_i^t \langle \epsilon_i^e \rangle_+}{\epsilon_{eq}^2}, \\
\text{else } \tilde{d}(\epsilon^e) &= 0.
\end{align*}$$

The parameters Dc and αc are determined by two expressions similar to those of Dt and αt; εD₀ represents the initial threshold from which damage grows. Dt and Dc are the tensile and compressive parts of the damage. At, Ac, Bt, Bc are four constants. The weights αt and αc are computed from the elastic strain tensor. They are defined as functions of the principal values of the strains εt and εe due to positive and negative effective stresses respectively (see [28]). In uniaxial tension, αt = 1 and αc = 0. In uniaxial compression, αt = 0 and αc = 1.

## 2.3. Model calibration

The complete elasto-plastic damage model contains 17 parameters, including the elastic ones. Parametric studies have been performed in order to check their sensitivity and to devise a calibration procedure. Following Crouch and Tahar [5], and from the sensitivity analysis, 4 parameters remain at a fixed value: α = 0.5, γ = 0.99, Ah = 7 × 10⁻⁵ and k₀ = 0.1.
$ε_{D_0}$, $A_t$, $B_t$ and $r_t$ are calibrated from a tensile test. $ε_{D_0}$ controls the peak stress (see [28]), $B_t$ and $A_t$ controls the shape of the softening regime. $r_t$ controls the amount of plastic strains and the residual stress for large strains. Cyclic compression provides $A_c$, $B_c$, $p$, $r_c$ and $B_h$. From the plot of the effective stress versus the plastic longitudinal strain (one point is extracted for each unloading cycle), the coefficients in the plastic model are obtained ($p$, $r_c$, $B_h$). The damage part is obtained from the plot of the equivalent strain versus damage. Finally, $A$ and $C_h$ are obtained from a confined cyclic compression test. $A$ is computed from the threshold of plastic strain growth, and $C_h$ is calibrated from the kinetics of plastic strain growth (same as in uniaxial compression).

# 3. Algorithmic aspects

## 3.1. Plastic process

As opposed to the damage constitutive equations which are explicit once the elastic strains are given, Eq. (2) in the plastic part of the model are path dependent and require an iterative process to be solved. These relations are integrated with the backward Euler scheme (between loading steps $t_{n+1}$ and $t_n$) according to a closest point projection. It yields the non-linear system of equations [41] in the effective stress space

$$ 
\begin{aligned}
&^{n+1}\boldsymbol{\sigma}' + \lambda \boldsymbol{E}m^{(n+1}\boldsymbol{\sigma}',^{n+1}k_h) = \boldsymbol{E}(^{n+1}\boldsymbol{\varepsilon} - {}^n\boldsymbol{\varepsilon}^p), \\
&^{n+1}k_h - \lambda h({}^{n+1}\boldsymbol{\sigma}',^{n+1}k_h) = {}^nk_h, \\
&F({}^{n+1}\boldsymbol{\sigma}',^{n+1}k_h) = 0,
\end{aligned} 
\quad (20) $$

where ${}^n\boldsymbol{\varepsilon}^p$, ${}^nk_h$ and the total strains ${}^{n+1}\boldsymbol{\varepsilon}$ are known at steps $n$ and $n+1$, respectively. The unknowns of the local problem are ${}^{n+1}\boldsymbol{\sigma}'$, ${}^{n+1}k_h$ and the incremental multiplier denoted as $\lambda$ in this section. An iterative Newton-Raphson method is used to solve this non-linear system of equations. The jacobian $J$ of the residual can be written as

$$ J = \begin{pmatrix} I_{\dim(\boldsymbol{\sigma}')} + \lambda \boldsymbol{E} \frac{\partial \boldsymbol{m}}{\partial \boldsymbol{\sigma}'} & \lambda \boldsymbol{E} \frac{\partial \boldsymbol{m}}{\partial k_h} & \boldsymbol{E}m \\ -\lambda \frac{\partial h}{\partial \boldsymbol{\sigma}'} & 1 - \lambda \frac{\partial h}{\partial k_h} & -h \\ \boldsymbol{m}^T & \left(\frac{\partial F}{\partial k_h}\right)^T & 0 \end{pmatrix}, \quad (21) $$

where superscript T denotes for the transpose operator. The derivatives of the flow vectors are computed according to a numerical differentiation technique described in [35]. Note that the following initial values are prescribed for the three unknowns (elastic predictor):

$$ 
\begin{aligned}
&^{n+1} \boldsymbol{\sigma}' = {}^n \boldsymbol{\sigma}' + \boldsymbol{E}(^{n+1} \boldsymbol{\varepsilon} - {}^n \boldsymbol{\varepsilon}), \\
&^{n+1} k_h = {}^n k_h, \\
&\lambda = 0.
\end{aligned} 
\quad (22) $$

## 3.2. Plastic-damage process

The full integration of the constitutive law is done in two main steps as depicted in Fig. 4. The state at step $n$ (${}^n\boldsymbol{\varepsilon}$, $^{n}\boldsymbol{\sigma}$) and the total strain ${}^{n+1}\boldsymbol{\varepsilon}$ at step ($n+1$) and global iteration $k$ are known. An effective stress ${}^{n+1}\boldsymbol{\sigma}'$ is then computed as sketched above. The elastic plastic strain decomposition ${}^{n+1}\boldsymbol{\varepsilon}^e$, ${}^{n+1}\boldsymbol{\varepsilon}^p$ follows and damage is evaluated (see Section 2.2) from the elastic strain explicitly. Finally, the total stress ${}^{n+1}\boldsymbol{\sigma}$ is computed from the value of damage and of the effective stress, see Eq. (1).

For the global equilibrium, solved according to a Newton-Raphson algorithm, a consistent tangent operator is computed. It derives from the elasto-plastic operator obtained from Eq. (1)

$$ \frac{\partial \boldsymbol{\sigma}}{\partial \boldsymbol{\varepsilon}} = (1 - \boldsymbol{D}) \frac{\partial \boldsymbol{\sigma}'}{\partial \boldsymbol{\varepsilon}} - \frac{\partial \boldsymbol{D}}{\partial \boldsymbol{\varepsilon}} \boldsymbol{\sigma}', \quad (23) $$

where the elasto-plastic tangent operator appears in the first term on the right hand side of this equation. Since damage depends on the elastic strain only, $\frac{\partial \boldsymbol{D}}{\partial \boldsymbol{\varepsilon}}$ can be written as

$$ \frac{\partial \boldsymbol{D}}{\partial \boldsymbol{\varepsilon}} = \frac{\partial \boldsymbol{D}}{\partial \boldsymbol{\varepsilon}^e} \frac{\partial \boldsymbol{\varepsilon}^e}{\partial \boldsymbol{\varepsilon}}. \quad (24) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. General algorithm for the elastic plastic damage model.</figcaption></figure>

The derivatives of the elastic strain with respect to the total strain are thus computed using Eq. (2)

$$ \frac{\partial \boldsymbol{\epsilon}^e}{\partial \boldsymbol{\epsilon}} = \boldsymbol{E}^{-1} \frac{\partial \boldsymbol{\sigma}'}{\partial \boldsymbol{\epsilon}}. \qquad (25) $$

For the derivatives of the effective stress with respect to the strains, Eqs. (2) are linearized following [35]. After simplifications, we have

$$ \left( \frac{d\boldsymbol{\sigma}'}{dk_h} \right) = \left[ A^{-1} - A^{-1} \begin{bmatrix} n\boldsymbol{m} \\ -n h \end{bmatrix} \frac{\left[ n \left( \frac{\partial F}{\partial \boldsymbol{\sigma}'} \right) \quad n \left( \frac{\partial F}{\partial k_h} \right) \right] A^{-1}}{\left[ n \left( \frac{\partial F}{\partial \boldsymbol{\sigma}'} \right) \quad n \left( \frac{\partial F}{\partial k_h} \right) \right] A^{-1} \begin{pmatrix} n\boldsymbol{m} \\ -n h \end{pmatrix}} \right] \begin{bmatrix} d\boldsymbol{\epsilon} \\ 0 \end{bmatrix} \qquad (26) $$

with

$$ A = {}_n\left( \begin{array}{cc} E^{-1} + \lambda \left( \frac{\partial \boldsymbol{m}}{\partial \boldsymbol{\sigma}'} \right) & \lambda \left( \frac{\partial \boldsymbol{m}}{\partial k_h} \right) \\ -\lambda \left( \frac{\partial h}{\partial \boldsymbol{\sigma}'} \right) & 1 - \lambda \left( \frac{\partial h}{\partial k_h} \right) \end{array} \right). \qquad (27) $$

$\frac{\partial \boldsymbol{\sigma}'}{\partial \boldsymbol{\epsilon}}$ is then computed by extracting the first 6 × 6 block of the matrix in Eq. (26).

In order to calculate the derivatives of the damage with respect to the elastic strains, a radial loading is assumed, same as in [28]

$$ d\alpha_t = d\alpha_c = 0 \qquad (28) $$

then

$$ \frac{\partial D}{\partial \boldsymbol{\epsilon}^e} = \left( \alpha_t \frac{\partial D_t(\boldsymbol{\epsilon}^e)}{\partial \varepsilon_{eq}} + \alpha_c \frac{\partial D_c(\boldsymbol{\epsilon}^e)}{\partial \varepsilon_{eq}} \right) \cdot \frac{\partial \varepsilon_{eq}}{\partial \boldsymbol{\epsilon}^e}. \qquad (29) $$

From [23], a simplified expression of the derivative of the damage with respect to the elastic strains can be given

$$ \frac{\partial D}{\partial \varepsilon_{kl}^e} = - \left[ \alpha_t \left( \frac{\varepsilon_{D_0} (1 - A_t)}{\varepsilon_{eq}^2} + \frac{A_t B_t}{\exp(B_t (\varepsilon_{eq} - \varepsilon_{D_0}))} \right) + \alpha_c \left( \frac{\varepsilon_{D_0} (1 - A_c)}{\varepsilon_{eq}^2} + \frac{A_c B_c}{\exp(B_c (\varepsilon_{eq} - \varepsilon_{D_0}))} \right) \right] \cdot \frac{\varepsilon_{kl}^{e+}}{\varepsilon_{eq}} \cdot \beta_{kl}, \qquad (30) $$

where $\varepsilon^{e+}$ is the positive elastic stain (spectral decomposition) and $\beta_{kl} = 2 * (1 - \delta_{kl})$ with $\delta_{kl}$ the Kronecker symbol.

Combining Eqs. (23)–(26) and (30), we obtain a consistent tangent operator for the elastic plastic damage formulation.
## 4. Model validation and comparison with a damage model

The constitutive law is examined first using three “elementary” test cases: simple tension, cyclic compression and triaxial loading with confinement. Next, the performance of the model is evaluated in a simple structural problem: a column made of a steel tube filled with concrete. Note that the model parameters will be adjusted to each formulation of concrete as the different tests have not been performed on the same material.

### 4.1. Uniaxial tension test

The numerical response of the elastic plastic damage law is first compared with experiments from [13]. The parameters chosen for this simulation are shown in Table 1. Fig. 5(a) gives the axial stress-strain curve.

For the sake of comparisons, a damage model is also considered (Fig. 5(b)). In this case, plastic strains are equal to zero; thus the elastic strain is equal to the total strain ($\epsilon = \epsilon^c$). The response coincides with the original damage model developed in [28]. The coefficients used for this problem are reported in Table 2. The two sets of parameters in the damage and plastic damage models are not the same because one has to take into account the plasticity effect in the second.

When the development of damage is predominant, both models are able to reproduce the monotonic experiment. Nevertheless, with the elastic plastic damage formulation, a constant stress value is obtained for a certain strain ($\epsilon \approx 0.00036$). At this loading level, the plastic limit surface has been reached. The effective response is perfectly plastic and no further damage can evolve due to a constant elastic strain. Even if this constant stress value is not dramatic for this example, it represents one limit of the proposed formulation. One improvement could be to propose another evolution law for the hardening variable, without a limit surface (or assuming that the limit surface is reached upon an infinite strain). It would enable the plastic loading surface to expand indefinitely, and therefore the elastic strain to increase under monotonically increasing applied strain.

### 4.2. Cyclic uniaxial compression

Cyclic compression is the second elementary test used to highlight the interest of the model. Experimental results are taken from [44]. Fig. 6(a) illustrates the numerical response with the pure damage law (without plasticity). The coefficients chosen for this simulation are reported in Table 3. They are different from those in Section 4.1 because the concrete is different too. With this type of relation, a zero stress corresponds to a zero strain. No irreversible effect is simulated.

<table><caption>Table 1
Damage and plastic parameters for simple tension test with the elastic plastic damage model</caption><thead><tr><th>E (Pa)</th><th>ν</th><th>A<sub>c</sub></th><th>A<sub>t</sub></th><th>B<sub>c</sub></th><th>B<sub>t</sub></th><th>ε<sub>D<sub>0</sub></sub></th><th>r<sub>c</sub></th><th>r<sub>t</sub></th><th>p</th><th>B<sub>h</sub></th><th>C<sub>h</sub></th><th>A</th></tr></thead><tbody><tr><td>3.125 × 10<sup>10</sup></td><td>0.2</td><td>2.75</td><td>1.1</td><td>2470</td><td>13,000</td><td>1.1 × 10<sup>-4</sup></td><td>120 × 10<sup>6</sup></td><td>11 × 10<sup>6</sup></td><td>0.4</td><td>2 × 10<sup>-2</sup></td><td>3 × 10<sup>-6</sup></td><td>-1</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Simulation of a uniaxial tension test [13]: (a) elastic plastic damage formulation; (b) elastic damage model.</figcaption></figure>

<table><caption>Table 2
Damage parameters for simple tension test with the elastic damage model</caption><thead><tr><th>E (Pa)</th><th>ν</th><th>A<sub>c</sub></th><th>A<sub>t</sub></th><th>B<sub>c</sub></th><th>B<sub>t</sub></th><th>ε<sub>D<sub>0</sub></sub></th></tr></thead><tbody><tr><td>3.125 × 10<sup>10</sup></td><td>0.2</td><td>2.75</td><td>0.88</td><td>2470</td><td>8000</td><td>0.6 × 10<sup>-4</sup></td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Simulation of a cyclic compression test [44]. Elastic damage response.</figcaption></figure>

<table><caption>Table 3
Damage parameters for cyclic compression test with the elastic damage model</caption><thead><tr><th>E (Pa)</th><th>ν</th><th>A<sub>c</sub></th><th>A<sub>t</sub></th><th>B<sub>c</sub></th><th>B<sub>t</sub></th><th>ε<sub>D0</sub></th></tr></thead><tbody><tr><td>3.2 × 10<sup>10</sup></td><td>0.2</td><td>1.15</td><td>1.0</td><td>1391</td><td>10,000</td><td>0.9375 × 10<sup>-4</sup></td></tr></tbody></table>

The response of the elastic plastic damage model is given in Fig. 7a, and the coefficients are listed in Table 4. $r_c$, $p$ and $B_h$ are the most influent constants for this test. They quantify the evolution of the plastic strains, especially at the unloading points. With this model, damage induces the global softening behaviour while plasticity reproduces quantitatively the evolution of the irreversible strains. Experimental and numerical unloading slopes are now similar, contrary to the standard damage formulation response. This difference is essential when an accurate value of the damage needs to be estimated. The elastic damage model overestimates $D$ whereas the full constitutive law provides more realistic results.

Figs. 6(b) and 7(b) illustrate the differences between both the damage and the plastic damage models in terms of volumetric responses. Notice that the isotropic damage law induces a contractant response (negative volumetric strains). The introduction of plasticity generates a change in the volumetric response from contractant to dilatant, a phenomenon which is experimentally observed (see [40] for example).

Plasticity combined with the evolution of isotropic damage plays thus a key role in the numerical simulation of a cyclic compression test. The development of irreversible strains during loading is quantitatively reproduced and the softening behaviour fits well. Moreover, the volumetric response, that was completely misevaluated by the damage model, is correctly simulated by the full formulation. This misevaluation is also due to the fact that damage is assumed to be isotropic. It constrains the Poisson ratio to remain constant. A more accurate volumetric response could be achieved with anisotropic

<figure><img src="image_8.png" /><figcaption>Fig. 7. Simulation of a cyclic compression test [44]. Elastic plastic damage response.</figcaption></figure>

<table><caption>Table 4
Damage and plastic parameters for cyclic compression test with the elastic plastic damage model</caption><thead><tr><th>E (Pa)</th><th>ν</th><th>A<sub>c</sub></th><th>A<sub>t</sub></th><th>B<sub>c</sub></th><th>B<sub>t</sub></th><th>ε<sub>D0</sub></th><th>r<sub>c</sub></th><th>r<sub>t</sub></th><th>p</th><th>B<sub>h</sub></th><th>C<sub>h</sub></th><th>A</th></tr></thead><tbody><tr><td>3.125 × 10<sup>10</sup></td><td>0.2</td><td>2.75</td><td>1.1</td><td>2470</td><td>15,000</td><td>1 × 10<sup>-4</sup></td><td>120 × 10<sup>6</sup></td><td>10 × 10<sup>6</sup></td><td>0.4</td><td>2 × 10<sup>-2</sup></td><td>3 × 10<sup>-6</sup></td><td>-1</td></tr></tbody></table>
damage approaches, nevertheless irreversible strains would have to be captured with some sort of plastic based mechanism. The resulting constitutive relations would be much more complex (coupled anisotropic damage and plasticity).

### 4.3. Triaxial test with confinement

The experimental results from [40] are used to evaluate the ability of the constitutive law to reproduce triaxial tests after hydrostatic loading. A vertical displacement is applied on a concrete cylinder (150 × 300 mm specimens) after an initial confinement. Numerical results are compared with experiments for different levels of hydrostatic pressures (P = 0, 1.5, 4.5, 9, 30 and 60 MPa).

Fig. 8(a) gives the axial response (the confinement phase is not represented) for the first four pressures with the full formulation using the coefficients of Table 5. Again some parameters are different from those listed in Sections 4.1 or 4.2 or because the concrete properties are not similar. The peaks are quantitatively well reproduced (except for 1.5 MPa). The global evolutions are also correct: the maximum axial stress increases with the pressure and softening is less and less significant. When the initial confinement takes higher values, damage plays a minor role and plasticity becomes predominant.

Fig. 8(b) presents the axial curves for 30 and 60 MPa for the same model parameters. Experimental results and simulations are in agreement. The decrease of the initial stiffness upon confinement is reproduced by the constitutive law especially. For 30 MPa for instance, non-linearity has already been initiated in the hydrostatic load regime, before the application of the vertical displacement begins. The introduction of plasticity and the characteristic shape of the yield surface (closed function along the volumetric invariant) provide a good description of these experiments.

Fig. 9(a) shows a comparison between the elastic plastic damage formulation and the damage model for the high confinement pressures. The standard damage law fails to reproduce the decrease of the initial stiffness. Note that as soon as the pressure takes important values (30 or 60 MPa in the figure), the damage model clearly overestimates the experimental behaviour. In fact, this is due to the definition of the equivalent strain $\varepsilon_{eq}$ (16) that characterises the material extension during loading in the damage model. When the hydrostatic pressure is applied, the sample is not subjected to tensile strains, the equivalent strain keeps a zero value and the material response is elastic, if no plastic effects are taken into account. This is acceptable for low confinement pressures; but it is no longer reasonable when higher levels are considered. Fig. 9(b) shows the two responses, according to the damage and plastic damage models, during the application of the hydrostatic pressure. As expected, the response is linear with the damage model. Fig. 10 shows the evolution of the computed transverse strains for low confinement pressures with the full plastic damage formulation and a comparison with experiment. Even though the transverse strains are underestimated at the beginning of the loading, the comparison with experiments is globally correct.

A better response might be achieved using a non-associated flow rule probably. Nevertheless, this improvement would require an increasing number of parameters, a more complex numerical implementation and has not been considered to be necessary: the associated plastic model seems indeed sufficient in our case to reproduce the volumetric behaviour globally

<figure><img src="image_9.png" /><figcaption>Fig. 8. Triaxial test with increasing confinement. (a) Axial stress–strain curves for low hydrostatic pressures. Straight lines (black markers) correspond to simulation, dotted lines (white markers) to experiment (1) 0 MPa, (2) 1.5 MPa, (3) 4.5 MPa, (4) 9 MPa. (b) Axial stress–strain curves for high hydrostatic pressures.</figcaption></figure>

<table><caption>Table 5
Damage and plastic parameters for triaxial confinement test with the elastic plastic damage model</caption><thead><tr><th>E (Pa)</th><th>ν</th><th>A<sub>c</sub></th><th>A<sub>t</sub></th><th>B<sub>c</sub></th><th>B<sub>t</sub></th><th>ε<sub>D<sub>0</sub></sub></th><th>r<sub>c</sub></th><th>r<sub>t</sub></th><th>p</th><th>B<sub>h</sub></th><th>C<sub>h</sub></th><th>A</th></tr></thead><tbody><tr><td>2.73 × 10<sup>10</sup></td><td>0.2</td><td>2.75</td><td>1.1</td><td>2600</td><td>15,000</td><td>1 × 10<sup>-4</sup></td><td>85 × 10<sup>6</sup></td><td>9 × 10<sup>6</sup></td><td>0.4</td><td>3 × 10<sup>-3</sup></td><td>3 × 10<sup>-6</sup></td><td>-0.5</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison between elastic damage and elastic plastic damage constitutive laws: (a) triaxial tests with increasing confinement. Axial stress–strain curves for high hydrostatic pressures; comparison between the elastic damage model (dam) and the elastic plastic damage formulation (plas); (b) hydrostatic confinement. Comparison between elastic damage model (pl dam) and elastic damage constitutive law (dam).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Triaxial test with increasing confinement. Axial stress versus transverse strain for 1.5, 4.5 and 9 MPa.</figcaption></figure>

and especially the change from a contractant to a dilatant evolution (see Section 4.2). It appears that the overall response is not too dilatant, as opposed to classical results in standard associated plasticity (at least for the experiments investigated in this paper). This may be due to the combination of damage with plasticity, which seems to avoid excessive dilatancy.

## 4.4. Circular concrete-filled tube column

Composite structures are widely used in building and other civil engineering structures, in seismic zones especially. Among them, concrete-filled steel tubes (CFT) present some interesting resisting characteristics such as high ductility and improved strength. They have additional advantages compared to pure steel or simple concrete members: concrete increases the stability of the column while the steel tube subjects the core concrete to a triaxial state of stress and induces passive confinement [26]. Various experimental and numerical studies have been recently carried out on such structures ([12,22] for example).

To highlight the interest of our constitutive law, the behaviour of a circular CFT is going to be simulated. The dimensions of the sample and the mechanical properties experimentally reported in [45] and measured on non-wrapped specimens, are listed in Table 6. The steel–concrete interface is assumed to be perfect. For the considered compressive strength $f'_c$, Giakoumelis and Lam [12] demonstrated, with greased and non-greased cylinders, that the steel–concrete interface has little influence on the global behaviour.

A vertical displacement (yielding uniaxial compression) is applied on both steel and concrete cross-sections at the same time. Two simulations are proposed. One using the standard elastic plastic damage formulation (Table 7) and another one employing the elastic damage constitutive law (Table 8). The steel tube is modelled with a Von Mises plasticity model with the yield stress $\sigma_e = 279.9$ MPa and the tangent hardening modulus $E_t = 2500$ MPa. One fourth of the cylinder is meshed with 99 eight node cubic elements.

Fig. 11(a) and (b) provide a comparison between the simulations with the two models and with the experiment (here localisation of damage is not expected and comparisons with the experiment are more relevant). With the elastic plastic
<table><caption>Table 6
Geometry and mechanical properties of the circular CFT as experimentally reported in [45]</caption><thead><tr><th rowspan="2"></th><th colspan="2">Geometry (mm)</th><th rowspan="2">L</th><th colspan="2">Steel properties</th><th colspan="2">Concrete properties</th></tr><tr><th>D</th><th>e</th><th>E<sub>s</sub> (Pa)</th><th>σ<sub>e</sub> (MPa)</th><th>E (Pa)</th><th>f'<sub>c</sub> (MPa)</th></tr></thead><tbody><tr><td>L</td><td>150</td><td>4</td><td>450</td><td>21 × 10<sup>10</sup></td><td>279.9</td><td>2.18 × 10<sup>10</sup></td><td>22</td></tr></tbody></table>

<table><caption>Table 7
Damage and plastic parameters for the CFT</caption><thead><tr><th>E (Pa)</th><th>ν</th><th>A<sub>c</sub></th><th>A<sub>t</sub></th><th>B<sub>c</sub></th><th>B<sub>t</sub></th><th>ε<sub>D<sub>0</sub></sub></th><th>r<sub>c</sub></th><th>r<sub>t</sub></th><th>p</th><th>B<sub>h</sub></th><th>C<sub>h</sub></th><th>A</th></tr></thead><tbody><tr><td>2.18 × 10<sup>10</sup></td><td>0.2</td><td>2.3</td><td>1.1</td><td>2700</td><td>15,000</td><td>0.95 × 10<sup>-4</sup></td><td>85 × 10<sup>6</sup></td><td>9 × 10<sup>6</sup></td><td>0.3</td><td>5 × 10<sup>-2</sup></td><td>3 × 10<sup>-6</sup></td><td>-0.5</td></tr></tbody></table>

<table><caption>Table 8
Damage parameters for the CFT</caption><thead><tr><th>E (Pa)</th><th>ν</th><th>A<sub>c</sub></th><th>A<sub>t</sub></th><th>B<sub>c</sub></th><th>B<sub>t</sub></th><th>ε<sub>D<sub>0</sub></sub></th></tr></thead><tbody><tr><td>2.18 × 10<sup>10</sup></td><td>0.2</td><td>1.15</td><td>0.8</td><td>1391</td><td>10,000</td><td>0.95 × 10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Simulation of a CFT: (a) evolution of the axial force as a function of the axial strain. Comparison between the elastic plastic damage formulation and the experiment; (b) comparison between the elastic plastic damage and the elastic damage models.</figcaption></figure>

damage relation, numerical and experimental axial forces are in agreement for a given axial strain. The “confinement” effect is highlighted with an increase of the maximal strength. On the contrary, the damage model underestimates the global response of the column. A softening branch appears, which is not observed during experiments. Note that the overall elastic stiffness does not fit exactly, same as for the simulations performed in [45]. This is probably due to experimental differences between the elastic mechanical properties measured on wrapped and non-wrapped samples due to incomplete hydration of concrete for example, see Kwon and Spacone [22] for details.

Experimentally, no confinement effect is noticed at the beginning of the loading. The transverse strain in concrete is lower than in steel due to differences in the Poisson ratio (0.2 and 0.3, respectively). According to the computation, concrete is thus under lateral tension because the steel concrete bond is perfect, but this regime seems to have little influence on the structural response. As the axial load increases, plasticity is responsible for the change of the apparent Poisson ratio. Lateral expansion in concrete becomes gradually higher than in the steel tube. A radial pressure develops at the interface and passive confinement appears. The evolution of the radial stress at the interface as a function of the axial strain is provided in Fig. 12(a) using the full formulation. Fig. 12(b) illustrates the evolution of the transverse strains $ε_c$ and $ε_s$ in concrete and in steel, respectively. As observed experimentally, concrete is first subjected to lateral extension then to lateral compression. The change of sign occurs immediately when $ε_c$ becomes higher than $ε_s$. On the contrary, with the simple elastic damage constitutive law, as $ε_c$ is always lower than $ε_s$, concrete is only subjected to lateral extension (Fig. 13). No passive confinement is observed and that is why the peak in the axial load is so small and inadequate compared with experiment. The study of the volumetric behaviour yields the same conclusions (Fig. 14). The change from a contractant to a dilatant response obtained with the introduction of plasticity is a direct consequence of the increase of the concrete transverse strains. With the elastic damage relation, the volumetric response is always contractant.
<figure><img src="image_13.png" /></figure>

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

<figure><img src="image_15.png" /><figcaption>Fig. 12. Transverse behaviour using the elastic plastic damage model: (a) evolution of the concrete transverse stress as a function of the axial strain; (b) evolution of the transverse concrete and steel transverse strains as a function of the axial strain and (c) zoom on the first part of the curve.</figcaption></figure>

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

<figure><img src="image_17.png" /><figcaption>Fig. 13. Transverse behaviour using the elastic damage model: (a) evolution of the transverse stress as a function of the axial strain; (b) evolution of the concrete and steel transverse strains as a function of the axial strain.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Evolution of the volumetric strain as a function of the axial load for both approaches.</figcaption></figure>
For concrete-filled steel tubes subjected to axial compressive loading, a correct description of the volumetric response is required in order to capture the change of apparent Poisson ratio of concrete that induces passive confinement. The elastic plastic damage constitutive law may also be required for pre-stressed structures where the material is subjected to biaxial confinement for instance.

## 5. Conclusions

An elastic plastic damage constitutive law has been presented and compared with a damage model on different types of simple loadings and structural test cases. In this model, isotropic damage is responsible for the softening response and the decrease in the elastic stiffness, while hardening plasticity accounts for the development of irreversible strains and volumetric compressive behaviour.

In compressive tests, improvements compared with the standard elastic isotropic damage relation have been achieved. Irreversible strains and the volumetric response of the material in uniaxial compression have been quantitatively well reproduced. The contractant-dilatant volumetric response is due to plasticity. With a scalar isotropic damage model, it cannot be captured. In triaxial tests, the decrease of the initial slope and brittle-ductile transition are also captured with the plastic damage formulation, as opposed to the damage formulation which fails for moderate to high confinement levels. Finally, the example of a concrete-filled steel tube has highlighted the importance of the volumetric response of concrete in order to capture the appropriate structural response with such a passive confinement.

Tension dominant applications have not been considered in this paper. Due to strain and damage localisation it is well established that simulated responses are expected to be mesh dependent (see for instance [37] or 17). Qualitative comparisons between the damage and plastic damage approaches have been carried out by Jason [16]. They show that the introduction of plasticity does not change the failure mode of the structure compared to usual results with a simple damage model because in tension dominated situations the damage part of the model is the essential ingredient of the constitutive response. The plastic damage model has the same modelling capabilities as the damage model in a type of structural analysis where it is established that the second performs very well. In the case of regular meshes, it is possible to adjust the softening response of the model to fit experimental data on structural applications (provided the discretisation is such that mesh alignment is consistent with the crack path). In a more general context, however, regularisation techniques should be included in the formulation to avoid numerical problems [37,34,38,1].

The next step is the durability analysis of containment vessels and the implementation of this constitutive relation in representative model structures [11]. This has been performed in Jason et al. [18]. The hydro-mechanical problem, from which the leakage rate upon damage for such structures can be computed, is currently being investigated.

## Acknowledgements

Partial financial supports from EDF and from the EU through MAECENAS project (FIS5-2001-00100) are gratefully acknowledged. The authors particularly thank R. Crouch (University of Sheffield) for his help in the design and the numerical development of the plasticity model. The authors would like to thank EDF for scientific support toward the developments in the FE code "Code_Aster" and the Spanish-French sponsorship HF2003-0041.

## References

[1] D. Addessi, S. Marfia, E. Sacco, A plastic nonlocal damage model, Comput. Methods Appl. Mech. Engrg. 191 (2002) 1291-1310.

[2] F. Armero, S. Oller, A general framework for continuum damage models. I. Infinitesimal plastic damage models in stress space, Int. J. Solids Struct. 37 (2000) 7409-7436.

[3] S.K. Bhowmik, J.H. Long, A general formulation for the cross-sections of yield surfaces in octahedral planes, in: G.N. Pande, J. Middleton (Eds.), Numeta 90, Elsevier, 1990, pp. 795-803.

[4] G. Chatzigeorgiou, V. Picandet, A. Khelidj, G. Pijaudier-Cabot, Coupling between progressive damage and permeability of concrete: analysis with a discrete model, Int. J. Numer. Anal. Methods Geomech. 29 (2005) 1005-1018.

[5] R.S. Crouch, B. Tahar, Application of a stress return algorithm for elasto-plastic hardening-softening models with high yield surface curvature, in: Proceedings of European Congress on Computational Methods in Applied Sciences and Engineering, Barcelona, 2000.

[6] D.C. Drucker, W. Prager, Soil mechanics and plasticity analysis of limit design, Quart. J. Appl. Mech. 10 (1952) 157-162.

[7] G. Etse, K.J. Willam, Fracture energy formulation for inelastic behavior of plain concrete, J. Engrg. Mech. ASCE 106 (9) (1994) 1013-1203.

[8] R. Faria, J. Oliver, M. Cervera, A strain-based plastic viscous damage model for massive concrete structures, Int. J. Solids Struct. 35 (1998) 1533-1558.

[9] S. Fichant, C. La Borderie, G. Pijaudier-Cabot, Isotropic and anisotropic descriptions of damage in concrete structures, Int. J. Mech. Cohesive Frictional Mater. 4 (1999) 339-359.

[10] F. Gatuingt, G. Pijaudier-Cabot, Coupled damage and plasticity modelling in transient dynamic analysis of concrete, Int. J. Numer. Anal. Methods Geomech. 26 (2002) 1-24.

[11] S. Ghavamian, Evaluation tests on models of non linear behaviour of cracking concrete using three dimensional modeling, Benchmark EDF R&D, CR-MMN 99/232, 1999.
[12] G. Giakoumelis, D. Lam, Axial capacity of circular concrete-filled tube columns, J. Construc. Steel Res. 60 (2004) 1049–1068.

[13] V.S. Gopalaratnam, S.P. Shah, Softening response of plain concrete in direct tension, ACI J. (May–June) (1985) 310–323.

[14] P. Grassl, K. Lundgren, K. Gylltoft, Concrete in compression: a plasticity theory with a novel hardening law, Int. J. Solids Struct. 39 (2002) 5205–5223.

[15] N.R. Hansen, H.L. Schreyer, Thermodynamically consistent theories for elastoplasticity coupled with damage, in: J.W. Ju et al. (Eds.), Damage Mechanics and Localisation, AMD vol. 142/MD vol. 34, ASME Publications, 1992, pp. 53–67.

[16] L. Jason, Relation endomanagement—perméabilité pour les bétons—application au calcul de structures, Ph.D. Dissertation, Ecole Centrale de Nantes, France, 2004.

[17] L. Jason, S. Ghavamian, G. Pijaudier-Cabot, A. Huerta, Benchmarks for the validation of a non local damage model, Rev. Franc. Génie Civ. 8 (2004) 303–328.

[18] L. Jason, G. Pijaudier-Cabot, A. Huerta, S. Ghavamian, Damage and plasticity for concrete behaviour, in: P. Neittaanmäki, T. Rossi, S. Korotov, E. Onate, J. Périaux, D. Knörzer (Eds.), Proc. ECCOMAS 2004, Jyväskylä, Finland, 2004.

[19] A.D. Jefferson, Craft—a plastic damage contact model for concrete—I. Model theory and thermodynamic considerations, Int. J. Solids Struct. 40 (2003) 5973–5999.

[20] J.W. Ju, Energy-based coupled elastoplastic damage models at finite strains, J. Engrg. Mech. ASCE 115 (1989) 2508–2526.

[21] W. Kratzig, R. Polling, An elasto-plastic damage model for reinforced concrete with minimum number of material parameters, Comput. Struct. 82 (2004) 1201–1215.

[22] M. Kwon, E. Spacone, Three-dimensional finite element analyses of reinforced concrete columns, Comput. Struct. 80 (2002) 199–212.

[23] P. Ladevèze, On an anisotropic damage theory, in: A. Boehler et al. (Eds.), Proceedings Failure Criteria of Structured Media, in Failure Criteria of Structured Media, Balkema, 1993, pp. 355–363.

[24] J. Lee, G.L. Fenves, Plastic damage models for cyclic loading of concrete structures, J. Engrg. Mech. 124 (8) (1998) 892–900.

[25] J. Lemaitre, J.L. Chaboche, Mécanique des Matériaux Solides, Dunod, Paris, 1984.

[26] Z. Lu, Y. Chen, Numerical simulation and research on hysteretic model for concrete-filled square tube beam-columns considering local buckling, Steel Struct. 2 (2002) 41–50.

[27] B. Luccioni, S. Oller, R. Danesi, Coupled plastic damaged model, Comput. Methods Appl. Mech. Engrg. 129 (1996) 81–89.

[28] J. Mazars, Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure, Ph.D. Thesis, Université Paris VI, 1984.

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

[30] Ph. Menétrey, K.J. Willam, Triaxial failure criterion for concrete and its generalization, ACI Struct. J. 92 (3) (1995) 311–318.

[31] G. Meschke, R. Lackner, Anisotropic modelling of cracked concrete based on plasticity-damage theory, in: D.R.J. Owen et al. (Eds.), Computational Plasticity, CIMNE, Barcelona, 1997, pp. 1543–1550.

[32] W. Nechnech, F. Meftah, J.M. Reynouard, An elasto-plastic damage model for plain concrete subjected to high temperatures, Engrg. Struct. 24 (2002) 597–611.

[33] S. Oller, E. Onate, J. Oliver, J. Lubliner, Finite element nonlinear analysis of concrete structures using a plastic damage model, Engrg. Fract. Mech. 35 (1990) 219–231.

[34] R.H.J. Peerlings, R. de Borst, W.A.M. Brekelmans, J.H.P. De Vree, Gradient enhanced damage for quasi brittle materials, Int. J. Numer. Methods Engrg. 39 (1996) 3391–3403.

[35] A. Pérez-Foguet, A. Rodriguez-Ferran, A. Huerta, Numerical differentiation for non trivial consistent tangent matrices: an application to the MRS-Lade model, Int. J. Numer. Methods Engrg. 48 (2000) 159–184.

[36] V. Picandet, A. Khelidj, G. Bastian, Effect of axial compressive damage on gas permeability of ordinary and high performance concrete, Cem. Concr. Res. 31 (2001) 1525–1532.

[37] G. Pijaudier-Cabot, Z.P. Bazant, Non local damage theory, J. Engrg. Mech. ASCE 113 (1987) 1512–1533.

[38] A. Rodriguez-Ferran, I. Morata, A. Huerta, Efficient and reliable nonlocal damage models, Comput. Methods Appl. Mech. Engrg. 193 (2004) 3431–3455.

[39] M.R. Salari, S. Saeb, K.J. Willam, S.J. Patchet, R.C. Carrasco, A coupled elastoplastic damage model for geomaterials, Comput. Methods Appl. Mech. Engrg. 193 (2004) 2625–2643.

[40] D. Sfer, I. Carol, R. Gettu, G. Etse, Study of the behavior of concrete under triaxial compression, J. Engrg. Mech. 128 (2) (2002) 156–163.

[41] J.C. Simo, T.J.R. Hugues, Computational Inelasticity, Springer-Verlag, Berlin, 1998.

[42] J.C. Simo, J.W. Ju, Strain and stress based continuum damage models—I. Formulation, Int. J. Solids Struct. 23 (1987) 821–840.

[43] J.C. Simo, J.W. Ju, Strain and stress based continuum damage models—II. Computational aspects, Int. J. Solids Struct. 23 (1987) 841–869.

[44] B.P. Sinha, K.H. Gerstle, L.G. Tulin, Stress strain relations for concrete under cyclic loading, J. Amer. Concr. Inst. (February) (1964) 195–211.

[45] K.A.S. Susantha, H. Ge, T. Usami, Uniaxial stress–strain relationship of concrete confined by various shaped steel tubes, Engrg. Struct. 23 (2001) 1331–1347.

[46] K. Willam, E.P. Warnke, Constitutive model for the triaxial behaviour of concrete, in: IABSE Seminar on Concrete Structures Subjected to Triaxial Stresses, III-1, 1974.

[47] S. Yazdani, H.L. Schreyer, Combined plasticity and damage mechanics model for plain concrete, J. Engrg. Mech. ASCE 116 (1990) 1435–1450.