/MAT/LAW124 (CDPM2)

Block Format Keyword A concrete material law accounting for plasticity, damage, and strain rate effect.

A Hillerborg regularization method is also available to avoid mesh size dependency in tension. Only a few parameters are needed to use this material law, which makes it user-friendly.

Format

(1)	(3)	(5)	(6)	(7)	(8)	(9)	(10)
/MAT/LAW124/`mat_ID`/`unit_ID` or /MAT/CDPM2/`mat_ID`/`unit_ID`
`mat_title`
$ρ_{i}$ $ρ_{i}$
`E`	$ν$ $ν$		`IDEL`		`IRATE`	`FCUT`
`ECC`	$Q_{h 0}$ $Q_{h 0}$	$f_{t}$ $f_{t}$		$f_{c}$ $f_{c}$		`HP`
`AH`	`BH`	`CH`		`DH`
`AS`	`BS`	`DF`			`DFLAG`	`DTYPE`	`IREG`
`WF`	`WF1`	`FT1`		`EFC`

Definition

Field	Contents	SI Unit Example
`mat_ID`	Material identifier. (Integer, maximum 10 digits)
`unit_ID`	(Optional) Unit identifier. (Integer, maximum 10 digits)
`mat_title`	Material title. (Character, maximum 100 characters)
$ρ_{i}$ $ρ_{i}$	Initial density. (Real)	$[\frac{kg}{m^{3}}]$ $[\frac{kg}{m^{3}}]$
`E`	Young’s modulus. (Real)	$[Pa]$ $[Pa]$
$ν$ $ν$	Poisson’s ratio. (Real)
`IDEL`	Element deletion flag. = 1 (Default) Not activated = 2 Activated (Integer)
`IRATE`	Rate dependency flag. = 1 (Default) Not activated = 2 Activated (Integer)
`FCUT`	Strain rate filtering frequency Default = 10 kHz (Real)	$[Hz]$ $[Hz]$
`ECC`	Eccentricity. Default in Comment 5 (Real)
$Q_{h 0}$ $Q_{h 0}$	Initial hardening. Default = 0.3 (Real)
$f_{t}$ $f_{t}$	Strength limit in tension. (Real)	$[Pa]$ $[Pa]$
$f_{c}$ $f_{c}$	Strength limit in compression. (Real)	$[Pa]$ $[Pa]$
`HP`	Hardening modulus. (Real)
`AH`	Compressive damage 1^st ductility parameter. Default = 0.08 (Real)
`BH`	Compressive damage 2^nd ductility parameter. Default = 0.003 (Real)
`CH`	Compressive damage 3^rd ductility parameter. Default = 2.0 (Real)
`DH`	Compressive damage 4^th ductility parameter. Default = 1.0E-6 (Real)
`AS`	First ductility measure parameter. Default = 15.0 (Real)
`BS`	Second ductility measure parameter. Default = 1.0 (Real)
`DF`	Dilation coefficient. Default = 0.85 (Real)
`DFLAG`	Damage model type. = 1 (Default) Non-symmetric damage = 2 Isotropic = 3 Multiplicative (Integer)
`DTYPE`	Tensile damage shape. = 1 Linear = 2 (Default) Bi-linear = 3 Exponential (Integer)
`IREG`	Regularization flag. = 1 Not activated = 2 (Default) Activated (Integer)
`WF`	Tensile inelastic displacement/strain at failure. Dimension depends on `IREG` parameter value. 3 (Real)	$[m]$ $[m]$
`WF1`	Tensile inelastic displacement/strain at softening slope change (bilinear damage only). Dimension depends on `IREG` parameter value. 3 Default = 0.15*`WF` (Real)	$[m]$ $[m]$
`FT1`	Tensile strength at `WF1`. Default = 0.3*`FT` (Real)	$[Pa]$ $[Pa]$
`EFC`	Compressive inelastic strain close to failure. Default = 1.0E-4 (Real)

▸Example

#RADIOSS STARTER
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/UNIT/1
unit for mat
                  Mg                  mm                   s
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/MAT/LAW124/1/1
Concrete CDPM2
#        Init. dens.
              2.3E-9                  
#                  E                  NU                IDEL               IRATE                FCUT
               28000                0.19                   0                   1                   0
#                ECC                 QH0                  FT                  FC                  HP
                   0                   0                 3.5                33.6                 0.5
#                 AH                  BH                  CH                  DH
                   0                   0                   0                   0
#                 AS                  BS                  DF               DFLAG     DTYPE      Ireg
                   0                   0                   0                   1         1         1
#                 WF                 WF1                 FT1                 EFC     
               0.006                   0                   0              0.0005  
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#enddata
/END
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

Comments

The CDPM2 material law is a user-friendly concrete material law, which considers several phenomena. Only a few parameters are mandatory to easily use this constitutive model. For this law, the elastic behavior is supposed to be isotropic. The plastic behavior is then characterized by the following yield function (Figure 1):(1)
$f_{p} = {⎧ ⎨ ⎩ [1 - q_{h 1} (κ_{p})] {(\frac{¯ ρ}{\sqrt{6} f_{c}} + \frac{{¯ σ}_{v}}{f_{c}})}^{2} + \sqrt{\frac{3}{2}} \frac{¯ ρ}{f_{c}} ⎫ ⎬ ⎭}^{2} + m_{0} q_{h 1}^{2} (κ_{p}) q_{h 2} (κ_{p}) [\frac{¯ ρ}{\sqrt{6} f_{c}} r (cos ¯ θ) + \frac{{¯ σ}_{v}}{f_{c}}] - q_{h 1}^{2} (κ_{p}) q_{h 2}^{2} (κ_{p})$ $f_{p} = {\{[1 - q_{h 1} (κ_{p})] {(\frac{\bar{ρ}}{\sqrt{6} f_{c}} + \frac{{\bar{σ}}_{v}}{f_{c}})}^{2} + \sqrt{\frac{3}{2}} \frac{\bar{ρ}}{f_{c}}\}}^{2} + m_{0} q_{h 1}^{2} (κ_{p}) q_{h 2} (κ_{p}) [\frac{\bar{ρ}}{\sqrt{6} f_{c}} r (\cos \bar{θ}) + \frac{{\bar{σ}}_{v}}{f_{c}}] - q_{h 1}^{2} (κ_{p}) q_{h 2}^{2} (κ_{p})$

Where the Haigh-Westergaard coordinates in stresses space are considered:

${\bar{σ}}_{v} = \frac{tr (σ)}{3}$ , $\bar{ρ} = \sqrt{2 J_{2}} = \sqrt{s : s}$ , $\bar{θ} = \frac{1}{3} \arccos (\frac{3 \sqrt{3}}{2} \frac{J_{3}}{J_{2}^{3 / 2}})$

Where,

$κ_{p}$

Hardening variable

$f_{c}$

Limit strength in compression

$f_{t}$

Limit strength in tension

$m_{0}$

A parameter considering the effect of eccentricity $e_{c c}$

(2)
$m_{0} = \frac{3 (f_{c}^{2} - f_{t}^{2})}{f_{c} f_{t}} \frac{e_{c c}}{e_{c c} + 1}$

Figure 1. CDPM2 model yield function shape (from Grassl)

Then $q_{h 1}$ and $q_{h 2}$ are the two hardening functions defined with (Figure 2):(3)
$q_{h 1} (κ_{p}) = \{\begin{matrix} q_{h 0} + (1 - q_{h 0}) (κ_{p}^{3} - 3 κ_{p}^{2} + 2 κ_{p}) & i f & κ_{p} < 1 \\ 1 & i f & κ_{p} \geq 1 \end{matrix}$
(4)
$q_{h 2} (κ_{p}) = \{\begin{matrix} 1 & i f & κ_{p} < 1 \\ 1 + H_{p} (κ_{p} - 1) & i f & κ_{p} \geq 1 \end{matrix}$

Figure 2. Hardening functions shape (from Grassl)

Here, $q_{h 0}$ is the initial hardening defined so that $0 < q_{h 0} < 1$ . $H_{p}$ is the hardening modulus whose recommended values are 0.01 without strain rate dependency (IRATE = 1), and 0.5 with rate effects (IRATE = 2). The evolution of the internal variable $κ_{p}$ is detailed below.

The William-Warnke function is used to develop the deviatoric section shape between tension and compression (Figure 1).(5)
$r (\cos \bar{θ}) = \frac{4 (1 - e_{c c}^{2}) \cos^{2} \bar{θ} + {(2 e_{c c} - 1)}^{2}}{2 (1 - e_{c c}^{2}) \cos \bar{θ} + (2 e_{c c} - 1) \sqrt{4 (1 - e_{c c}^{2}) \cos^{2} \bar{θ} + 5 e_{c c}^{2} - 4 e_{c c}}}$
The plastic strain evolution is defined by a non-associated plastic potential using the following equation:(6)
$g_{p} = {\{[1 - q_{h 1} (κ_{p})] {(\frac{\bar{ρ}}{\sqrt{6} f_{c}} + \frac{{\bar{σ}}_{v}}{f_{c}})}^{2} + \sqrt{\frac{3}{2}} \frac{\bar{ρ}}{f_{c}}\}}^{2} + q_{h 1}^{2} (κ_{p}) [\frac{m_{0} \bar{ρ}}{\sqrt{6} f_{c}} + \frac{m_{g}}{f_{c}}]$

With

$m_{g} = A_{g} B_{g} f_{c} \exp \frac{{\bar{σ}}_{v} - q_{h 2} \frac{f_{t}}{3}}{B_{g} f_{c}}$

$A_{g} = \frac{3 f_{t} q_{h 2}}{f_{c}} + \frac{m_{0}}{2}$

$B_{g} = \frac{(\frac{q_{h 2}}{3}) (1 + \frac{f_{t}}{f_{c}})}{\ln A_{g} - \ln (2 D_{f} - 1) - \ln (3 q_{h 2} + \frac{m_{0}}{2}) + \ln (D_{f} + 1)}$

Where, $D_{f}$ is the dilation parameter.

This plastic potential is used to compute the evolution of the plastic strain tensor and, thus, the evolution of the internal variable $κ_{p}$ as follow:(7)
${\dot{κ}}_{p} = \frac{‖{\dot{ε}}_{p}‖}{x_{h} ({\bar{σ}}_{v})} (2 \cos^{2} \bar{θ})$

with $x_{h} = \{\begin{matrix} A_{h} - (A_{h} - B_{h}) \exp (\frac{- R_{h} ({\bar{σ}}_{v})}{C_{h}}) & if & R_{h} ({\bar{σ}}_{v}) \geq 0 \\ E_{h} \exp (\frac{R_{h} ({\bar{σ}}_{v})}{F_{h}}) + D_{h} & if & R_{h} ({\bar{σ}}_{v}) < 0 \end{matrix}$

Where, $‖{\dot{ε}}_{p}‖$ is the Euclidian norm of the increment of the plastic strain tensor.
The CDPM2 model considers an unsymmetrical damage evolution between tension and compression. These variables are respectively denoted by $ω_{t}$ and $ω_{c}$ . The damage variables evolution is triggered by a strain criterion defined by:(8)
$ε_{e q} = \frac{ε_{0} m_{0}}{2} (\frac{\bar{ρ}}{\sqrt{6} f_{c}} r (\cos \bar{θ}) + \frac{{\bar{σ}}_{v}}{f_{c}}) + \sqrt{\frac{ε_{0}^{2} m_{0}^{2}}{4} {(\frac{\bar{ρ}}{\sqrt{6} f_{c}} r (\cos \bar{θ}) + \frac{{\bar{σ}}_{v}}{f_{c}})}^{2} + \frac{3 ε_{0}^{2} {\bar{ρ}}^{2}}{2 f_{c}^{2}}} \geq ε_{0} = \frac{f_{t}}{E}$
When this criterion is reached, the damage history variables to the corresponding loading case (tension or compression) are updated:
- In Tension:
  $κ_{d t 2}^{n} = κ_{d t 2}^{n - 1} + \frac{\max (ε_{e q} - κ_{d t}^{n - 1}, 0)}{x_{s}}$ , $κ_{d t 1}^{n} = κ_{d t 1}^{n - 1} + \frac{Δ ε_{p}}{x_{s}}$ , and $κ_{d t}^{n} = ε_{e q}^{n}$
- In Compression:
  $κ_{d c 2}^{n} = κ_{d c 2}^{n - 1} + \frac{α_{c} \max (ε_{e q} - κ_{d c}^{n - 1}, 0)}{x_{s}}$ , $κ_{d c 1}^{n} = κ_{d c 1}^{n - 1} + α_{c} β_{c} \frac{Δ ε_{p}}{x_{s}}$ and $κ_{d c}^{n} = α_{c} ε_{e q}^{n}$
  
  With,
  
  $α_{c} = \sum_{i} \frac{{〈σ_{p i}〉}_{-} ({〈σ_{p i}〉}_{+} + {〈σ_{p i}〉}_{-})}{{‖σ_{p}‖}^{2}}$ , $β_{c} = \frac{f_{t} q_{h 2} (κ_{p}) \sqrt{2 / 3}}{\bar{ρ} \sqrt{1 + 2 D_{f}^{2}}}$ , $x_{s} = 1 + (A_{s} - 1) R_{s}^{B S}$ with $R_{s} = \{\begin{matrix} - \frac{\sqrt{6} {\bar{σ}}_{v}}{\bar{ρ}} & if & {\bar{σ}}_{v} \leq 0 \\ 0 & if & {\bar{σ}}_{v} > 0 \end{matrix}$
The inelastic strains can be obtained from damage history variable using the following equations:

$ε_{i n e l}^{t} = κ_{d t 1} + ω_{t} κ_{d t 2}$ and $ε_{i n e l}^{c} = κ_{d c 1} + ω_{c} κ_{d c 2}$

The damage history variable finally enables to update the corresponding damage variables.
Regarding tensile damage, three different evolution shapes are available depending on the DTYPE parameter value:
- DTYPE = 1: Linear damage where $f_{t}$ is the strength limit at the beginning of damage, and $w_{f}$ is the failure displacement for which the stiffness becomes null (Figure 3).(9)
  $ω_{t} = \frac{E κ_{d t} w_{f} - f_{t} w_{f} + f_{t} κ_{d t 1} h}{E κ_{d t} w_{f} - f_{t} h κ_{d t 2}}$
  
  Figure 3. Uniaxial tension test on a single unit element using linear damage evolution . with $f_{t}$ =3.5 MPa and $w_{f}$ =0.002 mm
- DTYPE = 2: Bilinear damage which is very similar to the linear damage apart from the use of the couple of values $w_{f 1}$ and $f_{t 1}$ which define the coordinates of the points where damage evolution changes its slope (Figure 4).(10)
  $ω_{t} = \{\begin{matrix} \frac{E κ_{d t} w_{f 1} - f_{t} w_{f 1} - (f_{t 1} - f_{t}) κ_{d t 1} h}{E κ_{d t} w_{f 1} + (f_{t 1} - f_{t}) h κ_{d t 2}} & if & \begin{matrix} h ε_{i n e l}^{t} > 0 & and & h ε_{i n e l}^{t} < w_{f 1} \end{matrix} \\ \frac{E κ_{d t} (w_{f} - w_{f 1}) - f_{t 1} (w_{f} - w_{f 1}) + f_{t 1} κ_{d t 1} h - f_{t 1} w_{f 1}}{E κ_{d t} (w_{f} - w_{f 1}) - f_{t 1} h κ_{d t 2}} & if & \begin{matrix} h ε_{i n e l}^{t} > w_{f 1} & and & h ε_{i n e l}^{t} < w_{f} \end{matrix} \end{matrix}$
  
  Figure 4. Uniaxial tension test on a single unit element using bilinear damage evolution . with $f_{t 1}$ =3.5 MPa , $f_{t}$ =1.5 MPa, $w_{f 1}$ =0.00075 mm and $w_{f}$ =0.002 mm
- DTYPE = 3: Exponential damage where the displacement threshold $w_{f}$ corresponds to the meeting point between uniaxial strain axis, and tangent curve to the beginning of stress softening (Figure 5).(11)
  $ω_{t} = 1 - \exp (- \frac{h ε_{i n e l}^{t}}{w_{f}})$
  
  Figure 5. Uniaxial tension test on a single unit element using exponential damage evolution . with $f_{t}$ =3.5 MPa and $w_{f}$ =0.002 mm
In these different equations, $h$ is a parameter that can be used to avoid mesh size dependency. When IREG = 1, no regularization method is used, and $h$ is set to 1. In that case, the critical value $w_{f}$ becomes a critical strain with no dimension. Otherwise, if IREG = 2, the Hillerborg’s regularization method ² (also called Crack Brand method ³) is used, and $h$ equals the initial element size. Then, the critical value $w_{f}$ becomes a critical displacement homogeneous to a displacement. Hillerborg’s regularization method is to ensure that the tensile fracture energy denoted $G_{f}^{t}$ remains constant no matter what element size is used (Figure 6).

Figure 6. Uniaxial tension test with bilinear damage on two different mesh sizes with IREQ = 0 (left); IREQ = 1 (right)

Regarding compression damage, only the exponential evolution shape is available without any regularization method (Figure 7). Mesh size dependency is assumed to be less sensitive.(12)
$ω_{c} = 1 - \exp (- \frac{ε_{i n e l}^{c}}{e_{f c}})$

Figure 7. Uniaxial compression test on a single unit element . with $e_{f c} = 5 × 10^{(- 4)}$
The effect of damage on stress computation will depend on the DFLAG parameter value:
- DFLAG = 1: Non-symmetric softening which considers the crack closure effect when switching from tension to compression, recovering the initial stiffness. On the opposite, a switch from tension to compression re-opens the already existing cracks. (Figure 8).(13)
  $σ = (1 - ω_{t}) σ_{e f f}^{t} + (1 - ω_{c}) σ_{e f f}^{c}$
  
  Where, $σ_{e f f}^{t}$ and $σ_{e f f}^{c}$ are respectively the tension and compression part of the undamaged (effective) stress tensor.
  
  Figure 8. Loading/unloading uniaxial test with non-symmetrical damage softening
- DFLAG = 2: Isotropic softening that considers only the effect of tensile damage in both tension and compression. Then crack closure is not considered. No changes of stiffness are observed when switching from tension to compression or the opposite (Figure 9). Tensile damage is also less likely to evolve in compression.(14)
  $σ = (1 - ω_{t}) σ_{e f f}$
  
  Figure 9. Loading/unloading uniaxial test with isotropic damage softening
- DFLAG = 3: Multiplicative softening where the effect of both tension and compression damage are considered and cumulated on the behavior (Figure 10).(15)
  $σ = (1 - (1 - ω_{t}) (1 - ω_{c})) σ_{e f f}$
  
  Figure 10. Loading/unloading uniaxial test with multiplicative damage softening
The last phenomena considered by CDPM2 model is the strain rate dependency. At a high strain rate, the concrete is more likely to have a larger tension or compression strength limit. This is introduced by the following equations:
$f_{t}^{r a t e} = α_{r a t e} f_{t}$ , $f_{c}^{r a t e} = α_{r a t e} f_{c}$ and $f_{t 1}^{r a t e} = α_{r a t e} f_{t 1}$

The dynamic increase factor (DIF) $α_{r a t e}$ is computed with:(16)
$α_{r a t e} = (1 - α_{c}) α_{r a t e}^{t} + α_{c} α_{r a t e}^{c}$

Where, $α_{c}$ is the compression factor defined above in damage history variables equations in Comment 3.

There is then a different strain rate dependency between tension and compression as concrete is more sensitive to strain rate effect in tension than in compression. The two dynamic increased factor for both tension and compression are computed with:

$α_{r a t e}^{t} = \{\begin{matrix} 1 & for & {\dot{ε}}_{\max} \leq 30 \times 10^{- 6} s^{- 1} \\ {(\frac{\dot{ε}}{{\dot{ε}}_{t 0}})}^{δ_{s}} & for & 30 \times 10^{- 6} s^{- 1} < {\dot{ε}}_{\max} \leq 1 s^{- 1} \\ β_{s} {(\frac{\dot{ε}}{{\dot{ε}}_{t 0}})}^{\frac{1}{3}} & for & 1 s^{- 1} \leq {\dot{ε}}_{\max} \end{matrix}$ with $δ_{s} = \frac{1}{1 + 8 \frac{f_{c}}{f_{c 0}}}$ , $\log β_{s} = 6 δ_{s} - 2$

$α_{r a t e}^{c} = \{\begin{matrix} 1 & for & {\dot{ε}}_{\max} \leq 30 \times 10^{- 6} s^{- 1} \\ {(\frac{\dot{ε}}{{\dot{ε}}_{c 0}})}^{1.026 α_{s}} & for & 30 \times 10^{- 6} s^{- 1} < {\dot{ε}}_{\max} \leq 30 s^{- 1} \\ γ_{s} {(\frac{\dot{ε}}{{\dot{ε}}_{c 0}})}^{\frac{1}{3}} & for & {30s}^{- 1} \leq {\dot{ε}}_{\max} \end{matrix}$ with $α_{s} = \frac{1}{5 + 9 \frac{f_{c}}{f_{c 0}}}$ , $\log γ_{s} = 6.56 α_{s} - 2$

The equivalent deviatoric strain rate $\dot{ε}$ is used to compute the DIF in the equations above.

No parameters need to be identified for strain rate effect. You only need to set the flag IRATE to 2. Figure 11 shows the expected tendency of strain rate effect on the CDPM2 behavior. By increasing the strength limits in tension/compression, the dissipated energy during failure is also affected which is often observed experimentally.
Figure 11. Uniaxial tests in tension/compression with strain rate dependency (IRATE = 1)
Default value of eccentricity can be obtained with:
$ε_{i} = \frac{f_{t} ({(1.16 f_{c})}^{2} - f_{c}^{2})}{1.16 f_{c} (f_{c}^{2} - f_{t}^{2})}$ and $e_{c c} = \frac{1 + ε_{i}}{2 - ε_{i}}$
Global damage variable can be output using /ANIM/BRICK/DAMG or /HED/SOLID/DAMG.

Peter Grassl, Dimitrios Xenos, Ulrika Nyström, Rasmus Rempling, Kent Gylltoft, CDPM2: A damage-plasticity approach to modelling the failure of concrete, International Journal of Solids and Structures, Volume 50, Issue 24, 2013, Pages 3805-3816, ISSN 0020-7683

A. Hillerborg, M. Modéer, P.-E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, Cement and Concrete Research, Volume 6, Issue 6, 1976, Pages 773-781, ISSN 0008-8846

Bažant, Z.P., Oh, B.H. Crack band theory for fracture of concrete. Mat. Constr. 16, 155–177 (1983)