김경종
17 KiB
Isotropic hardening is generally considered to be a suitable model for problems in which the plastic straining goes well beyond the incipient yield state where the Bauschinger effect is noticeable ( Rice, 1975). Therefore, this hardening theory is used for such applications as dynamic problems involving finite strains and manufacturing processes--any process involving large plastic strain and in which the plastic strain does not continuously reverse direction sharply.
Some cases, such as low-cycle fatigue situations, involve relatively low amplitude strain cycling. In these cases it becomes important to model the Bauschinger effect. Kinematic hardening is the simplest theory that does this. ABAQUS offers a linear kinematic and a nonlinear isotropic/kinematic hardening model for such cases. These models are described in ``Models for metals subjected to cyclic loading,'' Section 4.3.5.
4.3.2 Isotropic elasto-plasticity
This material model is very commonly used for metal plasticity calculations, either as a rate-dependent or as a rate-independent model, and has a particularly simple form. Because of this simplicity the algebraic equations associated with integrating the model are easily developed in terms of a single variable, and the material stiffness matrix can be written explicitly. This results in particularly efficient code. In this section these equations are developed.
For simplicity of notation all quantities not explicitly associated with a time point are assumed to be evaluated at the end of the increment.
The Mises yield function with associated flow means that there is no volumetric plastic strain; since the elastic bulk modulus is quite large, the volume change will be small. Thus, we can define the volume strain as
\varepsilon_ {v o l} = \mathrm{trace} (\pmb {\varepsilon});
and, hence, the deviatoric strain is
\mathbf {e} = \pmb {\varepsilon} - \frac {1}{3} \varepsilon_ {v o l} \mathbf {I}.
Material model definition
The strain rate decomposition is
d \varepsilon = d \varepsilon^ {e l} + d \varepsilon^ {p l}.
Using the standard definition of corotational measures, this can be written in integrated form as
Equation 4.3.2-1
\pmb {\varepsilon} = \pmb {\varepsilon} ^ {e l} + \pmb {\varepsilon} ^ {p l}.
The elasticity is linear and isotropic and, therefore, can be written in terms of two
Mechanical Constitutive Theories
temperature-dependent material parameters. For the purpose of this development it is most appropriate to choose these parameters as the bulk modulus, K, and the shear modulus, G. These are computed readily from the user's input of Young's modulus, E, and Poisson's ratio, º, as
K = \frac {E}{3 (1 - 2 \nu)}
and
G = \frac {E}{2 (1 + \nu)}.
Then, the elasticity can be written in volumetric and deviatoric components as follows.
Volumetric:
Equation 4.3.2-2
p = - K \varepsilon_ {v o l},
where
p = - \frac {1}{3} \operatorname{trace} (\pmb {\sigma})
is the equivalent pressure stress.
Deviatoric:
Equation 4.3.2-3
\mathbf {S} = 2 G \mathbf {e} ^ {e l},
where S is the deviatoric stress,
\mathbf {S} = \boldsymbol {\sigma} + p \mathbf {I}.
The flow rule is
Equation 4.3.2-4
d \mathbf {e} ^ {p l} = d \bar {e} ^ {p l} \mathbf {n},
where
\mathbf {n} = \frac {3}{2} \frac {\mathbf {S}}{q},
q = \sqrt {\frac {3}{2} \mathbf {S} : \mathbf {S}},
and d \bar { e } ^ { p l } is the (scalar) equivalent plastic strain rate.
The plasticity requires that the material satisfy a uniaxial-stress plastic-strain strain-rate relationship. If the material is rate independent, this is the yield condition:
Equation 4.3.2-5
q = \sigma^ {0},
where \sigma ^ { 0 } ( \bar { e } ^ { p l } , \theta ) is the yield stress and is defined by the user as a function of equivalent plastic strain ( \bar { e } ^ { p l } ) and temperature (µ).
If the material is rate dependent, the relationship is the uniaxial flow rate definition:
\dot {\bar {e}} ^ {p l} = h (q, \bar {e} ^ {p l}, \theta),
where h is a known function. For example, the *RATE DEPENDENT option offers an overstress power law model of the form
\dot {\bar {e}} ^ {p l} = D \left(\frac {q}{\sigma^ {0}} - 1\right) ^ {p},
where D ( \theta ) and p ( \theta ) are temperature-dependent material parameters that are defined on the *RATE DEPENDENT option and \sigma ^ { 0 } ( \bar { e } ^ { p l } , \theta ) is the static yield stress.
Integrating this relationship by the backward Euler method gives
Equation 4.3.2-6
\Delta \bar {e} ^ {p l} = \Delta t h (q, \bar {e} ^ {p l}, \theta).
This equation can be inverted (numerically, if necessary) to give q as a function of \bar { e } ^ { p l } at the end of the increment.
Thus, both the rate-independent model and the integrated rate-dependent model give the general uniaxial form
Equation 4.3.2-7
q = \overline {{\sigma}} (\bar {e} ^ {p l}),
where { \overline { { \sigma } } } = \sigma ^ { 0 } for the rate-independent model, and \overline { \sigma } is obtained by inversion of Equation 4.3.2-6 for the rate-dependent model.
Equation 4.3.2-1 to Equation 4.3.2-7 define the material behavior. In any increment when plastic flow is occurring (which is determined by evaluating q based on purely elastic response and finding that its
Mechanical Constitutive Theories
value exceeds ¾0), these equations must be integrated and solved for the state at the end of the increment. As in the general discussion in ``Metal plasticity models,'' Section 4.3.1, the integration is done by applying the backward Euler method to the flow rule (Equation 4.3.2-4), giving
Equation 4.3.2-8
\Delta \mathbf {e} ^ {p l} = \Delta \bar {e} ^ {p l} \mathbf {n}.
Combining this with the deviatoric elasticity (Equation 4.3.2-3) and the integrated strain rate decomposition (Equation 4.3.2-1) gives
Equation 4.3.2-9
\mathbf {S} = 2 G (\mathbf {e} ^ {e l} \big | _ {t} + \Delta \mathbf {e} - \Delta \bar {e} ^ {p l} \mathbf {n}).
Then, using the integrated flow rule (Equation 4.3.2-8), together with the Mises definition of the flow direction, n (in Equation 4.3.2-4), this becomes
(1 + \frac {3 G}{q} \Delta \bar {e} ^ {p l}) \mathbf {S} = 2 G (\mathbf {e} ^ {e l} \big | _ {t} + \Delta \mathbf {e}).
For simplicity of notation we write
\hat {\mathbf {e}} = \left. \mathbf {e} ^ {e l} \right| _ {t} + \Delta \mathbf {e},
so that this equation is
Equation 4.3.2-10
(1 + \frac {3 G}{q} \Delta \bar {e} ^ {p l}) \mathbf {S} = 2 G \hat {\mathbf {e}}.
Taking the inner product of this equation with itself gives
Equation 4.3.2-11
q + 3 G \Delta \bar {e} ^ {p l} = 3 G \tilde {e},
where
\tilde {e} = \sqrt {\frac {2}{3} \hat {\mathbf {e}} : \hat {\mathbf {e}}}.
The Mises equivalent stress, q, must satisfy the uniaxial form defined in Equation 4.3.2-7, so that from Equation 4.3.2-11,
3 G (\tilde {e} - \Delta \bar {e} ^ {p l}) - \overline {{\sigma}} = 0.
Equation 4.3.2-12
This is a nonlinear equation for \Delta \bar { e } ^ { p l } in the general case when \overline { \sigma } depends on the equivalent plastic strain (that is, when the material is rate-dependent, or when there is nonzero work hardening). (It is linear in \Delta e ^ { p l } for rate-independent perfect plasticity.) We solve it by Newton's method:
c ^ {p l} = \frac {3 G (\tilde {e} - \Delta \bar {e} ^ {p l}) - \overline {{\sigma}}}{3 G + H}, \quad \mathrm{where} \quad H = \frac {d \overline {{\sigma}}}{d \bar {e} ^ {p l}},
\Delta \bar {e} ^ {p l} = \Delta \bar {e} ^ {p l} + c ^ {p l},
and we iterate until convergence is achieved.
Once \Delta \bar { e } ^ { p l } is known, the solution is fully defined: Using Equation 4.3.2-5,
q = \overline {{\sigma}};
and so, from Equation 4.3.2-10,
\mathbf {S} = \frac {2 G}{1 + \frac {3 G}{q} \Delta \bar {e} ^ {p l}} \hat {\mathbf {e}}.
From Equation 4.3.2-4,
\mathbf {n} = \frac {3}{2} \frac {\mathbf {S}}{q};
and thus, from Equation 4.3.2-6,
\Delta \mathbf {e} ^ {p l} = \Delta \bar {e} ^ {p l} \mathbf {n}.
For cases where three direct strain components are provided by the kinematic solution (that is, all but plane stress and uniaxial stress cases), Equation 4.3.2-2 defines
p = - K \varepsilon_ {v o l},
so that the solution is then fully defined.
Plane stress
For plane stress \varepsilon _ { 3 3 } is not defined by the kinematics but by the plane stress constraint
\sigma_ {3 3} = 0.
This additional equation (or equivalently p = S _ { 3 3 } ) must be solved along with the yield condition and Equation 4.3.2-9. The predicted third strain component
\hat {\varepsilon} _ {3 3} = - \frac {\nu}{1 - \nu} (\hat {\varepsilon} _ {1 1} + \hat {\varepsilon} _ {2 2}),
where
\hat {\varepsilon} _ {1 1} = \varepsilon_ {1 1} ^ {e l} \big | _ {t} + \Delta \varepsilon_ {1 1}
and
\hat {\varepsilon} _ {2 2} = \varepsilon_ {2 2} ^ {e l} \big | _ {t} + \Delta \varepsilon_ {2 2},
serves as an initial guess toward the final value of \hat { \varepsilon } _ { 3 3 } that enables (with the correct plastic straining) the plane stress condition to be satisfied.
Uniaxial stress
For cases with only one direct strain component defined by the kinematic solution (uniaxial deformation), we require
p = S _ {2 2} = S _ {3 3} (= - \frac {1}{2} S _ {1 1}),
so that
Equation 4.3.2-13
\sigma_ {1 1} = \frac {3}{2} S _ {1 1}.
Material stiffness
For this simple plasticity model the material stiffness matrix can be derived without the need for matrix inversion (as was needed in the general case described in ``Integration of plasticity models,'' Section 4.2.2), as follows.
Taking the variation of Equation 4.3.2-10 with respect to all quantities at the end of the increment gives
Equation 4.3.2-14
(1 + \frac {3 G}{q} \Delta \bar {e} ^ {p l}) \partial \mathbf {S} + \mathbf {S} \frac {3 G}{q} (\partial \bar {e} ^ {p l} - \frac {\Delta \bar {e} ^ {p l}}{q} \partial q) = 2 G \partial \hat {\mathbf {e}}.
Now, from Equation 4.3.2-5,
\partial q = H \partial \bar {e} ^ {p l},
and, from Equation 4.3.2-11,
\partial q + 3 G \partial \bar {e} ^ {p l} = 3 G \partial \tilde {e}.
Combining these last two results,
\partial \bar {e} ^ {p l} = \frac {1}{1 + B} \partial \tilde {e},
where
B = \frac {H}{3 G}.
From the definition of e~ (see Equation 4.3.2-11),
\partial \tilde {e} = \frac {2}{3 \tilde {e}} \hat {\mathbf {e}}: \partial \hat {\mathbf {e}},
and, hence,
\partial \bar {e} ^ {p l} = \frac {2}{3} \frac {1}{\tilde {e} (1 + B)} \hat {\mathbf {e}}: \partial \hat {\mathbf {e}}.
Combining these results with Equation 4.3.2-14 gives
\partial \mathbf {S} = \left[ Q \mathfrak {I} - R \mathbf {S} \mathbf {S} \right]: \partial \hat {\mathbf {e}},
where \begin{array} { r } { Q = \frac { 2 } { 3 } q / { \tilde { e } } . } \end{array} , = is the fourth-order unit tensor, and
R = \frac {1}{q \tilde {e}} \frac {(1 - \Delta \bar {e} ^ {p l} H / q)}{(1 + B)}.
For all cases where three direct strains are defined by the kinematic solution, the material stiffness is completed by
\partial p = - K \mathbf {I}: \partial \varepsilon ,
so since
\pmb {\sigma} = \mathbf {S} - p \mathbf {I}
and
\partial \hat {\mathbf {e}} = \left(\mathfrak {I} - \frac {1}{3} \mathbf {I I}\right): \partial \pmb {\varepsilon},
we have
\partial \pmb {\sigma} = \left[ Q \widetilde {\mathfrak {S}} + \left(K - \frac {1}{3} Q\right) \mathbf {I I} - R \mathbf {S S} \right]: \partial \pmb {\varepsilon}.
Plane stress
For the plane stress case, the material stiffness matrix is found by imposing
\partial \sigma_ {3 3} = 0
on the general material stiffness matrix obtained for the plane strain case.
Uniaxial stress
For the uniaxial stress case the material stiffness matrix is available directly from the variation of Equation 4.3.2-13 as
\partial \sigma_ {1 1} = [ \frac {3}{2} Q - R \sigma_ {1 1} ^ {2} ] \partial \varepsilon_ {1 1}.
4.3.3 Stress potentials for anisotropic metal plasticity
The metal plasticity models in ABAQUS use the Mises stress potential for isotropic behavior and the Hill stress potentials for anisotropic behavior. Both of these potentials depend only on the deviatoric stress, so the plastic part of the response is incompressible. This means that, in cases where the plastic flow dominates the response (such as limit load calculations or metal forming problems), except for plane stress problems, the finite elements should be chosen so that they can accommodate the incompressible flow. Usually the reduced integration elements are used for this purpose: in ABAQUS/Standard the "hybrid" elements can also be used, at higher cost. The fully integrated first-order continuum elements in ABAQUS/Standard use selectively reduced integration, whereby the volumetric strain is calculated at the centroid of the element only. Those elements that are described in ``Solid isoparametric quadrilaterals and hexahedra, '' Section 3.2.4, are also suitable for such problems.
The Mises stress potential is
f (\pmb {\sigma}) = q,
where
q = \sqrt {\frac {3}{2} \mathbf {S} : \mathbf {S}},
in which S is the deviatoric stress:
\mathbf {S} = \pmb {\sigma} - \frac {1}{3} \mathrm{trace} (\pmb {\sigma}) \mathbf {I} = \pmb {\sigma} - \frac {1}{3} \mathbf {I} \mathbf {I}: (\pmb {\sigma}).
The potential is a circle in the plane normal to the hydrostatic axis in principal stress space. For this function,
\frac {\partial f}{\partial \pmb {\sigma}} = \frac {1}{q} \frac {3}{2} \mathbf {S},
and
\frac {\partial^ {2} f}{\partial \pmb {\sigma} \partial \pmb {\sigma}} = \frac {1}{q} \left(\frac {3}{2} \mathfrak {I} - \frac {1}{2} \mathbf {I} \mathbf {I} - \frac {\partial f}{\partial \pmb {\sigma}} \frac {\partial f}{\partial \pmb {\sigma}}\right)
in which is the fourth-order unit tensor.
Hill's stress function is a simple extension of the Mises function to allow anisotropic behavior. The function is
f (\pmb {\sigma}) = \sqrt {F (\sigma_ {y} - \sigma_ {z}) ^ {2} + G (\sigma_ {z} - \sigma_ {x}) ^ {2} + H (\sigma_ {x} - \sigma_ {y}) ^ {2} + 2 L \tau_ {y z} ^ {2} + 2 M \tau_ {z x} ^ {2} + 2 N \tau_ {x y} ^ {2}},
in terms of rectangular Cartesian stress components, where F; G; H; L; M; N are constants obtained by tests of the material in different orientations. They are defined as
F = \frac {\sigma_ {0} ^ {2}}{2} \left(\frac {1}{\bar {\sigma} _ {2 2} ^ {2}} + \frac {1}{\bar {\sigma} _ {3 3} ^ {2}} - \frac {1}{\bar {\sigma} _ {1 1} ^ {2}}\right),
G = \frac {\sigma_ {0} ^ {2}}{2} \left(\frac {1}{\bar {\sigma} _ {3 3} ^ {2}} + \frac {1}{\bar {\sigma} _ {1 1} ^ {2}} - \frac {1}{\bar {\sigma} _ {2 2} ^ {2}}\right),
H = \frac {\sigma_ {0} ^ {2}}{2} \left(\frac {1}{\bar {\sigma} _ {1 1} ^ {2}} + \frac {1}{\bar {\sigma} _ {2 2} ^ {2}} - \frac {1}{\bar {\sigma} _ {3 3} ^ {2}}\right),
L = \frac {3}{2} \left(\frac {\tau_ {0}}{\bar {\tau} _ {2 3}}\right) ^ {2},
M = \frac {3}{2} \big (\frac {\tau_ {0}}{\bar {\tau} _ {1 3}} \big) ^ {2},
N = \frac {3}{2} \left(\frac {\tau_ {0}}{\bar {\tau} _ {1 2}}\right) ^ {2},
where \sigma _ { 0 } , \bar { \sigma } _ { 1 1 } , \bar { \sigma } _ { 2 2 } , \bar { \sigma } _ { 3 3 } , \bar { \tau } _ { 1 2 } , \bar { \tau } _ { 2 3 } , \bar { \tau } _ { 1 3 } are specified by the user and \tau _ { 0 } = \sigma _ { 0 } / \sqrt { 3 } . . ¾¹ and \bar { \tau } are the values of stress that make the potential equal to \sigma _ { 0 } if only one stress component is nonzero.
For this function
\frac {\partial f}{\partial \pmb {\sigma}} = \frac {1}{f} \mathbf {b},
where
\mathbf {b} = \left[ \begin{array}{c} - G (\sigma_ {z} - \sigma_ {x}) + H (\sigma_ {x} - \sigma_ {y}) \\ F (\sigma_ {y} - \sigma_ {z}) - H (\sigma_ {x} - \sigma_ {y}) \\ - F (\sigma_ {y} - \sigma_ {z}) + G (\sigma_ {z} - \sigma_ {x}) \\ 2 N \tau_ {x y} \\ 2 M \tau_ {z x} \\ 2 L \tau_ {y z} \end{array} \right].
In addition,
\frac {\partial^ {2} f}{\partial \pmb {\sigma} \partial \pmb {\sigma}} = \frac {1}{f} \left(\frac {\partial \mathbf {b}}{\partial \pmb {\sigma}} - \frac {1}{f ^ {2}} \mathbf {b} \mathbf {b}\right),
where
\frac {\partial \mathbf {b}}{\partial \pmb {\sigma}} = \left[ \begin{array}{c c c c c c} G + H & - H & - G & 0 & 0 & 0 \\ - H & F + H & - F & 0 & 0 & 0 \\ - G & - F & F + G & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 N & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 M & 0 \\ 0 & 0 & 0 & 0 & 0 & 2 L \end{array} \right].
4.3.4 Rate-dependent metal plasticity (creep)
The rate-dependent plasticity (creep) models provided in ABAQUS/Standard are used to model inelastic straining of materials that are rate sensitive. High-temperature "creep" in structures is one important class of examples of the application of such a material model. Because such problems generally involve relatively small amounts of inelastic straining (otherwise the structure is not a suitable design), the explicit, forward Euler, method is often satisfactory as an integrator for the flow rule. This method is only conditionally stable, but the stability limit is usually sufficiently large