An inelastic constitutive model for concrete

The model uses the classical concepts of plasticity theory: a strain rate decomposition into elastic and inelastic strain rates, elasticity, yield, flow, and hardening.

We begin with a strain rate decomposition:

where $d\mathit{\epsilon }$ is the total mechanical strain rate, $d{\mathit{\epsilon }}^{el}$ is the elastic strain rate (which includes crack detection strains—this elastic strain will be further decomposed when we describe the elasticity), and $d{\mathit{\epsilon }}_c^{pl}$ is the plastic strain rate associated with the “compression” surface.

We assume that the elastic part of the strain is always small, so this equation can be integrated as

The “compression” surface is

where p is the effective pressure stress, defined as

and q is the Mises equivalent deviatoric stress:

where $\mathbf{S}=\mathit{\sigma }+p\mathbf{I}$ are the deviatoric stress components; $a_0$ is a constant, which is chosen from the ratio of the ultimate stress reached in biaxial compression to the ultimate stress reached in uniaxial compression; and ${\tau }_c\left({\lambda }_c\right)$ is a hardening parameter (${\tau }_c$ is the size of the yield surface on the q-axis at $p=0$, so that ${\tau }_c$ is the yield stress in a state of pure shear stress when all components of $\mathit{\sigma }$ are zero except ${\sigma }_{12}={\sigma }_{21}={\tau }_c$). The hardening is measured by the value of ${\lambda }_c$: the ${\tau }_c\left({\lambda }_c\right)$ relationship is user-specified.

This simple surface is a straight line in (p–q) space and provides a good match to experimental data over a fairly wide range of pressure stress values (up to four to five times the maximum compressive stress that can be carried by the concrete in uniaxial compression). This form of the surface means that, as the hardening (${\lambda }_c$) changes, the surfaces in (p–q) space are similar, so—for example—the ratio of flow stress in biaxial loading to flow stress in uniaxial loading is the same at all flow stress levels. This does not appear to be contradicted by any experimental data, and it means that only one constant ($a_0$) is needed to define the shape of the surface.

The value of $a_0$ is established from the user's data as follows. In uniaxial compression $p=\frac{1}{3}{\sigma }_c$ and $q={\sigma }_c$, where ${\sigma }_c$ is the stress magnitude. Therefore, on $f_c=0$,

In biaxial compression $p=2/3 {\sigma }_{bc}$ and $q={\sigma }_{bc}$, where ${\sigma }_{bc}$ is the magnitude of each nonzero principal stress. Therefore, on $f_c=0$,

The value of ${\sigma }_{bc}^u/{\sigma }_c^u=r_{bc}^{\sigma }$ is specified by the user as part of the failure surface data (typically $r_{bc}^{\sigma }\approx 1.16$). $a_0$ can be calculated from Equation 4 and Equation 5 as

The “compression” surface is shown in Figure 2 and Figure 3.

The user defines the hardening by specifying the magnitude of the stress, $\left|{\sigma }_{11}\right|$, in a uniaxial compression test as a function of the inelastic strain magnitude, $\left|{\epsilon }_{11}\right|$. These data are used to define the ${\tau }_c\left({\lambda }_c\right)$ relationship, as follows.

In uniaxial compression $p=\frac{1}{3}{\sigma }_c$ and $q={\sigma }_c$, where ${\sigma }_c$ is the stress magnitude. During active plastic loading $f_c=0$, so by using the definition of $f_c$ (Equation 4), we obtain ${\tau }_c$ immediately as

The model uses associated flow, so if $f_c=0$ and $d{\lambda }_c>0$,

otherwise, $d{\mathit{\epsilon }}_c^{pl}=0$.

In the definition of $d{\mathit{\epsilon }}_c^{pl}$, $c_0$ is a constant that is chosen so that the ratio of ${\epsilon }_{11}^{pl}$ in a monotonically loaded biaxial compression test to ${\epsilon }_{11}^{pl}$ in a monotonically loaded uniaxial compression test is $r_{bc}^{\epsilon }$, a value specified by the user as part of the failure surface data (typically $r_{bc}^{\epsilon }\approx 1.28$). The equation defining $c_0$ from $r_{bc}^{\epsilon }$ and the other constants in the compression surface is derived next.

The gradient of the flow potential for the compressive surface is

Since

and

then

In uniaxial compression $p=\frac{1}{3}{\sigma }_c$, $q={\sigma }_c$, and $S_{11}=-\frac{2}{3}{\sigma }_c$, so Equation 7 defines

This equation can be integrated immediately to give

so ${\lambda }_c$ is known from ${\left({\epsilon }_c^{pl}\right)}_{11}^c$ and the constants $a_0$ and $c_0$. Equation 6 and Equation 9, therefore, define the ${\tau }_c\left({\lambda }_c\right)$ relationship from the concrete model input data once $c_0$ is known.

The constant $c_0$ is calculated from the user's definition of $r_{bc}^{\epsilon }$, the ratio of ${\left({\epsilon }_c^{pl}\right)}_{11}^{bc}$ to ${\left({\epsilon }_c^{pl}\right)}_{11}^c$, the total plastic strain components that would occur in monotonically loaded biaxial and uniaxial compression tests. In biaxial compression when both nonzero principal stresses have the magnitude ${\sigma }_{bc}$, $p=\frac{2}{3}{\sigma }_{bc}=\frac{2}{3}{r}_{bc}^{\sigma }{\sigma }_c$, $q={\sigma }_{bc}=r_{bc}^{\sigma }{\sigma }_c$, and $S_{11}=-\frac{1}{3}{r}_{bc}^{\sigma }{\sigma }_c$, so the flow rule gives

Using this equation and Equation 8 then defines $c_0$ from $r_{bc}^{\epsilon }$ and the other constants as

Cracking dominates the material behavior when the state of stress is predominantly tensile. The model uses a “crack detection” plasticity surface in stress space to determine when cracking takes place and the orientation of the cracking. Damaged elasticity is then used to describe the postfailure behavior of the concrete with open cracks.

Numerically we use the “crack detection” plasticity model for the increment in which cracking takes place and subsequently use damaged elasticity once the crack's presence and orientation have been detected. As a result there is at least one increment in which we calculate crack detection “plastic” strains. Since these are really just the outcome of a numerical device to treat cracking, they are recast as elastic strains in the direction of cracking and as plastic strains in the other directions. (This means that we retain the stresses calculated for equilibrium purposes, as well as the strain decomposition of Equation 1.)

The basis of the postcracked behavior is the brittle fracture concept of Hilleborg (1976). We assume that the fracture energy required to form a unit area of crack surface, $G_f$, is a material property. This value can be calculated from measuring the tensile stress as a function of the crack opening displacement (Figure 4), as

Figure 4. Cracking behavior based on fracture energy.

Typical values of $G_f$ range from 40 N/m (0.22 lb/in) for a typical construction concrete (with ${\sigma }_c^u\approx$ 20 MPa, 2850 lb/in²) to 120 N/m (0.67 lb/in) for a high strength concrete (with ${\sigma }_c^u\approx$ 40 MPa, 5700 lb/in²).

The implication of assuming that $G_f$ is a material property is that, when the elastic part of the displacement, $u^{el},$ is eliminated, the relationship between the stress and the remaining part of the displacement, $u^{cr}=u-u^{el},$ is fixed, regardless of the specimen size. For example, consider a specimen developing a single crack across its section as tensile displacement is applied to it: $u^{cr}$ is the displacement across the crack and is not changed by using a longer or shorter specimen in the test (so long as the specimen is significantly longer than the width of the crack band, which will typically be of the order of the aggregate size). Thus, one important part of the cracked concrete's tensile behavior is defined in terms of a stress/displacement relationship. In the finite element implementation of this model we must, therefore, compute the relative displacement at an integration point to provide $u^{cr}$. We do this in Abaqus by multiplying the strain by a characteristic length associated with the integration point. The characteristic crack length is based on the element geometry and formulation: it is a typical length of a line across an element for a first-order element; it is half of the same typical length for a second-order element. For beams and trusses it is a characteristic length along the element axis. For membranes and shells it is a characteristic length in the reference surface. For axisymmetric elements it is a characteristic length in the r–z plane only. For cohesive elements it is equal to the constitutive thickness. This definition of the characteristic length is used because we do not necessarily know in which direction the concrete will crack and so cannot choose the length measure in any particular direction. Thus, if there are elements in the model that have large aspect ratios, the model will likely provide different results if it is loaded in different directions and cracking occurs in such elements. This effect should be considered by the user in defining values for the material properties.

In reinforced concrete the $\sigma$–$u^{cr}$ relationship must also represent the action of the bond between the concrete and the rebar as the concrete cracks. We assume this is accommodated by increasing the value of $G_f$ based on comparisons with experiments on reinforced material.

We first describe the crack detection plasticity model and then discuss the damaged elasticity.

We decompose the elastic strain rate of Equation 1 as

where $d{\mathit{\epsilon }}^{el}$ is the total mechanical strain rate for the crack detection problem, $d{\mathit{\epsilon }}_d^{el}$ is the elastic strain rate, and $d{\mathit{\epsilon }}_t^{pl}$ is the plastic strain rate associated with the crack detection surface.

The crack detection surface is the Coulomb line

where ${\sigma }_t^u$ is the failure stress in uniaxial tension and $b_0$ is a constant that is defined from the value of the tensile failure stress, ${\sigma }_I$, in a state of biaxial stress when the other nonzero principal stress, ${\sigma }_{II}$, is at the uniaxial compression ultimate stress value, ${\sigma }_c^u$. ${\sigma }_t\left({\lambda }_t\right)$ is a hardening parameter (${\sigma }_t$ is the equivalent uniaxial tensile stress). The hardening is measured by ${\lambda }_t$, with the ${\sigma }_t\left({\lambda }_t\right)$ relationship defined from the user-specified tension stiffening data (see Figure 5).

Figure 5. “Tension stiffening” model.

The stress measures $\widehat{p}$ and $\widehat{q}$ are defined in the same way as p and q, except that all stress components ${\sigma }_{\alpha \beta }$ associated with open cracks (that is, if $\alpha$ or $\beta$ is a crack direction in which the direct strain ${\epsilon }_{\overline{\alpha \alpha }}^{el}>0$ or ${\epsilon }_{\overline{\beta \beta }}^{el}>0$) are not included in these measures: they are invariants in subspaces of the stress space.

This surface has a simple mathematical form but matches plane stress data quite well. The hardening is introduced in the particular form shown in Equation 11 so that, as ${\sigma }_t\to 0$, the surface becomes $q-3p=0$, which in (p–q) space is the cone containing the principal axes of stress. This means that, as the tension stiffening is exhausted in a plane stress test, the stress point will drop back onto the nearest principal stress axis.

The value of $b_0$ is obtained as follows. The failure surface data include a definition of f, a ratio that states that, in a plane stress test cracking would occur when one principal stress has the value ${\sigma }_I=-{\sigma }_c^u$ (${\sigma }_c^u$ is the magnitude of the ultimate stress in uniaxial compression) and the other nonzero principal stress has the value ${\sigma }_{II}=f{\sigma }_t^u$. Another value provided as part of the failure surface data is $r_t^{\sigma }$, which defines

Cracking would, therefore, occur at the point with principal stresses $-{\sigma }_c^u$, $f{r}_t^{\sigma }{\sigma }_c^u$, and 0. For these values

and

Therefore, with ${\sigma }_t={\sigma }_t^u$,

so

The crack detection surface is shown in Figure 2 and Figure 3.

The crack detection model uses the assumption of associated flow, so if $f_t=0$ and $d{\lambda }_t>0$,

otherwise, $d{\mathit{\epsilon }}_t^{pl}=0$.

The user defines the tension stiffening behavior by specifying the magnitude of the stress, ${\sigma }_t$, in a uniaxial tension test as a function of the inelastic strain. (When the fracture energy concept is used to define the postfailure behavior, “strain” is now defined as $u^{cr}/c$, where c is the characteristic length associated with the integration point.) The ${\sigma }_t\left({\lambda }_t\right)$ relationship is defined as follows.

Using the definition of $f_t$, Equation 11, in the flow rule above we write

In uniaxial tension $S_{11}=\frac{2}{3}{\sigma }_t$ and $q={\sigma }_t$. Therefore, in uniaxial tension,

and, hence,

Upon integration Equation 13 gives ${\lambda }_t$ from ${\left(d{\epsilon }_t^{pl}\right)}_{11}$; and, therefore, the ${\sigma }_t\left({\lambda }_t\right)$ relationship is obtained from the tension stiffening input data.

Following crack detection we use damaged elasticity to model the failed material. The elasticity is written in the form

where $\mathbf{D}$ is the elastic stiffness matrix for the concrete.

Let $\alpha$ represent a cracked direction, with corresponding direct stress ${\sigma }_{\overline{\alpha \alpha }}$ and direct elastic strain ${\epsilon }_{\overline{\alpha \alpha }}^{el}$. In these expressions and in the remainder of this section, no summation is implied by repeated indices with a bar over them. If the fracture energy concept is used, the strains are related to the user-specified stress/displacement definition for the tension stiffening behavior by $\epsilon =u/c$, where c is the characteristic length associated with the integration point.

Then, in $\mathbf{D}$, $D_{\overline{\alpha \alpha }\beta \gamma }$ is the usual elasticity of the concrete if ${\epsilon }_{\overline{\alpha \alpha }}\le 0$. If ${\epsilon }_{\overline{\alpha \alpha }}^{\mathrm{open}}>{\epsilon }_{\overline{\alpha \alpha }}>0$,

where ${\sigma }_{\overline{\alpha \alpha }}^{\mathrm{open}}$ is the stress corresponding to ${\epsilon }_{\overline{\alpha \alpha }}^{\mathrm{open}}$ (as defined in the tension stiffening data), and

If ${\epsilon }_{\overline{\alpha \alpha }}={\epsilon }_{\overline{\alpha \alpha }}^{\mathrm{open}}$,

which follows from the tension stiffening data.

We also assume no Poisson effect for open cracks: for ${\epsilon }_{\overline{\alpha \alpha }}^{el}>0$, $D_{\overline{\alpha \alpha }\beta \gamma }=0$ for $\beta \ne \alpha$, $\gamma \ne \alpha$.

The shear terms in the elasticity associated with existing crack directions are

where $\widehat{G}={\varrho }^{\mathrm{close}}G\text{ for }{\epsilon }_{\overline{\alpha \alpha }}\le 0,$ and $\widehat{G}={\varrho }^{\mathrm{open}}G\text{ for }{\epsilon }_{\overline{\alpha \alpha }}>0.$ In these expressions G is the elastic shear modulus, ${\varrho }^{\mathrm{close}}$ is a constant defined by the user as part of the shear retention data (see Figure 6),

Figure 6. Shear retention.

and ${\varrho }^{\mathrm{open}}$ is a linear function of ${\overline{\epsilon }}^{el}$, ${\varrho }^{\mathrm{open}}=\left(1-{\overline{\epsilon }}_{\overline{\alpha \alpha }}^{el}/{\epsilon }^{\max }\right)$ and is also defined by the user in the shear retention data. Here ${\overline{\epsilon }}^{el}=⟨{\epsilon }_{\overline{\alpha \alpha }}^{el}⟩+⟨{\epsilon }_{\overline{\beta \beta }}^{el}⟩$, where $⟨$ and $⟩$ are Macauley brackets, defining

for any function f.

As soon as the crack detection surface ($f_t$) has been activated, we assume that cracking has occurred. The crack direction, ${\mathbf{n}}_{\alpha }$, is taken to be the direction of that part of the maximum principal plastic strain increment conjugate to the crack detection surface, $\mathrm{\Delta }{\mathit{\epsilon }}_t^{pl}$, that is orthogonal to the directions of any existing cracks at the same point. This crack orientation is stored for subsequent calculations, which are done for convenience in a local coordinate system oriented so that one of the coordinate directions is the crack direction, ${\mathbf{n}}_{\alpha }$. Cracking is irrecoverable in the sense that, once a crack has occurred at a point, it remains throughout the rest of the calculation. Following crack detection, the crack affects the calculations by damaging the elasticity, as defined above. In addition, if the elastic strain across a crack is tensile, the invariants used in the crack detection surface are defined in the stress sub-space in which all stress components associated with the open crack direction are neglected, as described in the section above on yield. This implies that no more than three cracks can occur at any point (two in a plane stress case, one in a uniaxial stress case).

The model is integrated using the backward Euler method generally used with the plasticity models in Abaqus. A material Jacobian consistent with this integration operator is used for the equilibrium iterations.