Concrete and Rock Materials

In Radioss these materials can be used to represent rock or concrete materials.

These materials use a Drücker–Prager yield criterion¹, which is a pressure-dependent model for determining whether a material has failed or undergone plastic yielding.

Concrete Material (/MAT/LAW10 and /MAT/LAW21)

Drücker-Prager Yield Criteria

The material has failed or undergone plastic yielding is determined by pressure using:(1)

F = \underset{J_{2} \begin{matrix}  \end{matrix} p a r t}{\underset{︸}{J_{2}}} - \underset{I_{1} \begin{matrix}  \end{matrix} p a r t}{\underset{︸}{(A_{0} + A_{1} P + A_{2} P^{2})}}

Where,

$J_{2}$: Second stress invariant (von Mises stress) of the deviatoric part of the stress and $P = - \frac{I_{1}}{3}$ .
$I_{1}$: First stress invariant (hydrostatic pressure).
$I_{1} = σ_{1} + σ_{2} + σ_{3} = - 3 P$
$J_{2} = \frac{1}{6} [{(σ_{1} - σ_{2})}^{2} + {(σ_{2} - σ_{3})}^{2} + {(σ_{3} - σ_{1})}^{2}] = \frac{1}{3} σ_{V M}^{2}$: In a uniaxial test.

A polynomial equation is used to describe the pressure

A_{0} + A_{1} P + A_{2} P^{2}

at the Drücker-Prager yield surface of the material: (2)

σ_{V M} = \sqrt{3 (A_{0} + A_{1} P + A_{2} P^{2})}

The constants of the polynomial

A_{0}, A_{1}, A_{2}

are determined by:

If $F < 0$ , $J_{2} < A_{0} + A_{1} P + A_{2} P^{2}$ the material is under yield surface and is in the elastic region.
If $F = 0$ , $J_{2} = A_{0} + A_{1} P + A_{2} P^{2}$ and the material is at the yield surface.
If $F > 0$ , $J_{2} > A_{0} + A_{1} P + A_{2} P^{2}$ and the material is past the yield surface and has failed.
If $A_{1} = A_{2} = 0$ , $σ_{V M} = \sqrt{3 J_{2}} = \sqrt{3 A_{0}}$ , which is the von Mises criterion.

Figure 2.

Pressure Computation

In LAW10, a polynomial equation with input parameters

C_{0}, C_{1}, C_{2}, C_{3}

is used to describe the pressure. The pressure can be plotted as a function of volumetric strain.(3)

μ = \frac{ρ}{ρ_{0}} - 1

If $P_{e x t} = 0$ , the pressure is $P = Δ P$ and the pressure limit is $P_{\min} = Δ P_{\min}$ .

Figure 3. Pressure curve without external pressure
If $P_{e x t} \neq 0$ , the pressure is shifted by $P_{e x t}$ , then $P = P_{e x t} + Δ P$ and the pressure limit is $P_{\min} = P_{e x t} + Δ P_{\min}$ .

Figure 4. Pressure curve with external pressure

Here, (4)

Δ P = {\begin{matrix} \max {Δ P_{\min}, C_{0} + C_{1} μ + C_{2} μ^{2} + C_{3} μ^{3}} \\ \max {Δ P_{\min}, C_{0} + C_{1} μ} \end{matrix} \begin{matrix} i f & μ \geq 0 \begin{matrix}  \end{matrix} c o m p r e s s i o n \\ i f & μ < 0 \begin{matrix}  \end{matrix} t r a c t i o n \end{matrix}

In traction or tension the pressure is linear and limited by $Δ P_{\min}$ .
In compression the pressure is nonlinear also limited by $Δ P_{\min}$ .

The only difference between the material laws is that in LAW10 the material constants $C_{0}, C_{1}, C_{2}, C_{3}$ are used to describe the pressure versus volumetric strain ( $P - μ$ curve). In LAW21 you can describe this curve via function input fct_ID_f.

Load and Unload

In LAW10 and LAW21 different loading and unloading paths of the

P - μ

curve can be considered by using the parameters

μ_{\max}

and B.

In Tension ( $μ < 0$ )
- For LAW10, linear loading and unloading with $P = C_{1} μ$ (Figure 3).
- For LAW21, loading is defined using the input function fct_ID_f and linear unloading with $P = K_{t} μ$ .
In Compression ( $μ > 0$ ), for both LAW10 and LAW21:
- If neither B and $μ_{\max}$ are defined, the loading and unloading path are identical.
  
  Figure 5. Identical loading and unloading for LAW10 and LAW21
- If either B or $μ_{\max}$ is defined:
  1. If only B is defined, $μ_{\max}$ is the volumetric strain where the tangent of $P - μ$ curve is equal to B with $B = {\frac{d P}{d μ} |}_{μ_{\max}}$ .
  2. If only $μ_{\max}$ is defined, then B is the tangent of $P - μ$ curve at $μ_{\max}$ . The loading and unloading in compression is:
    - If $μ > μ_{\max}$ , loading and unloading path are identical.
    - If $μ < μ_{\max}$ , loading and unloading path are different, it is linear unloading with slope B.
      
      Figure 6. Different loading and unloading treatment for LAW10 and LAW21

Concrete Material (/MAT/LAW24)

LAW24 uses a Drücker-Prager criteria with or without a cap in yield to model a reinforced concrete material. This material law assumes that the two failure mechanisms of the concrete material are tensile cracking and compressive crushing.

Concrete Tensile Behavior

In LAW24, the options H_t, D_sup, and

ε_{\max}

can be used to describe tensile cracking and failure in tension.

In the initial very small elastic phase, the material has an elastic modulus E_c.

Once tensile strength, f_t is reached, the concrete starts to soften with the slope H_t. The maximum damage factor, D_sup, is significant because it enables the modeling of residual stiffness during and after a crack.

The residual stiffness is computed as:(5)

E = (1 - D_{\sup}) \cdot E_{c}

When there is crack closure, the concrete becomes elastic again, and the damage factor (for each direction) is conserved.

The bearing capacity of concrete in tensile is much lower than in compression. It is normally considered elastic when in tension.

It is recommended to choose a D_sup value close to 1 (default is 0.99999) in order to minimize the current stiffness at the end of the damage and consequently avoid residual stress in tension, which can become very high if the element is highly deformed due to tension. This will happen if the force causing the damage remains.

It is possible to adjust the D_sup (and H_t) in order to simulate and fit the behavior of concrete reinforced by fibers. The concrete material fails once it reaches the total failure strain $ε_{\max}$ .

Concrete Yield Surface in Compression

For concrete, the yield surface is the beginning of the plastic hardening zone which is between the failure surface, $r_{f}$ , and the yield surface.

The yield surface is assumed to be the same as the failure surface in the tension zone. In compression, the yield surface is a scaled down failure surface using the factor

k (σ_{m}, k_{0})

. The yield in LAW24 for concrete is:(6)

f = \underset{\begin{array}{l} J_{2} \\ p a r t \end{array}}{\underset{︸}{r}} - \underset{I_{1} \begin{matrix}  \end{matrix} p a r t}{\underset{︸}{k (σ_{m}, k_{0}) \cdot r_{f}}} = 0

For I_cap =0 or 1 (without a cap in yield) the yield curve is:

Figure 9. Drücker-Prager Criteria without a Cap in Yield
For I_cap =2 (with cap in yield) the yield is:

Figure 10. Drücker-Prager Criteria with Cap in Yield

$r < k (σ_{m}, k_{0}) \cdot r_{f}$ (green area in Figure 10)

The material is under yield in the elastic phase.

$r \geq r_{f}$ (red area in Figure 10)

The material has failed.

$k (σ_{m}, k_{0}) \cdot r_{f} < r < r_{f}$ (yellow area in Figure 10)

The material is above yield and below the failure surface which is the plastic hardening phase.

The input parameter $ρ_{t}$ is the hydrostatic failure pressure in a uniaxial tension test and $ρ_{c}$ is the hydrostatic pressure by failure in a uniaxial compression test.

The scale factor

k (σ_{m}, k_{0})

is a function of mean stress

σ_{m}

and can be described as:

When $σ_{m} \geq ρ_{t}$ (in tension) the scale factor $k (σ_{m}, k_{0}) = 1$ . In this case, the yield surface equals the failure surface, $r = r_{f}$ .

Figure 11. $k$ Function in the Tension Zone
In the tension-compression region, $ρ_{t} > σ_{m} \geq ρ_{c}$ , then
$k (σ_{m}, k_{0}) = 1 + \frac{(1 - k_{0}) \cdot [ρ_{t} (2 ρ_{c} - ρ_{t}) - 2 ρ_{c} σ_{m} + σ_{m}^{2}]}{{(ρ_{c} - ρ_{t})}^{2}}$ with $k_{y} \leq k_{0} \leq 1$

Figure 12. $k$ Function in the Compression-Tension Mixed Zone
The rest of the curve depends on the I_cap option and the different scale factors $k (σ_{m}, k_{0})$ used.
- For I_cap =0 or 1 and $σ_{m} < ρ_{c}$ (in compression), then $k (σ_{m}, k_{0}) = k_{y}$
  
  Figure 13. $k$ Function in the Compression Zone
- For I_cap =2 (with cap in yield) and $ρ_{c} < σ_{m} < f_{k}$ (in compression), then $k (σ_{m}, k_{0}) = k_{y}$
  
  Figure 14. $k$ Function in Drücker-Prager Criteria without a Cap
- In $f_{k} < σ_{m} < f_{0}$ (in cap zone)
  $k (σ_{m}, k_{0}) = k_{0} [1 - {(\frac{σ_{m} - f_{k}}{f_{0} - f_{k}})}^{2}]$ with $0 \leq k_{0} \leq k_{y}$
  
  Figure 15. $k$ Function in Drücker-Prager Criteria with a Cap

The material constant $k_{y}$ should be $0 \leq k_{y} \leq 1$ . A higher value of $k_{y}$ results in a higher yield surface.

For example, if I_cap =2 (yield with cap), the difference of yield surface between

k_{y} = 0.8

and

k_{y} = 0.6

(Figure 16). The default value of

k_{y}

in LAW24 is 0.5.

Concrete Plastic Flow Rule in Compression

A non-associated plastic flow rule is used in LAW24. The plastic flow rule is:(7)

g = α I_{1} + \sqrt{J_{2}}

Where,

$α$: Plastic dilatancy.
$α = \frac{\partial g}{\partial I_{1}}$: Governs the volumetric plastic flow.
$I_{1}$: First stress invariant (hydrostatic pressure).

Experimentally,

α

is a linear function of

k_{0}

:(8)

α = \frac{(1 - k_{0}) α_{y} + (k_{0} - K_{y}) α_{f}}{1 - K_{y}}

If $k_{0} = K_{y}$: then, $α = α_{y}$ which means the material is in yield.
If $k_{0} < K_{y}$: then, $α$ becomes negative is the cap region.
If $k_{0} = 1$: then, $α = α_{f}$ which means the material has failed.

The values of $α_{y}, α_{f}$ are used to describe the material beyond yield, but before failure. It is recommended to use -0.2 and -0.1 for $α_{y}, α_{f}$ in LAW24. If very small values of $α_{y}, α_{f}$ are used, there is no volumetric plasticity (no cap region).

Concrete Crushing Failure in Compression

Failure surface is given by: (9)

f = r - r_{f} (σ_{m}, θ) = 0

Where,

r = \sqrt{2 J_{2}} / f_{c}

σ_{m} = I_{1} / 3 f_{c}

and

θ

is Lode angle, such as:(10)

\cos 3 θ = \frac{J_{3}}{2} {(\frac{3}{J_{2}})}^{3 / 2}

An Ottosen surface is built to design this surface using: (11)

r_{f} (σ_{m}, θ) = \frac{1}{a} (- b + \sqrt{b^{2} - a (σ_{m} - c)})

Where,

a

b_{c}

b_{t}

and

c

are 4 values which shape the surface and (12)

b (b_{c}, b_{t}, θ) = \frac{1}{2} [b_{c} (1 - \cos 3 θ) + b_{t} (1 + \cos 3 θ)]

For concrete, the compression failure curve

r_{f}

can be defined with a strength of:

$f_{t}$: Uniaxial tension (triaxiality is 1/3)
$f_{c}$: Uniaxial compression (triaxiality is -1/3)
$f_{b}$: Biaxial compression (triaxiality is -2/3)
$f_{2}$: Confined compression strength (tri-axial test)
$s_{0}$: Under confined pressure

The best way to fully determine the 3D failure envelope is to get experimental data for all of these values,

f_{c}, f_{t}, f_{b}, f_{2}, s_{0}

which are schematically illustrated in Figure 18.

Table 1. Input from the 4 Experimental Tests
Load Type	Surface Point	Default Input	Criteria $r$	Pressure $σ_{m}$	Lode Angle $θ$
Compression	$(f_{c}, 0, 0)$	Mandatory	$r = \sqrt{2 / 3}$	$σ_{m} = - 1 / 3$	$\cos θ = - 1$
Direct Tensile	$(f_{t}, 0, 0)$	$f_{t} =0 .1 f_{c}$	$r = \sqrt{2 / 3} (\frac{f_{t}}{f_{c}})$	$σ_{m} = 1 / 3 (\frac{f_{t}}{f_{c}})$	$\cos θ = 1$
Biaxial Compression	$(f_{2}, s_{0}, s_{0})$	If `I`_cap = 1 $f_{2} =4 .0 f_{c}$ If `I`_cap = 2 $f_{2} =7 .0 f_{c}$ $s_{0} =1 .25 f_{c}$	$r = \sqrt{2 / 3} (\frac{f_{t}}{f_{c}}) σ_{m}$	$σ_{m} = 2 / 3 (\frac{f_{t}}{f_{c}})$	$\cos θ = 1$
Compression Strength under Confinement Pressure	$(f_{b}, f_{b}, 0)$	$f_{b} =1 .2 f_{c}$	$r = \sqrt{2 / 3} \frac{f_{2} - s_{0}}{f_{c}}$	$σ_{m} = \frac{f_{2} + 2 s_{0}}{3 f_{c}}$	$\cos θ = - 1$

Figure 19 and Figure 20 show the points that determine the failure surface.

From these plots that the failure envelope is not a convex surface. Figure 21 shows this behavior.

In this particular case, the compressive strength is changing but all other ratios are fixed

\frac{f_{t}}{f_{c}}, \frac{f_{b}}{f_{c}}, \frac{f_{2}}{f_{c}}, \frac{s_{0}}{f_{c}}

. This leads to an envelope scaling, as shown in Figure 24.

Here with same strength in LAW24, but different confined compression strength

f_{2}

f_{c}

and the ratios

\frac{f_{t}}{f_{c}}, and \frac{f_{b}}{f_{c}}

in the

r - σ_{m}

space (used to define the concrete failure) are:

Where the failure curve is defined using $r = \sqrt{2 J_{2}} = \sqrt{\frac{2}{3}} σ_{V M}$ and $σ_{m} = \frac{I_{1}}{3}$ is the mean stress (pressure), then $I_{1}$ and $J_{2}$ are the first and second stress invariants.

The material fails once it reaches the failure curve $r_{f}$ .

Concrete Reinforcement

In Radioss there are two different ways to simulate the reinforcement in concrete.

One way is to use beam or truss elements and connect them to the concrete with kinematic conditions.
Another way is to use the parameters in LAW24 along with the orthotropic solid property /PROP/TYPE6 to define the reinforced direction. Parameters $α_{1}, α_{2}, α_{3}$ in LAW24 are used to define the reinforcement cross-section area ratio to the whole concrete section area in direction 1, 2, 3.(13) $α_{i} = \frac{A r e a_{s t e e l}}{A r e a_{c o n c r e t e}}$

Where,

σ_{y}

is the yield stress of the reinforcement. If steel is used as a reinforcement, then

σ_{y}

is the yield stress of steel and

E_{t}

is the modulus of steel in the plastic phase.

Concrete Material (/MAT/LAW81)

LAW81 can be used to model rock or concrete materials.

Drücker-Prager Yield Criteria

LAW81 uses a Drücker–Prager yield criterion where the yield surface and the failure surface are the same. The yield criteria is:(14)

F = \underset{J_{2} part}{\underset{︸}{q}} - \underset{I_{1} part}{\underset{︸}{r_{c} (p) \cdot (p \tan ϕ + c)}} = 0

Where,

$q$: von Mises stress with $q = σ_{V M} = \sqrt{3 J_{2}}$
$p$: Pressure is defined as $p = \frac{1}{3} I_{1}$

The yield surface can be described in two parts:

The linear part ( $p \leq p_{a}$ ), where the scale function is $r_{c} (p) = 1$ which leads to the von Mises stress being linearly proportional to pressure:(15) $q = p \tan ϕ + c$
Where,

$c$

Cohesive and is the intercept of yield envelope with the shear strength.

If $c = 0$ , the material has no strength under tension.

$ϕ$

Angle of internal friction, which defines the slope of the yield envelope.

$c$ and $ϕ$ are also used to define the Mohr-Coulomb yield surface. The Drücker-Prage yield surface is a smooth version of the Mohr-Coulomb yield surface.

The second part ( $p_{a} < p < p_{b}$ ) of the yield surface simulates a cap limit. An increase of pressure in a rock or concrete material will increase the yield of the material; but, if pressure increases enough, then the rock or concrete material will be crushed. The Drücker-Prager model with the cap limit can be used to model this behavior. The cap limit defined in part and uses the scale function:(16) $r_{c} (p) = \sqrt{1 - {(\frac{p - p_{a}}{p_{b} - p_{a}})}^{2}}$
The von Mises stress is:(17) $q = \sqrt{1 - {(\frac{p - p_{a}}{p_{b} - p_{a}})}^{2}} \cdot (p \tan ϕ + c)$ Where,

$p_{b}$

Curve is defined using the fct_ID_Pb input

$p_{a}$

Computed by Radioss using the input $α$ ratio value.

$p_{a} = α \cdot p_{b}$ with $0 < α < 1$ .

Where, $p_{0}$ is the maximum point of yield curve, where $\frac{\partial F}{\partial p} (p_{0}) = 0$

If $p = p_{b}$ , then $r_{c} (p_{b}) = 0$ and the yield function is then,

$q = 0 \cdot (p \tan ϕ + c) = 0$ which means the material is crushed.

The input parameters

ϕ, c, p_{b}, α

need to determine for the Drücker–Prager yield surface. At least four tests are needed to fit these parameters. In the simplest case, uniaxial tension and uniaxial compression can be used to determine the linear part,

ϕ, and c

. To determine

p_{b}, and α

biaxial compression tests and compression/compression tests are needed (refer to CC00 and CC01 in RD-E: 4701 Concrete Validation with Kupfer Tests).

For most materials such as metal, the plastic strain increment could be considered normal to yield surface. However, if the plastic strain increment normal to yield surface is used for rock or concrete materials, the plastic volume expansion is overestimated. Therefore, a non-associated plastic flow rule is used in these materials. In LAW81 the plastic flow function

G

defined as:

$G = q - p \cdot \tan ψ = 0$ if $p \leq p_{a}$
$G = q - \tan ψ (p - \frac{{(p - p_{a})}^{2}}{2 (p_{0} - p_{a})}) = 0$ if $p_{a} < p \leq p_{0}$
$G = F$ if $p > p_{0}$

Since the pressure is

p_{0}

, the yield function

F

and plastic flow function

G

are the same and the following condition is fulfilled:(18)

G (p_{0}) = F (p_{0})

(19)

\frac{\partial G}{\partial p} |_{p_{0}} = \frac{\partial F}{\partial p} |_{p_{0}} = 0

The pressure

p_{0}

can be calculated using the yield surface where

\frac{\partial F}{\partial p} |_{p_{0}} = 0

. With

G

defined as:(20)

G = q - \tan ψ (p - \frac{{(p - p_{a})}^{2}}{2 (p_{0} - p_{a})}) = 0

The parameter

ψ

can be determined using the von Mises stress at pressure,

p_{0}

in the function.

¹ Han, D. J., and Wai-Fah Chen. "A nonuniform hardening plasticity model for concrete materials." Mechanics of materials 4, no. 3-4 (1985): 283-302