Damaged plasticity model for concrete and other quasi-brittle materials

This section describes the concrete damaged plasticity model provided in Abaqus for the analysis of concrete and other quasi-brittle materials.

The following topics are discussed:

ProductsAbaqus/StandardAbaqus/Explicit

The material library in Abaqus also includes other constitutive models for concrete based on the smeared crack approach. These are the smeared crack model in Abaqus/Standard, described in An inelastic constitutive model for concrete, and the brittle cracking model in Abaqus/Explicit, described in A cracking model for concrete and other brittle materials.

The concrete damaged plasticity model is primarily intended to provide a general capability for the analysis of concrete structures under cyclic and/or dynamic loading. The model is also suitable for the analysis of other quasi-brittle materials, such as rock, mortar and ceramics; but it is the behavior of concrete that is used in the remainder of this section to motivate different aspects of the constitutive theory. Under low confining pressures, concrete behaves in a brittle manner; the main failure mechanisms are cracking in tension and crushing in compression. The brittle behavior of concrete disappears when the confining pressure is sufficiently large to prevent crack propagation. In these circumstances failure is driven by the consolidation and collapse of the concrete microporous microstructure, leading to a macroscopic response that resembles that of a ductile material with work hardening.

Modeling the behavior of concrete under large hydrostatic pressures is out of the scope of the plastic-damage model considered here. The constitutive theory in this section aims to capture the effects of irreversible damage associated with the failure mechanisms that occur in concrete and other quasi-brittle materials under fairly low confining pressures (less than four or five times the ultimate compressive stress in uniaxial compression loading). These effects manifest themselves in the following macroscopic properties:

  • different yield strengths in tension and compression, with the initial yield stress in compression being a factor of 10 or more higher than the initial yield stress in tension;

  • softening behavior in tension as opposed to initial hardening followed by softening in compression;

  • different degradation of the elastic stiffness in tension and compression;

  • stiffness recovery effects during cyclic loading; and

  • rate sensitivity, especially an increase in the peak strength with strain rate.

The plastic-damage model in Abaqus is based on the models proposed by Lubliner et al. (1989) and by Lee and Fenves (1998). The model is described in the remainder of this section. An overview of the main ingredients of the model is given first, followed by a more detailed discussion of the different aspects of the constitutive model.

About the inviscid concrete damaged plasticity model

The main ingredients of the inviscid concrete damaged plasticity model are summarized below.

Strain rate decomposition

An additive strain rate decomposition is assumed for the rate-independent model:

ε˙=ε˙el+ε˙pl,

where ε˙ is the total strain rate, ε˙el is the elastic part of the strain rate, and ε˙pl is the plastic part of the strain rate.

Stress-strain relations

The stress-strain relations are governed by scalar damaged elasticity:

σ=(1-d)D0el:(ε-εpl)=Del:(ε-εpl),

where D0el is the initial (undamaged) elastic stiffness of the material; Del=(1-d)D0el is the degraded elastic stiffness; and d is the scalar stiffness degradation variable, which can take values in the range from zero (undamaged material) to one (fully damaged material). Damage associated with the failure mechanisms of the concrete (cracking and crushing) therefore results in a reduction in the elastic stiffness. Within the context of the scalar-damage theory, the stiffness degradation is isotropic and characterized by a single degradation variable, d. Following the usual notions of continuum damage mechanics, the effective stress is defined as

σ¯=defD0el:(ε-εpl).

The Cauchy stress is related to the effective stress through the scalar degradation relation:

σ=(1-d)σ¯.

For any given cross-section of the material, the factor (1-d) represents the ratio of the effective load-carrying area (i.e., the overall area minus the damaged area) to the overall section area. In the absence of damage, d=0, the effective stress σ¯ is equivalent to the Cauchy stress, σ. When damage occurs, however, the effective stress is more representative than the Cauchy stress because it is the effective stress area that is resisting the external loads. It is, therefore, convenient to formulate the plasticity problem in terms of the effective stress. As discussed later, the evolution of the degradation variable is governed by a set of hardening variables, ε~pl, and the effective stress; that is, d=d(σ¯,ε~pl).

Hardening variables

Damaged states in tension and compression are characterized independently by two hardening variables, ε~tpl and ε~cpl, which are referred to as equivalent plastic strains in tension and compression, respectively. The evolution of the hardening variables is given by an expression of the form

ε~pl=[ε~tplε~cpl];            ε~˙pl=h(σ¯,ε~pl)ε˙pl,

as described later in this section.

Microcracking and crushing in the concrete are represented by increasing values of the hardening variables. These variables control the evolution of the yield surface and the degradation of the elastic stiffness. They are also intimately related to the dissipated fracture energy required to generate micro-cracks.

Yield function

The yield function, F(σ¯,ε~pl), represents a surface in effective stress space, which determines the states of failure or damage. For the inviscid plastic-damage model

F(σ¯,ε~pl)0.

The specific form of the yield function is described later in this section.

Flow rule

Plastic flow is governed by a flow potential G according to the flow rule:

ε˙pl=λ˙G(σ¯)σ¯,

where λ˙ is the nonnegative plastic multiplier. The plastic potential is defined in the effective stress space. The specific form of the flow potential for the concrete damaged plasticity model is discussed later in this section. The model uses nonassociated plasticity, therefore requiring the solution of nonsymmetric equations.

Summary

In summary, the elastic-plastic response of the concrete damaged plasticity model is described in terms of the effective stress and the hardening variables:

(1)σ¯=D0el:(ε-εpl){σ¯|F(σ¯,ε~pl)0},ε~˙pl=h(σ¯,ε~pl)ε˙pl,ε˙pl=λ˙G(σ¯)σ¯,

where λ˙ and F obey the Kuhn-Tucker conditions: λ˙F=0;    λ˙0;    F0 The Cauchy stress is calculated in terms of the stiffness degradation variable, d(σ¯,ε~pl), and the effective stress as

(2)σ=(1-d)σ¯.

The constitutive relations for the elastic-plastic response, Equation 1, are decoupled from the stiffness degradation response, Equation 2, which makes the model attractive for an effective numerical implementation. The inviscid model summarized here can be extended easily to account for viscoplastic effects through the use of a viscoplastic regularization by permitting stresses to be outside the yield surface.

Damage and stiffness degradation

The evolution equations of the hardening variables ε~tpl and ε~cpl are conveniently formulated by considering uniaxial loading conditions first and then extended to multiaxial conditions.

Uniaxial conditions

It is assumed that the uniaxial stress-strain curves can be converted into stress versus plastic strain curves of the form

(3)σt=σt(ε~tpl,ε~˙tpl,θ,fi),σc=σc(ε~cpl,ε~˙cpl,θ,fi),

where the subscripts t and c refer to tension and compression, respectively; ε~˙tpl and ε~˙cpl are the equivalent plastic strain rates, ε~tpl=0tε~˙tpldt and ε~cpl=0tε~˙cpldt are the equivalent plastic strains, θ is the temperature, and fi,(i=1,2,) are other predefined field variables.

Under uniaxial loading conditions the effective plastic strain rates are given as

(4)ε~˙tpl=ε˙11pl,in uniaxial tension andε~˙cpl=-ε˙11pl,in uniaxial compression.

In the remainder of this section we adopt the convention that σc is a positive quantity representing the magnitude of the uniaxial compression stress; that is, σc=-σ11.

As shown in Figure 1, when the concrete specimen is unloaded from any point on the strain softening branch of the stress-strain curves, the unloading response is observed to be weakened: the elastic stiffness of the material appears to be damaged (or degraded). The degradation of the elastic stiffness is significantly different between tension and compression tests; in either case, the effect is more pronounced as the plastic strain increases. The degraded response of concrete is characterized by two independent uniaxial damage variables, dt and dc, which are assumed to be functions of the plastic strains, temperature, and field variables:

(5)dt=dt(ε~tpl,θ,fi),(0dt1),dc=dc(ε~cpl,θ,fi),(0dc1).
Figure 1. Response of concrete to uniaxial loading in tension (a) and compression (b).

The uniaxial degradation variables are increasing functions of the equivalent plastic strains. They can take values ranging from zero, for the undamaged material, to one, for the fully damaged material.

If E0 is the initial (undamaged) elastic stiffness of the material, the stress-strain relations under uniaxial tension and compression loading are, respectively:

σt=(1-dt)E0(εt-ε~tpl),σc=(1-dc)E0(εc-ε~cpl).

Under uniaxial loading cracks propagate in a direction transverse to the stress direction. The nucleation and propagation of cracks, therefore, causes a reduction of the available load-carrying area, which in turn leads to an increase in the effective stress. The effect is less pronounced under compressive loading since cracks run parallel to the loading direction; however, after a significant amount of crushing, the effective load-carrying area is also significantly reduced. The effective uniaxial cohesion stresses, σ¯t and σ¯c , are given as

σ¯t=σt(1-dt)=E0(εt-ε~tpl),σ¯c=σc(1-dc)=E0(εc-ε~cpl).

The effective uniaxial cohesion stresses determine the size of the yield (or failure) surface.

Uniaxial cyclic conditions

Under uniaxial cyclic loading conditions the degradation mechanisms are quite complex, involving the opening and closing of previously formed micro-cracks, as well as their interaction. Experimentally, it is observed that there is some recovery of the elastic stiffness as the load changes sign during a uniaxial cyclic test. The stiffness recovery effect, also known as the “unilateral effect,” is an important aspect of the concrete behavior under cyclic loading. The effect is usually more pronounced as the load changes from tension to compression, causing tensile cracks to close, which results in the recovery of the compressive stiffness.

The concrete damaged plasticity model assumes that the reduction of the elastic modulus is given in terms of a scalar degradation variable, d, as

E=(1-d)E0,

where E0 is the initial (undamaged) modulus of the material.

This expression holds both in the tensile (σ11>0) and compressive (σ11<0) sides of the cycle. The stiffness reduction variable, d, is a function of the stress state and the uniaxial damage variables, dt and dc. For the uniaxial cyclic conditions, Abaqus assumes that

(6)(1-d)=(1-stdc)(1-scdt),   0st,    sc1,

where st and sc are functions of the stress state that are introduced to represent stiffness recovery effects associated with stress reversals. They are defined according to

st=1-wtr*(σ¯11);   0wt1,sc=1-wc(1-r*(σ¯11));    0wc1,

where

r*(σ¯11)=H(σ¯11)={1ifσ¯11>00ifσ¯11<0.

The weight factors wt and wc, which are assumed to be material properties, control the recovery of the tensile and compressive stiffness upon load reversal. To illustrate this, consider the example in Figure 2, where the load changes from tension to compression. Assume that there was no previous compressive damage (crushing) in the material; that is, ε~cpl=0 and dc=0. Then

(1-d)=(1-scdt)=(1-(1-wc(1-r*))dt).

In tension (σ¯11>0), r*=1; thus, d=dt as expected. In compression (σ¯11<0), r*=0, and d=(1-wc)dt. If wc=1, then d=0; therefore, the material fully recovers the compressive stiffness (which in this case is the initial undamaged stiffness, E=E0). If, on the other hand, wc=0, then d=dt and there is no stiffness recovery. Intermediate values of wc result in partial recovery of the stiffness.

Figure 2. Illustration of the effect of the compression stiffness recovery parameter wc.

The evolution equations of the equivalent plastic strains are also generalized to the uniaxial cyclic conditions as

(7)ε~˙tpl=r*ε˙11pl,ε~˙cpl=-(1-r*)ε˙11pl,

which clearly reduces to Equation 4 during the tensile and compressive phases of the cycle.

Multiaxial conditions

The evolution equations for the hardening variables must be extended for the general multiaxial conditions. Based on Lee and Fenves (1998) we assume that the equivalent plastic strain rates are evaluated according to the expressions

(8)ε~˙tpl=defr(σ¯^)ε˙^maxpl,ε~˙cpl=def-(1-r(σ¯^))ε˙^minpl,

where ε˙^maxpl and ε˙^minpl are, respectively, the maximum and minimum eigenvalues of the plastic strain rate tensor ε˙pl and

r(σ¯^)=defi=13σ¯^ii=13|σ¯^i|;   0r(σ¯^)1

is a stress weight factor that is equal to one if all principal stresses σ¯^i,(i=1,2,3), are positive and equal to zero if they are negative. The Macauley bracket is defined by x=12(|x|+x). In uniaxial loading conditions Equation 8 reduces to the uniaxial definitions Equation 4 and Equation 7, since ε˙^maxpl=ε˙11pl in tension, and ε˙^minpl=ε˙11pl in compression.

If the eigenvalues of the plastic strain rate tensor (ε˙^i,    i=1,2,3) are ordered such that ε˙^maxpl=ε˙^1ε˙^2ε˙^3=ε˙^minpl, the evolution equation for general multiaxial stress conditions can be expressed in the following matrix form:

ε~˙pl=[ε~˙tplε~˙cpl]=h^(σ¯^,ε~pl)ε˙^pl,

where

h^(σ¯^,ε~pl)=[r(σ¯^)0000-(1-r(σ¯^))],

and

ε˙^pl=[ε˙^1ε˙^2ε˙^3].

Elastic stiffness degradation

The plastic-damage concrete model assumes that the elastic stiffness degradation is isotropic and characterized by a single scalar variable, d:

(9)Del=(1-d)D0el;   0d1.

The definition of the scalar degradation variable d must be consistent with the uniaxial monotonic responses (dt and dc), and it should also should capture the complexity associated with the degradation mechanisms under cyclic loading. For the general multiaxial stress conditions Abaqus assumes that

(10)(1-d)=(1-stdc)(1-scdt),   0st,    sc1,

similar to the uniaxial cyclic case, only that st and sc are now given in terms of the function r(σ¯^) as

st=1-wtr(σ¯^);   0wt1,sc=1-wc(1-r(σ¯^));    0wc1.

It can be easily verified that Equation 10 for the scalar degradation variable is consistent with the uniaxial response.

The experimental observation in most quasi-brittle materials, including concrete, is that the compressive stiffness is recovered upon crack closure as the load changes from tension to compression. On the other hand, the tensile stiffness is not recovered as the load changes from compression to tension once crushing micro-cracks have developed. This behavior, which corresponds to wt=0 and wc=1, is the default used by Abaqus. Figure 3 illustrates a uniaxial load cycle assuming the default behavior.

Figure 3. Uniaxial load cycle (tension-compression-tension) assuming default values for the stiffness recovery factors: wt=0 and wc=1.

Yield condition

The plastic-damage concrete model uses a yield condition based on the yield function proposed by Lubliner et al. (1989) and incorporates the modifications proposed by Lee and Fenves (1998) to account for different evolution of strength under tension and compression. In terms of effective stresses the yield function takes the form

(11)F(σ¯,ε~pl)=11-α(q¯-3αp¯+β(ε~pl)σ¯^max-γ-σ¯^max)-σ¯c(ε~cpl)0,

where α and γ are dimensionless material constants;

p¯=-13σ¯:I

is the effective hydrostatic pressure;

q¯=32S¯:S¯

is the Mises equivalent effective stress;

S¯=p¯I+σ¯,

is the deviatoric part of the effective stress tensor σ¯; and σ¯^max is the algebraically maximum eigenvalue of σ¯. The function β(ε~pl) is given as

β(ε~pl)=σ¯c(ε~cpl)σ¯t(ε~tpl)(1-α)-(1+α),

where σ¯t and σ¯c are the effective tensile and compressive cohesion stresses, respectively.

In biaxial compression, with σ¯^max=0, Equation 11 reduces to the well-known Drucker-Prager yield condition. The coefficient α can be determined from the initial equibiaxial and uniaxial compressive yield stress, σb0 and σc0, as

α=σb0-σc02σb0-σc0.

Typical experimental values of the ratio σb0/σc0 for concrete are in the range from 1.10 to 1.16, yielding values of α between 0.08 and 0.12 (Lubliner et al., 1989).

The coefficient γ enters the yield function only for stress states of triaxial compression, when σ¯^max<0. This coefficient can be determined by comparing the yield conditions along the tensile and compressive meridians. By definition, the tensile meridian (TM) is the locus of stress states satisfying the condition σ¯^max=σ¯^1>σ¯^2=σ¯^3 and the compressive meridian (CM) is the locus of stress states such that σ¯^max=σ¯^1=σ¯^2>σ¯^3, where σ¯^1, σ¯^2, and σ¯^3 are the eigenvalues of the effective stress tensor. It can be easily shown that (σ¯^max)TM=23q¯-p¯ and (σ¯^max)CM=13q¯-p¯, along the tensile and compressive meridians, respectively. With σ¯^max<0 the corresponding yield conditions are

(23γ+1)q¯-(γ+3α)p¯=(1-α)σ¯c,        (TM)
(13γ+1)q¯-(γ+3α)p¯=(1-α)σ¯c.                    (CM)

Let Kc=q¯(TM)/q¯(CM) for any given value of the hydrostatic pressure p¯ with σ¯^max<0; then

Kc=γ+32γ+3.

The fact that Kc is constant does not seem to be contradicted by experimental evidence (Lubliner et al., 1989). The coefficient γ is, therefore, evaluated as

γ=3(1-Kc)2Kc-1.

A value of Kc=23, which is typical for concrete, gives γ=3

If σ¯^max>0, the yield conditions along the tensile and compressive meridians reduce to

(23β+1)q¯-(β+3α)p¯=(1-α)σc¯,(TM)
(13β+1)q¯-(β+3α)p¯=(1-α)σ¯c.                    (CM)

Let Kt=q¯(TM)/q¯(CM) for any given value of the hydrostatic pressure p¯ with σ¯^max>0; then

Kt=β+32β+3.

Typical yield surfaces are shown in Figure 4 in the deviatoric plane and in Figure 5 for plane-stress conditions.

Figure 4. Yield surfaces in the deviatoric plane, corresponding to different values of Kc.

Figure 5. Yield surface in plane stress.

Flow rule

The plastic-damage model assumes nonassociated potential flow,

ε˙pl=λ˙G(σ¯)σ¯.

The flow potential G chosen for this model is the Drucker-Prager hyperbolic function:

(12)G=(ϵσt0tanψ)2+q¯2-p¯tanψ,

where ψ is the dilation angle measured in the p–q plane at high confining pressure; σt0 is the uniaxial tensile stress at failure; and ϵ is a parameter, referred to as the eccentricity, that defines the rate at which the function approaches the asymptote (the flow potential tends to a straight line as the eccentricity tends to zero). This flow potential, which is continuous and smooth, ensures that the flow direction is defined uniquely. The function asymptotically approaches the linear Drucker-Prager flow potential at high confining pressure stress and intersects the hydrostatic pressure axis at 90°. See Models for granular or polymer behavior for further discussion of this potential.

Because plastic flow is nonassociated, the use of the plastic-damage concrete model requires the solution of nonsymmetric equations.

Viscoplastic regularization

Material models exhibiting softening behavior and stiffness degradation often lead to severe convergence difficulties in implicit analysis programs. Some of these convergence difficulties can be overcome by using a viscoplastic regularization of the constitutive equations. The concrete damaged plasticity model can be regularized using viscoplasticity, therefore permitting stresses to be outside of the yield surface. We use a generalization of the Duvaut-Lions regularization, according to which the viscoplastic strain rate tensor, ε˙vpl, is defined as

(13)ε˙vpl=1μ(εpl-εvpl).

Here μ is the viscosity parameter representing the relaxation time of the viscoplastic system and εpl is the plastic strain evaluated in the inviscid backbone model.

Similarly, a viscous stiffness degradation variable, dv, for the viscoplastic system is defined as

(14)d˙v=1μ(d-dv),

where d is the degradation variable evaluated in the inviscid backbone model. The stress-strain relation of the viscoplastic model is given as

(15)σ=(1-dv)D0el:(ε-εvpl).

The solution of the viscoplastic system relaxes to that of the inviscid case as t/μ, where t represents time. Using the viscoplastic regularization with a small value for the viscosity parameter (small compared to the characteristic time increment) usually helps improve the rate of convergence of the model in the softening regime, without compromising results.

Integration of the model

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.