FLAC3D Theory and Background • Constitutive Models

Concrete Model

This concrete model is a plastic-damage model extended and modified from the work in Lublinear et al (1989) and Lee & Fenves (1998) using the concepts of fracture-energy-based damage and modulus degradation in continuum damage mechanics.

Formulations

Modulus Degradation

Based on the concept of the continuum damage theory , the stress \(\sigma\) is linked with the effective stress \(\bar{\sigma}\) [note: in this model, the concept of “effective stress” is a concept in damage theory, which is different from the concept of Terzaghi’s effective stress in fluid-machinal interaction. A variable with a bar on top is in terms of effective stress] by an isotropic damage factor \(D\):

(1)\[\sigma_{ij} = (1-D) \bar{\sigma}_{ij}\]

together with Hookie’s law \(\bar{\sigma}_{ij} = E_{ijkl} \epsilon_{kl}^e\) and additive strain decomposition \(\epsilon_{ij} = \epsilon_{ij}^e + \epsilon_{ij}^p\), it obtains that

(2)\[\sigma_{ij} = (1-D) \bar{\sigma}_{ij} = (1-D) E_{ijkl} (\epsilon_{kl} - \epsilon_{kl}^p)\]

In other words, the nominal elastic modulus is in degradation by a factor of \((1-D)\) from the elastic modulus if relating nominal stress \(\sigma\) and strain \(\epsilon\): \(E \to (1-D) E\).

Damage Based on Fracture-Energy

In this model, two damage variables, one for tensile damage and the other for compressive damage, are defined independently, and each variable is factorized into the effective-stress response and stiffness degradation response parts. Considering a uniaxial tensile or compressive stress state, the state variable \(\chi \in [c,t]\) is introduced. The state is uniaxial tensile for \(\chi = t\) and uniaxial compressive for \(\chi = c\).

The stress is assumed to be determined by the plastic strain. The relation between the uniaxial stress, denoted by \(\sigma_{\chi}\), and the corresponding scalar plastic strain, denoted by \(\epsilon_{\chi}^p\), is assumed as

(3)\[\sigma_{\chi} = f_{\chi 0} \left[(1+a_{\chi}) \exp{(-b_{\chi} \epsilon_{\chi}^p)} - a_{\chi} \exp{(-2b_{\chi} \epsilon_{\chi}^p)} \right]\]

where \(a_{\chi}\) and \(b_{\chi}\) are positive material constants to be determined, and \(f_{\chi 0}\) is the initial yield stress without damage.

Suppose that we have available experimentally derived stress-strain diagrams in uniaxial tension and compression (top ones in Figure 1 & 2), and that these can be converted into \(\sigma_{\chi}\)-\(\epsilon_{\chi}^p\) curves (bottom ones in Figure 1 & 2). Assume that the area under these curves are infinite and equal to \(g_{\chi}\), define

(4)\[\kappa_{\chi} = \frac{1}{g_{\chi}} \int_{0}^{\epsilon_{\chi}^p} {\sigma_{\chi} \, d \epsilon_{\chi}^p}\]

where apparently \(0 \le \kappa_{\chi} \le 1\). Stress in terms of \(\epsilon_{\chi}^p\), or \(\sigma_{\chi} = f(\epsilon_{\chi}^p)\), can rewritten (Lublinear et al, 1989) in terms of \(\kappa_{\chi}\), or

(5)\[\sigma_{\chi} = f(\kappa_{\chi}) = \frac{f_{\chi 0}}{a_{\chi}} \left[ (1 + a_{\chi}) \sqrt{\phi(\kappa_{\chi})} - \phi(\kappa_{\chi}) \right]\]

where \(\phi(\kappa_{\chi}) = 1 + a_{\chi}(2 + a_{\chi}) \kappa_{\chi}\).

Figure 1: Typical uniaxial compression curves: top - uniaxial compressive stress vs compressive strain; bottom - uniaxial compressive stress vs plastic compressive strain.

Figure 2: Typical uniaxial tension curves: top - uniaxial tensile stress vs tensile strain; bottom - uniaxial tensile stress vs plastic tensile strain.

For stress defined by Equation (3), the quantity \(g_{\chi}\), which denotes the dissipated energy density during the entire process of microcracking, has a closed-form expression of the area under the curve of A-B-D in Figure 1-(b) or 2-(b):

(6)\[g_{\chi} = \int_{0}^{+\infty} {\sigma_{\chi} \, d \epsilon_{\chi}^p} = \frac{f_{\chi 0}}{b_{\chi}} \left( 1 + \frac{a_{\chi}}{2} \right)\]

The initial slope of the curve \(\sigma_{\chi}\)-\(\epsilon_{\chi}^p\) is

(7)\[\left. \frac{d \sigma_{\chi}}{d \epsilon_{\chi}^p} \right|_{\epsilon_{\chi}^p = 0} = f_{\chi 0} b_{\chi} (a_{\chi} - 1)\]

Note that \(a_{\chi} > 1\) implies initial hardening (e.g., when \(\chi = c\)) and \(a_{\chi} < 1\) implies softening immediately after yielding (e.g., when \(\chi = t\)). When \(a_{\chi} > 1\), \(\sigma_{\chi}\) obtain a maximum value (see Point A in Figure 1)

(8)\[f_{\chi}^m = f_{\chi 0} \frac{(1+a_{\chi})^2}{4a_{\chi}}\]

at

(9)\[\epsilon_{\chi}^{p,m} = - \frac{\ln{\left(\frac{1+a_{\chi}}{2a_{\chi}}\right)}}{b_{\chi}}\]

For a curve with initial hardening after yielding which is typical in uniaxial compression, the material constant \(a_{\chi}\) can be back-calculated from Equation (8) by

(10)\[a_{\chi} = 2 (f_{\chi}^m/f_{\chi 0}) - 1 + 2 \sqrt{(f_{\chi}^m/f_{\chi 0})^2 - (f_{\chi}^m/f_{\chi 0})}\]

For a curve with immediate softening after yielding which is typical in uniaxial extension, a simple way is assuming \(a_t = 0\), which leads to

(11)\[\sigma_t = f_{t0} \exp{(-b_t \epsilon_t^p)}\]

and

(12)\[g_t = f_{t0} / b_t\]

Once \(a_{\chi}\) is determined or assumed, \(b_{\chi}\) can be obtained from Equation (6) provided \(g_{\chi}\) and \(f_{\chi 0}\) are known.

Because the capacity for the dissipated energy per unit volume cannot be given as a material property, \(g_{\chi}\) should be derived from other known material properties such as the fracture energy (Lee & Fenvas, 1998). Assuming that the fracture energy in the uniaxial tensile state \(G_{\chi}\) and its counterpart in the uniaxial compressive state are given as a material properties, \(g_{\chi}\) is the specific fracture energy normalized by the localization zone size, also referred to as the characteristic length \(l_{\chi}\), which leads to (Lubliner et al. 1989)

(13)\[g_{\chi} = \frac{G_{\chi}}{l_{\chi}}\]

To maintain objective results at the structural level, \(l_{\chi}\) should be an objective value or assumed as a material property.

Yield Function

The yield surface is schemed by Figure Figure #model-concrete-yd in the biaxial-stress (plane-stress) space. Note that we assume the principle stresses with the order \(\sigma_1 \le \sigma_2 \le \sigma_3\) and positive for tension. Due to symmetry, only A-B-C-D is considered.

For line A-B, a tension yield criteria is used:

(14)\[f^t (\bar{\sigma}) = \bar{\sigma}_3 - \bar{\sigma}_t = 0\]

For line B-C, a Mohr-Coulomb yield criteria is used:

(15)\[f^s (\bar{\sigma}) = - \bar{\sigma}_t \bar{\sigma}_1 + \bar{\sigma}_c \bar{\sigma}_3 - \bar{\sigma}_c \bar{\sigma}_t= 0\]

For curve C-D, a Drucker-Prager yield criteria is used:

(16)\[f^d (\bar{\sigma}) = \alpha \bar{I}_1 + \bar{q} - (1-\alpha) \bar{\sigma}_c = 0\]

where \(\bar{I}_1\) is the first stress invariant, \(\bar{q}\) is the equivalent shear stress, and

(17)\[\alpha = \frac{f_{b0}/f_{c0} - 1}{2(f_{b0}/f_{c0}) - 1}\]

where \(f_{c0}\) is the uniaxial yield compressive strength and \(f_{b0}\) is the biaxial yield compressive strength. Since usually \(f_{b0} \ge f_{c0}\), \(\alpha\) is in the range of \(0 \le \alpha < 0.5\).

Figure 3: Schematic of initial yield surface in the biaxial-stress (plane-stress) space.

Flow Rule

The plastic strain rate is evaluated by the flow rule, which is defined by a scalar plastic potential function \(g\). For a plastic potential in the effective stress space, the plastic strain is given by

(18)\[\Delta \epsilon_{ij}^p = \lambda \frac{\partial g}{\partial \bar{\sigma}_{ij}}\]

where \(\lambda\) is a non-negative function referred to as the plastic consistency parameter (plastic multiplier).

For yielding at line A-B, an associated potential function is used:

(19)\[g^t (\bar{\sigma}) = \bar{\sigma}_3 - \bar{\sigma}_t\]

For yielding at line B-C, a nonassociated Mohr-Coulomb type potential function is used:

(20)\[g^s (\bar{\sigma}) = - \bar{\sigma}_1 + N_{\psi} \bar{\sigma}_3\]

where \(N_{psi} = \min(\bar{\sigma}_c/\bar{\sigma}_t, (1+\sin{\psi})/(1-\sin{\psi})\) and \(\psi\) is the dilation angle which is an input material property.

For yielding at curve C-D, a nonassociated Drucker-Prager yield criteria is used:

(21)\[g^d (\bar{\sigma}) = \alpha_p \bar{I}_1 + \bar{q}\]

where \(\alpha_p = min(\alpha, 6\sin{\psi}/(3-\sin{\psi}))\).

Another internal variable set, besides the plastic strain, is needed to represent the damage states. The damage variable, \(\kappa_{\chi}\), is assumed to be the only necessary state variable and its evolution is expressed as

(22)\[\Delta \kappa_t = r(\bar{\sigma}) \frac{f_t}{g_t} \Delta \epsilon_3^p\]

(23)\[\Delta \kappa_c = (1-r(\bar{\sigma})) \frac{f_c}{g_c} \Delta \epsilon_1^p\]

where \(\epsilon_1^p\) is the minimum principal plastic strain, \(\epsilon_3^p\) is the maximum principal plastic strain; \(r\) is a weight function of the effective principal stresses, which ranges from 0 in pure compressive loading to 1 in pure tensile loading:

(24)\[r(\bar{\sigma}) = \frac{\sum_{i=1}^3 \langle \bar{\sigma}_i \rangle}{\sum_{i=1}^3 | \bar{\sigma}_i |}\]

Evolution of Damage

To model the different damage states for tensile and compressive loading, both \(\kappa_t\) and \(\kappa_c\) are used as independent damage variables. The evolution equations for the hardening and degradation variables are obtained from the factorization of the uniaxial stress strength function, \(f_{\chi}=f_{\chi}(\kappa_{\chi})\), such that

(25)\[f_{\chi} = (1 - D_{\chi}) \bar{f}_{\chi}\]

in which \(0 \le D_{\chi} < 1\) and \(\Delta D_{\chi} > 0\). \(\bar{\sigma}_{\chi} = \bar{f}_{\chi}\) describes the evolution of the cohesion variable used in the yield function, and \(D_{\chi} = D_{\chi}(\kappa_{\chi})\) defines the degradation damage as a function of the damage variable:

(26)\[D = 1 - (1- D_c)(1 - sD_t)\]

where

(27)\[s = s_0 + (1 - s_0) r(\bar{\sigma})\]

with \(0 \le s_0 \le 1\) and so \(s_0 \le s \le 1\) represents stiffness recovery.

The degradation is also assumed an exponent form as

(28)\[1 - D_{\chi} = \exp{( - d_{\chi} \epsilon_{\chi}^p)}\]

where \(d_{\chi}\) are material constants, which can be determined form the unloading-reloading (see B-C lines in Figure 1-a & 2-a).

Material Properties

In this model, two material properties, initial Young’ modulus and Poisson’s ratio, (\(E\), \(\nu\)), or initial bulk and shear modulus, (\(K\), \(G\)), are required to determine the elasticity.

To match the uniaxial compression and tension experimental curves, one straightforward approach is to specify the initial and peak (maximum) strengths (\(f_{\chi 0}\) and \(f_{\chi}^m\)), fracture energies (\(G_{\chi}\)), characteristic lengths (\(l_{\chi}\)), and damages at specific stresses (\(D_c^{peak}\) and \(D_t^{half}\)) in the uniaxial compression/tension states. These input will determine the parameters \(a_{\chi}\), \(b_{\chi}\), \(d_{\chi}\) and \(g_{\chi}\).

Once \(f_{\chi 0}\) and \(f_{\chi}^m\) are known, there are two cases:

• if the curve is hardening and then softening (usually in uniaxial compression experiments, see Figure 1), or in other words, \(f_{\chi}^m > f_{\chi 0}\), then \(a_{\chi}\) can be determined by Equation (10).

• if the curve is immediately softening (usually in uniaxial extension experiments, see Figure 2), or in other words, \(f_{\chi}^m = f_{\chi 0}\), then \(a_{\chi}\) will be assumed zero.

Once \(G_{\chi}\) and \(l_{\chi}\) are known, then \(g_{\chi} = G_{\chi}/l_{\chi}\), \(b_{\chi}\) can then be determined by Equation (6) so that \(b_{\chi}=f_{\chi 0}(1+a_{\chi}/2)/g_{\chi}\).

Equation (28) can determine \(d_{\chi}\):

• for uniaxial compression where the curve is hardening and then softening, then \(d_c\) can be calculated by the the damage at the peak stress,

(29)\[d_c = b_c \frac{\ln{(1-D_c^{peak})}}{\ln{\left( \frac{1+a_c}{2a_c} \right)}}\]

• for uniaxial extension where the curve is immediately softening, then \(d_t\) can be calculated by the the damage at the half peak/initial stress,

(30)\[d_t = b_t \frac{\ln{(1-D_t^{half})}}{-\ln{2}}\]

\(\alpha_p\) is a material property determining the dilation, with typical value \(\alpha_p=0.2\).

The biaxial compression initial strength \(f_{b0}\) is requited to calculate \(\alpha\) in Equation (17). Here the ratio \(f_{b0}/f_{c0}\) is instead as an input material property, with a typical value \(f_{b0}/f_{c0} = 1.05{\sim}1.2\).

The stiffness recovery parameter \(s_0\) is usually in the range of \(0 \le s_0 \le 0.2\).

References

Lee, J., & Fenves, G. L. (1998). Plastic-damage model for cyclic loading of concrete structures. Journal of engineering mechanics, 124(8), 892-900.

Lee, J., & Fenves, G. L. (1998). A plastic damage concrete model for earthquake analysis of dams. Earthquake engineering & structural dynamics, 27(9), 937-956.

Lubliner, J., Oliver, J., Oller, S., & Onate, E. (1989). A plastic-damage model for concrete. International Journal of solids and structures, 25(3), 299-326.

concrete Model Properties

Use the following keywords with the zone property command to set these properties of the concrete model.

concrete

bulk f

bulk modulus, \(K\).

compression-strength f

uniaxial compression peak/maximum strength, \(f_c^m\).

compression-initial f

uniaxial compression initial yield strength, \(f_{c0}\).

compression-energy-fracture f

fracture energy in uniaxial compression state, \(G_c\).

compression-length-reference f

reference (characteristic) length in uniaxial compression state, \(l_c\).

compression-d f

compression damage parameter, \(d_c\).

dilation f

dilation angle, \(\psi\).

poisson f

Poisson’s ratio, \(\nu\).

ratio-biaxial f

strength ratio of biaxial to uniaxial compression initial, usually in a range of (1,1.25], \(f_{b0}/f_{c0}\).

shear f

shear modulus, \(G\).

tension-strength f

uniaxial tension peak/maximum strength, \(f_t^m\).

tension-initial f

uniaxial tension initial yield strength, \(f_{t0}\). Usually, \(f_{t0} = f_t^m\).

tension-energy-fracture f

fracture energy in uniaxial compression state, \(G_t\).

tension-length-reference f

reference (characteristic) length in uniaxial tension state, \(l_t\).

tension-d f

tension damage parameter, \(d_t\).

tension-recovery f

stiffness recovery parameter for tension, typical value is 0.2, \(s_0\).

young f

Young’s modulus, \(E\).

compression f (r)

Current compression strength, \(\sigma_c\).

damage f (r)

Current total damage, \(D\).

damage-compression f (r)

Current damage in compression, \(D_c\).

damage-tension f (r)

Current damage in tension, \(D_t\).

strain-plastic-compression f (r)

Current plastic strain in compression, \(\epsilon^p_c\).

strain-plastic-tension f (r)

Current plastic strain in tension, \(\epsilon^p_t\).

tension f (r)

Current tension strength, \(\sigma_t\).

Notes

• Only one of the two options is required to define the elasticity: bulk modulus \(K\) and shear modulus \(G\), or, Young’s modulus \(E\) and Poisson’s ratio \(\nu\). When choosing the latter, Young’s modulus \(E\) must be assigned in advance of Poisson’s ratio \(\nu\).

Key

(a) Advanced property.

This property has a default value; simpler applications of the model do not need to provide a value for it.

(r) Read-only property.

This property cannot be set by the user. Instead, it can be listed, plotted, or accessed through FISH.