MultiPhysicsVault/.raw/AbaqusTheoriesManual/AbaqusTheoriesManual_061.md

The linear part of the behavior may be compressible or incompressible depending on the value of Poisson's ratio, but the nonlinear part of the behavior is incompressible (because the flow is normal to the Mises stress potential). Since the model will generally be used for cases when the deformation is dominated by plastic flow, the use of selectively reduced integration elements or "hybrid" (mixed formulation) elements is recommended with this material model (except in plane stress).

Strain energy density

The model is often used to obtain fully plastic solutions for fracture mechanics when the J -integral is needed. For evaluation of the J -integral the strain energy density is required. This is

W = \int \pmb {\sigma}: d \pmb {\varepsilon}.

From Equation 4.3.9-2 this may be obtained as

W = \frac {1 + \nu}{3 E} q ^ {2} + \frac {3}{2} \frac {1 - 2 \nu}{E} p ^ {2} + \frac {n}{n + 1} \frac {\alpha}{E \sigma_ {0} ^ {n - 1}} q ^ {n + 1}.

Stress solution

During an analysis at each integration point the latest estimate of the kinematic solution is provided to the constitutive routines, which must provide the corresponding stress tensor calculated for the material model being used. Since this material model is nonlinear, we solve for the stresses by using the methods described below.

The uniaxial case (the only nonzero stress component is one direct stress)

In this case

Equation 4.3.9-3

E \varepsilon = \sigma + \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} \sigma ,

where now q = \left| \sigma \right| .

We solve Equation 4.3.9-3 for \sigma using Newton's method. Writing c _ { \sigma } as the correction to \sigma , the Newton equations for Equation 4.3.9-3 are

\left[ 1 + n \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} \right] c _ {\sigma} = E \varepsilon - \sigma - \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} \sigma ,

\sigma = \sigma + c _ {\sigma}.

Mechanical Constitutive Theories

As an initial guess we use \sigma = E | \varepsilon | { \mathrm { ~ i f ~ } } E | \varepsilon | \leq \sigma _ { 0 } and \sigma = \pm { [ E | \varepsilon | \sigma _ { 0 } ^ { n - 1 } / \alpha ] } ^ { 1 / n } (with the sign chosen as the same sign as ") if E | \varepsilon | > \sigma _ { 0 } .

In this case the material stiffness matrix is

\frac {\partial \sigma}{\partial \varepsilon} = \frac {E}{1 + n \alpha (q / \sigma_ {0}) ^ {n - 1}}.

All strain components defined by the kinematic solution

When all strain components are defined kinematically (that is, all cases except uniaxial and plane stress), projecting Equation 4.3.9-2 onto I gives

p = - \frac {E}{3 (1 - 2 \nu)} \mathbf {I}: \varepsilon ,

a linear relationship for the volumetric behavior. Defining the deviatoric strain as

\mathbf {e} = \varepsilon - \frac {1}{3} \mathbf {I I}: \varepsilon

and using Equation 4.3.9-2, we find

E \mathbf {e} = \left(1 + \nu + \frac {3}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) \mathbf {S}.

Equation 4.3.9-4

Defining the equivalent deviatoric strain as

\bar {e} = \sqrt {\frac {2}{3} \mathbf {e} : \mathbf {e}}

Equation 4.3.9-5

and using Equation 4.3.9-4, we obtain the scalar equation

E \bar {e} = \frac {2}{3} (1 + \nu) q + \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} q.

Equation 4.3.9-6

This equation is solved for q using Newton's method:

\left[ \frac {2}{3} (1 + \nu) + n \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} \right] c _ {q} = E \bar {e} - \frac {2}{3} (1 + \nu) q - \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} q,

q = q + c _ {q}.

As for the uniaxial case we will use the starting guesses:

q = \frac {3}{2} \frac {E}{1 + \nu} \bar {e} \qquad \mathrm{if} \qquad \frac {3}{2} \frac {E}{1 + \nu} \bar {e} \leq \sigma_ {0},

and

q = \left(\frac {\sigma_ {0} ^ {n - 1} E \bar {e}}{\alpha}\right) ^ {1 / n} \qquad \mathrm{if} \qquad \frac {3}{2} \frac {E}{1 + \nu} \bar {e} > \sigma_ {0}.

Equation 4.3.9-6 can also be used to define

1 + \nu + \frac {3}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} = \frac {3}{2} E \frac {\bar {e}}{q},

so that Equation 4.3.9-4 becomes

Equation 4.3.9-7

\mathbf {S} = \frac {2}{3} \frac {q}{\bar {e}} \mathbf {e}.

Thus, once q is known, S is defined; and, hence, ¾ is known as

\pmb {\sigma} = \mathbf {S} - p \mathbf {I}.

The material stiffness is defined as follows. From Equation 4.3.9-7 we have

\partial \mathbf {S} = \frac {2}{3} \frac {q}{\bar {e}} \left[ \partial \mathbf {e} + \left(\frac {\partial q}{q} - \frac {\partial \bar {e}}{\bar {e}}\right) \mathbf {e} \right];

Equation 4.3.9-8

and Equation 4.3.9-6 gives

\partial q = \frac {E}{\frac {2}{3} (1 + \nu) + n \alpha (q / \sigma_ {0}) ^ {n - 1}} \partial \bar {e};

and from Equation 4.3.9-5,

\partial \bar {e} = \frac {2}{3} \frac {\mathbf {e}}{\bar {e}}: \partial \mathbf {e}.

Using these results, Equation 4.3.9-8 becomes

\partial \mathbf {S} = \frac {2}{3} \frac {q}{\bar {e}} \left[ \mathfrak {I} - \frac {2}{3 \bar {e}} \left\{\frac {1}{\bar {e}} - \frac {1}{\bar {e} + (n - 1) \alpha (q / \sigma_ {0}) ^ {n - 1} (q / E)} \right\} \mathbf {e} \mathbf {e} \right]: \partial \mathbf {e}.

Now

\partial \pmb {\sigma} = \partial \mathbf {S} - \partial p \mathbf {I},

\partial p = - \frac {E}{3 (1 - 2 \nu)} \mathbf {I}: \partial \pmb {\varepsilon},

and

\partial \mathbf {e} = \left(\mathfrak {I} - \frac {1}{3} \mathbf {I I}\right): \partial \pmb {\varepsilon}.

Combining these results, we obtain the material stiffness as

\partial \pmb {\sigma} = \left[ \frac {2}{3} \frac {q}{\bar {e}} \left(\mathfrak {F} - \frac {2}{3 \bar {e}} \left\{\frac {1}{\bar {e}} - \frac {1}{\bar {e} + (n - 1) \alpha \left(q / \sigma_ {0}\right) ^ {n - 1} (q / E)} \right\} \mathbf {e} \mathbf {e}\right) + \frac {1}{3} \left(\frac {E}{1 - 2 \nu} - \frac {2}{3} \frac {q}{\bar {e}}\right) \mathbf {I} \mathbf {I} \right]: \partial \pmb {\varepsilon}.

Plane stress

For this case we obtain a solution for ¾ assuming that the material is fully incompressible. We then use this solution as a starting guess for a Newton loop to find ¾. Finally, we compute the corresponding material stiffness matrix.

Incompressible approximation:

In this case \nu = 1 / 2 and \varepsilon _ { 3 3 } = - ( \varepsilon _ { 1 1 } + \varepsilon _ { 2 2 } ) , so that e¹ is known and Equation 4.3.9-6 can be solved as before for q . Equation 4.3.9-7 then defines S, and the plane stress constraint requires that \sigma _ { 3 3 } = 0 , so

p = S _ {3 3} = - \frac {2}{3} \frac {q}{\bar {e}} (\varepsilon_ {1 1} + \varepsilon_ {2 2});

and, hence, we obtain the initial estimate of the solution

\sigma_ {1 1} = \frac {2}{3} \frac {q}{\bar {e}} (2 \varepsilon_ {1 1} + \varepsilon_ {2 2});

\sigma_ {2 2} = \frac {2}{3} \frac {q}{\bar {e}} (\varepsilon_ {1 1} + 2 \varepsilon_ {2 2}),

and

\sigma_ {1 2} = \frac {q}{3 \bar {e}} \gamma_ {1 2}.

Newton solution for the actual stresses:

Equation 4.3.9-2 defines the stresses, where, for this case,

p = - \frac {1}{3} (\sigma_ {1 1} + \sigma_ {2 2}),

\mathbf {S} = \boldsymbol {\sigma} + p \mathbf {I},

and

q = \sqrt {\frac {3}{2} \mathbf {S} : \mathbf {S}}.

Therefore, Equation 4.3.9-2 becomes

Equation 4.3.9-9

E \pmb {\varepsilon} = \left(1 + \nu + \frac {3}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) \pmb {\sigma} + 3 \left(\nu + \frac {1}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) p \mathbf {I}.

These equations are solved using Newton's method:

\left[ \left(1 + \nu + \frac {3}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) \mathfrak {F} + \frac {9}{4} \alpha (n - 1) \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} \frac {\mathbf {S}}{q} \frac {\mathbf {S}}{q} - \left(\nu + \frac {1}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) \mathbf {I I} \right]: \mathbf {c} _ {\sigma}

= E \pmb {\varepsilon} - \left(1 + \nu + \frac {3}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) \pmb {\sigma} - 3 \left(\nu + \frac {1}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) p \mathbf {I},

\boldsymbol {\sigma} = \boldsymbol {\sigma} + \mathbf {c} _ {\sigma}.

The material stiffness is directly available from Equation 4.3.9-9 as

E \left[ \left(1 + \nu + \frac {3}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) \mathfrak {S} + \frac {9}{4} \alpha (n - 1) \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1} \frac {\mathbf {S}}{q} \frac {\mathbf {S}}{q} - \left(\nu + \frac {1}{2} \alpha \left(\frac {q}{\sigma_ {0}}\right) ^ {n - 1}\right) \mathbf {I I} \right] ^ {- 1}.

4.3.10 Heat generation caused by plastic straining

The *INELASTIC HEAT FRACTION option allows the introduction of a factor, ´, which defines heat generation caused by mechanical dissipation associated with plastic straining. This term can be introduced as a source of coupling for thermal-mechanical analysis. Such coupling might be important

in a simulation in which extensive inelastic deformation is occurring fairly rapidly in a material whose mechanical properties are temperature dependent. If the process is very slow, the heat generated by the plastic deformation has time to dissipate, and uncoupled, isothermal, analysis is sufficient to model the process. If the process is extremely rapid, the heat has no time to diffuse, and uncoupled, adiabatic, analysis (in which each integration point is treated as if it is thermally insulated from its neighbors) is sufficient. Fully coupled analysis is required for cases that lie far enough from both extremes. This section defines the heat generation term caused by inelastic straining and describes how this term contributes to the overall Newton solution scheme.

The model assumes that plastic straining gives rise to a heat flux per unit volume of

r ^ {p l} = \eta \pmb {\sigma}: \dot {\pmb {\varepsilon}} ^ {p l},

where r ^ { p l } is the heat flux that is added into the thermal energy balance, ´ is the factor defined on the *INELASTIC HEAT FRACTION option and is assumed to be a constant, ¾ is the stress, and \dot { \varepsilon } ^ { p l } is the rate of plastic straining. For all of the plasticity models in ABAQUS, the plastic strain increment is written from the flow potential as

\dot {\pmb {\varepsilon}} ^ {p l} = \dot {\varepsilon} ^ {p l} \mathbf {n},

where n is the flow direction (we assume { \bf n } = { \bf n } ( \pmb { \sigma } , \varepsilon ^ { p l } , \theta ) , where µ is the temperature) and \varepsilon ^ { p l } is a scalar measure of plastic straining that is used as a hardening parameter in the yield surface and flow potential definitions in some of the plasticity models. In this case we consider isotropic hardening theories only: ABAQUS provides only thermo-mechanical coupling for such models.

ABAQUS generally uses a backward Euler scheme to integrate the plastic strain, so r ^ { p l } at the end of the increment is approximated as

r ^ {p l} = \frac {1}{2 \Delta t} \eta \Delta \varepsilon^ {p l} \mathbf {n}: (\pmb {\sigma} + \pmb {\sigma} _ {t}),

where all quantities are evaluated at the end of the increment (at time t + \Delta t ) except \pmb { \sigma } _ { t } . This notation is adopted throughout the remainder of this section. This term is used as the contribution to the thermal energy balance equation.

When Newton's method is used to solve the nonlinear equations, the coupling term gives rise to three contributions to the Jacobian matrix for the Newton method:

\frac {\partial r ^ {p l}}{\partial \theta}, \frac {\partial r ^ {p l}}{\partial \varepsilon}

from the thermal energy balance equation, and

\frac {\partial \pmb {\sigma}}{\partial \theta}

Mechanical Constitutive Theories

from the mechanical equilibrium equation. The general form for these terms is now derived.

The mechanical constitutive model has the following general form. The elasticity defines the stresses by

Equation 4.3.10-1

\pmb {\sigma} = \frac {\partial W}{\partial \pmb {\varepsilon} ^ {e l}},

where W = W ( \pmb { \varepsilon } ^ { e l } , \theta ) is the strain energy density potential and \pmb \varepsilon ^ { e l } is the mechanical elastic strain. We implicitly assume that the elasticity is not fully incompressible, although the derivation is not significantly different if this is not the case, since the pressure stress will do no work in a fully incompressible material and so makes no contribution to the terms under discussion.

We assume that there is an additive strain rate decomposition that can be integrated to give

Equation 4.3.10-2

\varepsilon^ {e l} = \varepsilon - \varepsilon^ {t h} - \varepsilon^ {p l},

where " is the total strain and \varepsilon ^ { t h } is the strain caused by thermal expansion. In the constitutive models in ABAQUS \pmb { \varepsilon } ^ { t h } = \pmb { \varepsilon } ^ { t h } ( \theta ) only. This form of decomposition of the deformation depends on " being measured as the integrated rate of deformation and on the elastic and thermal strains being small: this is true for the standard plasticity models provided in the program.

The plastic flow definition is integrated by the backward Euler method to give

Equation 4.3.10-3

\Delta \pmb {\varepsilon} ^ {p l} = \Delta \varepsilon^ {p l} \mathbf {n}.

Finally, assuming there is a single active yield surface or a single active flow surface, rate-independent models introduce a yield surface constraint, while rate-dependent models provide an integrated flow rate constraint, both of which are incorporated in the general form

Equation 4.3.10-4

f (\boldsymbol {\sigma}) = \overline {{\sigma}},

where f ( \pmb { \sigma } ) is a scalar stress function (for example, the Mises or Hill stress function for simple metal plasticity models) and \overline { { \sigma } } = \sigma ^ { 0 } ( \varepsilon ^ { p l } , \theta ) is the yield stress for a rate-independent model, while \overline { { \sigma } } = B \sigma ^ { 0 } for a rate-dependent model, where

B = B (\frac {\Delta \varepsilon^ {p l}}{\Delta t}, \theta)

defines the rate effect from the average plastic strain rate over the increment. For example, by default, the *RATE DEPENDENT suboption defines

\dot {\varepsilon} ^ {p l} = D \left(\frac {\tilde {q}}{\sigma^ {0}} - 1\right) ^ {p},

where \tilde { q } is the Mises or Hill equivalent stress, and D(µ) and p ( \theta ) are material parameters. Using the average plastic strain rate over the increment in this expression defines

B = 1 + \left(\frac {\Delta \varepsilon^ {p l}}{D \Delta t}\right) ^ {1 / p}.

Equation 4.3.10-1 to Equation 4.3.10-4 are a general definition of all of the standard isotropic hardening plasticity models integrated by the backward Euler method.

We now take variations of these equations with respect to all quantities at the end of the increment:

\partial \pmb {\sigma} = \frac {\partial^ {2} W}{\partial \pmb {\varepsilon} ^ {e l} \partial \theta} \partial \theta + \frac {\partial^ {2} W}{\partial \pmb {\varepsilon} ^ {e l} \partial \pmb {\varepsilon} ^ {e l}}: \left(\partial \pmb {\varepsilon} - \frac {\partial \pmb {\varepsilon} ^ {t h}}{\partial \theta} \partial \theta - \partial \pmb {\varepsilon} ^ {p l}\right),

\partial \pmb {\varepsilon} ^ {p l} = \partial \varepsilon^ {p l} \mathbf {n} + \Delta \varepsilon^ {p l} \left(\frac {\partial \mathbf {n}}{\partial \pmb {\sigma}}: \partial \pmb {\sigma} + \frac {\partial \mathbf {n}}{\partial \varepsilon^ {p l}} \partial \varepsilon^ {p l} + \frac {\partial \mathbf {n}}{\partial \theta} \partial \theta\right),

and

\mathbf {m}: \partial \pmb {\sigma} = \frac {\partial \overline {{\sigma}}}{\partial \varepsilon^ {p l}} \partial \varepsilon^ {p l} + \frac {\partial \overline {{\sigma}}}{\partial \theta} \partial \theta ,

where

\mathbf {m} = \frac {\partial f}{\partial \pmb {\sigma}}.

For simplicity of notation we now define

\mathbf {D} ^ {e l} = \frac {\partial^ {2} W}{\partial \pmb {\varepsilon} ^ {e l} \partial \pmb {\varepsilon} ^ {e l}},

\hat {\mathbf {n}} = \mathbf {n} + \Delta \varepsilon^ {p l} \frac {\partial \mathbf {n}}{\partial \varepsilon^ {p l}},

d = \mathbf {m}: \mathbf {D} ^ {e l}: \hat {\mathbf {n}} + \frac {\partial \overline {{\boldsymbol {\sigma}}}}{\partial \varepsilon^ {p l}},

\mathbf {Z} = \mathbf {I} - \frac {1}{d} \hat {\mathbf {n}} \mathbf {m}: \mathbf {D} ^ {e l},

\mathbf {H} = \mathbf {I} + \Delta \varepsilon^ {p l} \mathbf {D} ^ {e l}: \mathbf {Z}: \frac {\partial \mathbf {n}}{\partial \pmb {\sigma}},

\mathbf {g} = \frac {\partial^ {2} W}{\partial \pmb {\varepsilon} ^ {e l} \partial \theta} - \mathbf {D} ^ {e l}: \left(\frac {\partial \pmb {\varepsilon} ^ {t h}}{\partial \theta} + \Delta \varepsilon^ {p l} \frac {\partial \mathbf {n}}{\partial \theta}\right),

\tilde {\mathbf {n}} = \mathbf {n} + (\pmb {\sigma} + \pmb {\sigma} _ {t}): \frac {\partial \mathbf {n}}{\partial \pmb {\sigma}},

\tilde {g} = \mathbf {m}: \mathbf {g} - \frac {\partial \overline {{\sigma}}}{\partial \theta},

\mathbf {h} = \mathbf {H} ^ {- 1}: \left\{\mathbf {g} - \frac {\tilde {g}}{d} \mathbf {D} ^ {e l}: \mathbf {n} \right\},

and

\mathbf {D} = \mathbf {H} ^ {- 1}: \mathbf {D} ^ {e l}: \mathbf {Z}.

These expressions allow us to write

\partial \boldsymbol {\sigma} = \mathbf {D}: \partial \boldsymbol {\varepsilon} + \mathbf {h} \partial \theta ,

and

\begin{array}{l} \partial r ^ {p l} = \frac {1}{2 \Delta t} \eta \left\{\Delta \varepsilon^ {p l} \tilde {\mathbf {n}}: \mathbf {D} + \frac {1}{d} (\pmb {\sigma} + \pmb {\sigma} _ {t}): \hat {\mathbf {n}} \mathbf {m}: \mathbf {D} ^ {e l}: \left[ \mathbf {I} - \Delta \varepsilon^ {p l} \frac {\partial \mathbf {n}}{\partial \pmb {\sigma}}: \mathbf {D} \right] \right\}: \partial \pmb {\varepsilon} \\ + \frac {1}{2 \Delta t} \eta \bigg \{\frac {1}{d} (\pmb {\sigma} + \pmb {\sigma} _ {t}): \hat {\mathbf {n}} \left(\tilde {g} - \Delta \varepsilon^ {p l} \mathbf {m}: \mathbf {D} ^ {e l}: \frac {\partial \mathbf {n}}{\partial \pmb {\sigma}}: \mathbf {h}\right) \\ \left. + \Delta \varepsilon^ {p l} \left(\boldsymbol {\sigma} + \boldsymbol {\sigma} _ {t}\right): \frac {\partial \mathbf {n}}{\partial \theta} + \Delta \varepsilon^ {p l} \left(\mathbf {n} + \left(\boldsymbol {\sigma} + \boldsymbol {\sigma} _ {t}\right): \frac {\partial \mathbf {n}}{\partial \boldsymbol {\sigma}}\right): \mathbf {h} \right\} \partial \theta . \\ \end{array}

4.4 Plasticity for non-metals

4.4.1 Porous elasticity

The porous elasticity model in ABAQUS/Standard is designed to be used in conjunction with plasticity

Mechanical Constitutive Theories

models that allow plastic volume changes as described in ``Models for granular or polymer behavior,'' Section 4.4.2, to ``Models for crushable foams,'' Section 4.4.6. Use of the porous elasticity model without one of these plasticity options is not recommended. The model is based on the experimental observation that in porous materials during elastic (recoverable) straining, the change in the void ratio--e--and the change in the logarithm of the equivalent pressure stress-- p--defined as

p = - \frac {1}{3} \mathrm{trace} \pmb {\sigma} = - \frac {1}{3} \pmb {\sigma}: \mathbf {I},

are linearly related, so that in rate form,

d e ^ {e l} = - \kappa d (\ln (p)),

where ∙ is a material parameter. In this form the material has zero tensile strength. If the tensile strength p _ { t } ^ { e l } is nonzero, the equivalent relation is

Equation 4.4.1-1

d e ^ {e l} = - \kappa d \left(\ln (p + p _ {t} ^ {e l})\right),

which includes the special case of zero tensile strength ( p _ { t } ^ { e l } = 0 ) . It can be shown that, if the compressibility of the solid material is neglected, the volume change of a material sample is

Equation 4.4.1-2

J = \frac {1 + e}{1 + e _ {0}},

where e _ { 0 } is the initial void ratio. If we define the elastic void ratio from the elastic volume change according to the relationship

Equation 4.4.1-3

J ^ {e l} = \frac {1 + e ^ {e l}}{1 + e _ {0}}

and then integrate the linear relation, the volumetric elasticity relationship is

Equation 4.4.1-4

\frac {\kappa}{1 + e _ {0}} \ln \left(\frac {p + p _ {t} ^ {e l}}{p _ {0} + p _ {t} ^ {e l}}\right) = 1 - J ^ {e l},

where p _ { 0 } is the initial pressure stress, prescribed by initial conditions. Note that for a zero tensile strength material it is required that p _ { 0 } > 0 . . This equation can be inverted to yield

Equation 4.4.1-5

p = - p _ {t} ^ {e l} + \left(p _ {0} + p _ {t} ^ {e l}\right) \exp \left[ \frac {1 + e _ {0}}{\kappa} (1 - \exp \varepsilon_ {v o l} ^ {e l}) \right],

as illustrated in Figure 4.4.1-1.