Damage-plastic model for concrete - ConcreteDPM

This model, developed by Grassl and Jirásek for failure of concrete under general triaxial stress, is described in detail in [8]. It belongs to the class of damage-plastic models with yield condition formulated in terms of the effective stress $\bar{\mbox{\boldmath $\sigma$}}=\mbox{\boldmath $D$}_e:(\mbox{\boldmath $\varepsilon$}-\mbox{\boldmath $\varepsilon$}_p)$. The stress-strain law is postulated in the form

$$\mbox{\boldmath$\sigma$} = (1-\omega)\bar{\mbox{\boldmath$\sigma$}}=(1-\omega)\mbox{\boldmath$D$}_e:(\mbox{\boldmath$\varepsilon$}-\mbox{\boldmath$\varepsilon$}_p)$$ (141)

where $D$$_e$ is the elastic stiffness tensor and $\omega$ is a scalar damage parameter. The plastic part of the model consists of a three-invariant yield condition, nonassociated flow rule and pressure-dependent hardening law. For simplicity, damage is assumed to be isotropic. In contrast to pure damage models with damage driven by the total strain, here the damage is linked to the evolution of plastic strain.

The yield surface is described in terms of the cylindrical coordinates in the principal effective stress space (Haigh-Westergaard coordinates), which are the volumetric effective stress $\bar{\sigma}_{\rm V} = {I_1(\bar{\mbox{\boldmath $\sigma$}})}/{3}$, the norm of the deviatoric effective stress $\bar{\rho} = \sqrt{2 J_2(\bar{\mbox{\boldmath $\sigma$}})}$, and the Lode angle $\theta$ defined by the relation

$$\cos 3 \theta = \frac{3 \sqrt{3}}{2}\frac{J_3}{J_2^{3/2}}$$ (142)

where $J_2$ and $J_3$ are the second and third deviatoric invariants. The yield function

$$f_{\rm p}(\bar\sigma_{\rm V},\bar{\rho},\bar{\theta};\kappa_{\rm p}) = \left(\left[1-q_{\rm {h}}(\kappa_{\rm p})\right]\left( \frac{\bar... ...{f}_c} \right)^2 + \sqrt{\frac{3}{2}} \frac {\bar{\rho}}{\bar{f}_c} \right)^2 +$$

$$+m_0 q_{\rm {h}}^2(\kappa_{\rm p}) \left(\frac{\bar{\rho}r(\bar{\... ... \frac{\bar{\sigma}_{\rm V}}{\bar{f}_c} \right) - q_{\rm {h}}^2(\kappa_{\rm p})$$ (143)

depends on the effective stress (which enters in the form of cylindrical coordinates) and on the hardening variable $\kappa_{\rm p}$ (which enters through a dimensionless variable $q_{\rm h}$). Parameter $\bar{f}_c$ is the uniaxial compressive strength. Note that, under uniaxial compression characterized by axial stress $\bar{\sigma}<0$, we have $\bar{\sigma}_{\rm V}=\bar{\sigma}/3$, $\bar{\rho}=-\sqrt{2/3}\,\bar{\sigma}$ and $\bar{\theta}=60^o$. The yield function then reduces to $f_{\rm p}=(\bar{\sigma}/\bar{f}_c)^2-q_{\rm {h}}^2$. This means that function $q_{\rm {h}}$ describes the evolution of the uniaxial compressive yield stress normalized by its maximum value, $\bar{f}_c$.

The evolution of the yield surface during hardening is presented in Fig. 8. The parabolic shape of the meridians (Fig. 8a) is controlled by the hardening variable $q_{\rm h}$ and the friction parameter $m_0$. The initial yield surface is closed, which allows modeling of compaction under highly confined compression. The initial and intermediate yield surfaces have two vertices on the hydrostatic axis but the ultimate yield surface has only one vertex on the tensile part of the hydrostatic axis and opens up along the compressive part of the hydrostatic axis. The deviatoric sections evolve as shown in Fig. 8b, and their final shape at full hardening is a rounded triangle at low confinement and almost circular at high confinement. The shape of the deviatoric section is controlled by the Willam-Warnke function

$$r(\theta)=\frac{4(1-e^2)\cos^2\theta+(2e-1)^2}{2(1-e^2)\cos\theta+(2e-1)\sqrt{4(1-e^2)\cos^2\theta+5e^2-4e}}$$ (144)

The eccentricity parameter $e$ that appears in this function, as well as the friction parameter $m_0$, are calibrated from the values of uniaxial and equibiaxial compressive strengths and uniaxial tensile strength.

Figure: Evolution of the yield surface during hardening: a) meridional section, b) deviatoric section for a constant volumetric effective stress of $\bar{\sigma}_{\rm V} = -\bar{f}_c/3$

The maximum size of the elastic domain is attained when the variable $q_{\rm h}$ is equal to one (which is its maximum value, as follows from the hardening law, to be specified in (149)). The yield surface is then described by the equation

$$f_{\rm p}\left(\bar{\sigma}_{\rm V},\bar{\rho},\bar{\theta};1\rig... ...r{f}_c}r(\bar{\theta}) + \frac{\bar{\sigma}_{\rm V}}{\bar{f}_c} \right) - 1 = 0$$ (145)

The flow rule

$$\dot{\mbox{\boldmath$\mbox{\boldmath$\varepsilon$}$}}_{\rm p} = \dot{\lambda} \frac{\partial g_{\rm p}}{\partial\bar{\sigma}}$$ (146)

is non-associated, which means that the yield function $f_{\rm p}$ and the plastic potential

$$g_{\rm p}(\bar{\sigma}_{\rm V},\bar{\rho};\kappa_{\rm p}) = \left(\left[1-q_{\rm {h}}(\kappa_{\rm p})\right]\left( \frac{\bar... ...{f}_c} \right)^2 + \sqrt{\frac{3}{2}} \frac {\bar{\rho}}{\bar{f}_c} \right)^2 +$$

$$+q_{\rm {h}}^2(\kappa_{\rm p}) \left( \frac{m_0 \bar{\rho}}{\sqrt{6}\bar{f}_c} + \frac{m_{\rm g}(\bar{\sigma}_{\rm V})}{\bar{f}_c} \right)$$ (147)

do not coincide and, therefore, the direction of the plastic flow $\partial g_{\rm p}/\partial\bar{\mbox{\boldmath $\sigma$}}$ is not normal to the yield surface. The ratio of the volumetric and the deviatoric parts of the flow direction is controled by function $m_{\rm g}$, which depends on the volumetric stress and is defined as

$$m_{\rm g}(\bar{\sigma}_{\rm V})=A_{\rm g}B_{\rm g}\bar{f}_c\exp{\frac{\bar{\sigma}_{\rm V} -\bar{f}_t/3}{B_{\rm g}\bar{f}_c}}$$ (148)

where $A_{\rm g}$ and $B_{\rm g}$ are model parameters that are determined from certain assumptions on the plastic flow in uniaxial tension and compression.

The dimensionless variable $q_{\rm h}$ that appears in the yield function (143) is a function of the hardening variable $\kappa_{\rm p}$. It controls the size and shape of the yield surface and, thereby, of the elastic domain. The hardening law is given by

$$q_{\rm h}(\kappa_{\rm p}) = \left\{ \begin{array}{ll} {q_{\rm h}}... ...kappa_{\rm p} < 1$} \\ 1 & \mbox{if $\kappa_{\rm p} \ge 1$} \end{array}\right.$$ (149)

The initial inclination of the hardening curve (at $\kappa_{\rm p}=0$) is positive and finite, and the inclination at peak (i.e., at $\kappa_{\rm p} = 1$) is zero.

The evolution law for the hardening variable,²

$$\dot{\kappa}_{\rm {p}} = \frac {\Vert \dot{\mbox{\boldmath$\varep... ... {p}}\Vert}{x_{\rm {h}}\left(\bar{\sigma}_{\rm V} \right)}(2\cos\bar{\theta})^2$$ (150)

sets the rate of the hardening variable equal to the norm of the plastic strain rate scaled by a hardening ductility measure

$$x_{\rm {h}}\left(\bar{\sigma}_{\rm V}\right) = \left\{ \begin{arr... ...+D_{\rm h} & \mbox{if $R_{\rm h}(\bar{\sigma}_{\rm V}) < 0$} \end{array}\right.$$ (151)

The dependence of the scaling factor $x_{\rm {h}}$ on the volumetric effective stress $\bar{\sigma}_{\rm V}$ is constructed such that the model response is more ductile under compression. The variable

$$R_{\rm h}(\bar{\sigma}_{\rm V}) = -\frac{\bar{\sigma}_{\rm V}}{\bar{f}_c}-\frac{1}{3}$$ (152)

is a linear function of the volumetric effective stress. Model parameters $A_{\rm h}, B_{\rm h}, C_{\rm h}$ and $D_{\rm h}$ are calibrated from the values of strain at peak stress under uniaxial tension, uniaxial compression and triaxial compression, whereas the parameters

$$E_{\rm h} = B_{\rm h} - D_{\rm h}$$ (153)

$$F_{\rm h} = \frac{\left(B_{\rm h} -D_{\rm h}\right)C_{\rm h}}{B_{\rm h} - A_{\rm h}}$$ (154)

are determined from the conditions of a smooth transition between the two parts of equation (151) at $R_{\rm h} = 0$.

For the present model, the evolution of damage starts after full saturation of plastic hardening, i.e., at $\kappa_{\rm p} = 1$. This greatly facilitates calibration of model parameters, because the strength envelope is fully controled by the plastic part of the model and damage affects only the softening behavior. In contrast to pure damage models, damage is assumed to be driven by the plastic strain, more specifically by its volumetric part, which is closely related to cracking. To slow down the evolution of damage under compressive stress states, the damage-driving variable $\kappa_{\rm d}$ is not set equal to the volumetric plastic strain, but it is defined incrementally by the rate equation

$$\dot{\kappa}_{\rm d} = \left\{ \begin{array}{ll} 0 & \mbox{if $\k... ...bar{\sigma}_{\rm V}\right) & \mbox{if $\kappa_{\rm p}\ge 1$} \end{array}\right.$$ (155)

where

$$x_{\rm s}(\bar{\sigma}_{\rm V}) = \left\{\begin{array}{ll} 1 + A_... ...m V})} & \mbox{if $R_{\rm s}(\bar{\sigma}_{\rm V}) \geq 1$} \end{array}\right .$$ (156)

is a softening ductility measure. Parameter $A_{\rm s}$ is determined from the softening response in uniaxial compression. The dimensionless variable $R_{\rm s}=\dot{\varepsilon}_{\rm pV}^{-}/\dot{\varepsilon}_{\rm {pV}}$ is defined as the ratio between the “negative” volumetric plastic strain rate

$$\dot{\varepsilon}_{\rm pV}^{-} = \displaystyle{\sum_{I = 1}^3 \langle-\dot{\varepsilon}_{\rm {p}I}\rangle}$$ (157)

and the total volumetric plastic strain rate $\dot{\varepsilon}_{\rm {pV}}$. This ratio depends only on the flow direction $\partial g_{\rm p}/\partial\bar{\mbox{\boldmath $\sigma$}}$, and thus $R_{\rm s}$ can be shown to be a unique function of the volumetric effective stress. In (157), $\dot{\varepsilon}_{\rm {p}I}$ are the principal components of the rate of plastic strains and $\left<\cdot\right>$ denotes the McAuley brackets (positive-part operator). For uniaxial tension, for instance, all three principal plastic strain rates are nonnegative, and so $\dot{\varepsilon}_{\rm pV}^{-} =0$, $R_{\rm s}=0$ and $x_{\rm s}=1$. This means that under uniaxial tensile loading we have $\kappa_{\rm d}=\kappa_{\rm p}-1$. On the other hand, under compressive stress states the negative principal plastic strain rates lead to a ductility measure $x_{\rm s}$ greater than one and the evolution of damage is slowed down. It should be emphasized that the flow rule for this specific model is constructed such that the volumetric part of plastic strain rate at the ultimate yield surface cannot be negative.

Table 40: Damage-plastic model for concrete - summary.

Description	Damage-plastic model for concrete
Record Format	ConcreteDPM d(rn) # E(rn) # n(rn) # tAlpha(rn) # ft(rn) # fc(rn) # wf(rn) # Gf(rn) # ecc(rn) # kinit(rn) # Ahard(rn) # Bhard(rn) # Chard(rn) # Dhard(rn) # Asoft(rn) # helem(rn) # href(rn) # dilation(rn) # yieldtol(rn) # newtoniter(in) #
Parameters	- d material density
	- E Young modulus
	- n Poisson ratio
	- tAlpha thermal dilatation coefficient
	- ft uniaxial tensile strength $f_t$
	- fc uniaxial compressive strength
	- wf parameter $w_f$ that controls the slope of the softening branch (serves for the evaluation of $\varepsilon$$_f=w_f/h$ to be used in (158))
	- Gf fracture energy, can be specified instead of wf, it is converted to $w_f=G_f/f_t$
	- ecc eccentricity parameter $e$ from (144), optional, default value 0.525
	- kinit parameter ${q_{\rm h}}_0$ from (149), optional, default value 0.1
	- Ahard parameter $A_{\rm h}$ from (151), optional, default value 0.08
	- Bhard parameter $B_{\rm h}$ from (151), optional, default value 0.003
	- Chard parameter $C_{\rm h}$ from (151), optional, default value 2
	- Dhard parameter $D_{\rm h}$ from (151), optional, default value $10^{-6}$
	- Asoft parameter $A_{\rm s}$ from (156), optional, default value 15
	- helem element size $h$, optional (if not specified, the actual element size is used)
	- href reference element size $h_{ref}$, optional (if not specified, the standard adjustment of the damage law is used)
	- dilation dilation factor (ratio between lateral and axial plastic strain rates in the softening regime under uniaxial compression), optional, default value -0.85
	- yieldtol tolerance for the implicit stress return algorithm, optional, default value $10^{-10}$
	- newtoniter maximum number of iterations in the implicit stress return algorithm, optional, default value 100
Supported modes	3dMat

The relation between the damage variable $\omega$ and the internal variable $\kappa_{\rm d}$ (maximum level of equivalent strain) is assumed to have the exponential form

$$\omega = 1-\exp\left(-\kappa_{\rm d}/\varepsilon_f\right)$$ (158)

where $\varepsilon _f$ is a parameter that controls the slope of the softening curve. In fact, equation (158) is used by the nonlocal version of the damage-plastic model, with $\kappa_{\rm d}$ replaced by its weighted spatial average (not yet available in the public version of OOFEM). For the local model, it is necessary to adjust softening according to the element size, otherwise the results would suffer by pathological mesh sensitivity. It is assumed that localization takes place at the peak of the stress-strain diagram, i.e., at the onset of damage. After that, the strain is decomposed into the distributed part, which corresponds to unloading from peak, and the localized part, which is added if the material is softening. The localized part of strain is transformed into an equivalent crack opening, $w$, which is under uniaxial tension linked to the stress by the exponential law

$$\sigma = (1-\omega) f_t = f_t \exp\left(-w/w_f\right)$$ (159)

Here, $f_t$ is the uniaxial tensile strength and $w_f$ is the characteristic crack opening, playing a similar role to $\varepsilon _f$. Under uniaxial tension, the localized strain can be expressed as the sum of the post-peak plastic strain (equal to variable $\kappa_{\rm d}$) and the unloaded part of elastic strain (equal to $\omega f_t/E$). Denoting the effective element size as $h$, we can write

$$w = h(\kappa_{\rm d}+\omega f_t/E)$$ (160)

and substituting this into (159), we obtain a nonlinear equation

$$1-\omega = \exp\left(-\frac{h}{w_f}\left(\kappa_{\rm d}+\omega f_t/E\right)\right)$$ (161)

from which the damage variable $\omega$ corresponding to the given internal variable $\kappa_{\rm d}$ can be computed by Newton iteration. The effective element size $h$ is obtained by projecting the element onto the direction of the maximum principal strain at the onset of cracking, and afterwards it is held fixed. The evaluation of $\omega$ from $\kappa_{\rm d}$ is no longer explicit, but the resulting load-displacement curve of a bar under uniaxial tension is totally independent of the mesh size. A simpler approach would be to use (158) with $\varepsilon_f = w_f/h$, but then the scaling would not be perfect and the shape of the load-displacement curve (and also the dissipated energy) would slightly depend on the mesh size. With the present approach, the energy per unit sectional area dissipated under uniaxial tension is exactly $G_f = w_f f_t$. The input parameter controling the damage law can be either the characteristic crack opening $w_f$, or the fracture energy $G_f$. If both are specified, $w_f$ is used and $G_f$ is ignored. If only $G_f$ is specified, $w_f$ is set to $G_f/f_t$.

The onset of damage corresponds to the peak of the stress-strain diagram under proportional loading, when the ratios of the stress components are fixed. This is the case e.g. for uniaxial tension, uniaxial compression, or shear under free expansion of the material (with zero normal stresses). However, for shear under confinement the shear stress can rise even after the onset of damage, due to increasing hydrostatic pressure, which increases the mobilized friction. It has been observed that the standard approach leads to strong sensitivity of the peak shear stress to the element size. To reduce this pathological effect, a modified approach has been implemented. The second-order work (product of stress increment and strain increment) is checked after each step and the element-size dependent adjustment of the damage law is applied only after the second-order work becomes negative. Up to this stage, the damage law corresponds to a fixed reference element size, which is independent of the actual size of the element. This size is set by the optional parameter href. If this parameter is not specified, the standard approach is used. For testing purposes, one can also specify the actual element size, helem, as a “material property”. If this parameter is not specified, the element size is computed for each element separately and represents its actual size.

The damage-plastic model contains 15 parameters, but only 6 of them need to be actually calibrated for different concrete types, namely Young's modulus $E$, Poisson's ratio $\nu$, tensile strength $f_{\rm t}$, compressive strength $f_{\rm c}$, parameter $w_f$ (or fracture energy $G_f$), and parameter $A_{\rm s}$ in the ductility measure (156) of the damage model. The remaining parameters can be set to their default values specified in [8].

The model parameters are summarized in Tab. 40. Note that it is possible to specify the “size” of finite element, $h$, which (if specified) replaces the actual element size in (161). The usual approach is to consider $h$ as the actual element size (evaluated automatically by OOFEM), in which case the optional parameter $h$ is missing (or is set to 0., which has the same effect in the code). However, for various studies of mesh sensitivity it is useful to have the option of specifying $h$ as an input “material” value.

If the element is too large, it may become too brittle and local snap-back occurs in the stress-strain diagram, which is not acceptable. In such a case, an error message is issued and the program execution is terminated. The maximum admissible element size

$$h_{\max} = \frac{EG_f}{f_t^2} = \frac{Ew_f}{f_t}$$ (162)

happens to be equal to Hillerborg's characteristic material length. For typical concretes it is in the order of a few hundred mm. If the condition $h<h_{\max}$ is violated, the mesh needs to be refined. Note that the effective element size $h$ is obtained by projecting the element. For instance, if the element is a cube of edge length 100 mm, its effective size in the direction of the body diagonal can be 173 mm.