김경종
32 KiB
define the reduced invariants
_ {0} ^ {t} J _ {1} = _ {0} ^ {t} I _ {1} (_ {0} ^ {t} I _ {3}) ^ {- 1 / 3}
{ } _ { 0 } ^ { \prime } J _ { 2 } = { } _ { 0 } ^ { \prime } I _ { 2 } ( { } _ { 0 } ^ { \prime } I _ { 3 } ) ^ { - 2 / 3 } \tag {6.202}
_ {0} ^ {I} J _ {3} = \left(_ {0} ^ {I} I _ {3}\right) ^ {1 / 2}
and we then use
_ {0} \overline {{W}} = C _ {1} (_ {0} ^ {t} J _ {1} - 3) + C _ {2} (_ {0} ^ {t} J _ {2} - 3) + \frac {1}{2} \kappa (_ {0} ^ {t} J _ {3} - 1) ^ {2} \tag {6.203}
where \kappa is the bulk modulus. We note that in this description
^ {t} \bar {p} = - \kappa \left(_ {0} ^ {t} J _ {3} - 1\right) \tag {6.204}
The relations in (6.203) and (6.204) are used to calculate, by chain differentiation, all required derivatives in the displacement/pressure formulation in Section 6.4.
The basic (three-term) Ogden material description uses the form
{ } _ { 0 } ^ { t } \tilde { W } = \sum _ { n = 1 } ^ { 3 } \frac { \mu _ { n } } { \alpha _ { n } } ( \lambda _ { 1 } ^ { \alpha _ { n } } + \lambda _ { 2 } ^ { \alpha _ { n } } + \lambda _ { 3 } ^ { \alpha _ { n } } - 3 ) ; \quad \lambda _ { 1 } \lambda _ { 2 } \lambda _ { 3 } = 1 \tag {6.205}
where the \lambda_{i} are the principal values of the stretch tensor \mathbf{\delta}_{0}\mathbf{U} and the \mu_{i} and \alpha_{i} are material constants. We note that \frac{1}{2}\sum_{n=1}^{3}\alpha_{n}\mu_{n} is the small strain shear modulus.
The material description is more effectively used with the principal values L_{i} of _{0}^{i}\mathbf{C} [instead of the principal values of _{0}^{i}\mathbf{U} , where _{0}^{i}\mathbf{C} = (_{0}^{i}\mathbf{U})^{2} ]. Then we have
{ } _ { 0 } ^ { t } \tilde { W } = \sum _ { n = 1 } ^ { 3 } \frac { \mu _ { n } } { \alpha _ { n } } \left( L _ { 1 } ^ { \alpha _ { n } / 2 } + L _ { 2 } ^ { \alpha _ { n } / 2 } + L _ { 3 } ^ { \alpha _ { n } / 2 } - 3 \right) ; \quad L _ { 1 } L _ { 2 } L _ { 3 } = 1 \tag {6.206}
As for the Mooney-Rivlin material model, we now assume that the rubberlike material is only almost incompressible and also replace the L_{i} with the terms L_{i}(L_{1}L_{2}L_{3})^{-1/3} so as to render the three terms under the summation sign unaffected by volumetric deformations. The modified Ogden material description is then given by
_ {0} ^ {t} \overline {{W}} = \sum_ {n = 1} ^ {3} \left\{\frac {\mu_ {n}}{\alpha_ {n}} \left[ \left(L _ {1} ^ {\alpha_ {n} / 2} + L _ {2} ^ {\alpha_ {n} / 2} + L _ {3} ^ {g _ {n} / 2}\right) \left(L _ {1} L _ {2} L _ {3}\right) ^ {- \alpha_ {n} / 6} - 3 \right] \right\} + \frac {1}{2} \kappa \left(_ {0} ^ {t} J _ {3} - 1\right) ^ {2} \tag {6.207}
Using (6.207) and (6.204), we can now directly obtain, by differentiation, all terms of the displacement/pressure formulation given in Section 6.4 (see T. Sussman and K. J. Bathe [B]).
In the foregoing presentation we implicitly assumed that the total Lagrangian formulation is used for the analysis. Since the strain energy densities with the material constants are defined for this formulation, it is most natural and effective to use a total Lagrangian solution. However, an updated Lagrangian solution, given the values of \overline{0W} above, could also be pursued and identical numerical results would be obtained if the transformation rules in Section 6.4.2 (or 6.6.1) are followed.
Finally, we should also note that the basic Mooney-Rivlin and Ogden material models can be generalized in a straightforward manner to higher-order models with more terms in the expressions for \tilde{W} (see Exercise 6.69), see also T. Sussman and K.J. Bathe [C] for another approach.
6.6.3 Inelastic Material Behavior; Elastoplasticity, Creep, and Viscoplasticity
A fundamental observation comparing elastic and inelastic analysis is that in elastic solutions the total stress can be evaluated from the total strain alone [as given in (6.184) and (6.191)], whereas in an inelastic response calculation the total stress at time t also depends on the stress and strain history. Typical inelastic phenomena are elastoplasticity, creep, and viscoplasticity, and a very large number of material models have been developed in order to characterize such a material response. Our objective is again not to summarize or survey the models available but rather to present some of the basic finite element procedures that are employed in inelastic response calculations. We are mainly concerned with the general approach followed in inelastic finite element analysis and some formulations and numerical procedures that are used efficiently.
In the incremental analysis of inelastic response, basically three kinematic conditions are encountered.
Small displacement and small strain conditions:
In this case a materially-nonlinear-only formulation is used that assumes infinitesimally small displacements and rotations and considers only material nonlinearities. As long as the material is elastic, the solution using this formulation is identical to the linear elastic solution discussed in Section 4.2.
Large displacements and rotations but small strains:
In this case the total Lagrangian formulation is employed effectively. As discussed in Section 6.2.3, the kinematic assumptions permit large displacements, large rotations, and large strains. However, assuming small strains, the material model used for materially-nonlinear-only analysis can be directly employed in the TL formulation for large displacement and large rotation analysis by simply substituting in the material characterization the second Piola-Kirchhoff stresses and Green-Lagrange strains for the small displacement engineering stress and strain measures (see Section 6.6.1, Fig. 6.7) — with the proviso that with small material strain hardening, a large strain analysis can be computationally more stable and be physically more accurate.
Large displacements and large strains:
In this case a total or updated Lagrangian formulation can be employed efficiently. However, the underlying constitutive formulations are more complex—although they are direct extensions of the materially-nonlinear-only and large displacement—small strain cases. There are various issues that need to be discussed specifically for the large strain analysis case, and for this reason we defer the presentation to Section 6.6.4.
Since the large displacement–small strain case represents, from the computational view, a simple extension of the materially-nonlinear-only case, we consider in this section only inelastic conditions in small displacements and small strains. However, Table 6.9 shows how in elastoplastic analysis some of the major equations, further discussed below, are simply used with second Piola-Kirchhoff stresses and Green-Lagrange strains to include large displacement effects.
TABLE 6.9 Continuum elastoplasticity formulations
| Materially-nonlinear-only formulation (infinitesimally small displacements): | |||
| $\int_{V} C_{ijrs}^{EP} e_{rs} \delta e_{ij} \, dV = {}^{t+\Delta t} \Re - \int_{V} {}^{t} \sigma_{ij} \, \delta e_{ij} \, dV$ | |||
| ${}^{t}F(^{t} \sigma_{ij}, {}^{t} \kappa, \ldots) = 0;$ | $de_{ij}^{P} = d\lambda \frac{\partial^{t} F}{\partial^{t} \sigma_{ij}};$ | $d\sigma_{ij} = C_{ijrs}^{E}(de_{rs} - de_{rs}^{P});$ | ${}^{t+dt} \sigma_{ij} = {}^{t} \sigma_{ij} + d\sigma_{ij}$ |
Total Lagrangian formulation (large displacements and large rotations, but small strains),
\delta \bar {\epsilon} ^ {P} = \int_ {0} ^ {1} [ 2 / 3 d _ {0} \epsilon_ {i j} ^ {P} d _ {0} \epsilon_ {i j} ^ {P} ] ^ {1 / 2} < 2 \%:
\int C _ {i j r s} ^ {E P} _ {0} e _ {r s} \delta_ {0} e _ {i j} d ^ {0} V + \int \delta S _ {i j} \delta_ {0} \eta_ {i j} d ^ {0} V = ^ {r + \Delta t} \mathcal {R} - \int \delta S _ {i j} \delta_ {0} e _ {i j} d ^ {0} V
{ } ^ { t } F ( { } _ { 0 } ^ { t } S _ { i j } , { } ^ { t } \kappa , \dots ) = 0 ; \quad d _ { 0 } \epsilon _ { i j } ^ { P } = d \lambda \frac { \partial ^ { t } F } { \partial _ { 0 } ^ { t } S _ { i j } } ; \quad d _ { 0 } S _ { i j } = { } _ { 0 } C _ { i j r s } ^ { E } ( d _ { 0 } \epsilon _ { r s } - d _ { 0 } \epsilon _ { r s } ^ { P } ) ; \quad { } ^ { t + d } { } _ { 0 } ^ { t } S _ { i j } = { } _ { 0 } ^ { t } S _ { i j } + d _ { 0 } S _ { i j }
Referring to the summary of computations in Table 6.8, let us assume that the solution has been obtained accurately to time t and that the total strains t+\Delta t e_{rs} corresponding to time t + \Delta t have been computed (we now omit any iteration superscript). Hence, we assume that all stresses, inelastic strains, and state variables for time t are known accurately. We then have two basic requirements for our solution scheme:
The computation of the stresses, inelastic strains, and state variables corresponding to the total strains at time t + \Delta t .
The computation of the tangent constitutive relation corresponding to the state evaluated above, which in materially-nonlinear-only analysis can be written as
C _ {i j r s} ^ {I N} = \frac {\partial^ {t + \Delta t} \sigma_ {i j}}{\partial^ {t + \Delta t} e _ {r s}} \tag {6.208}
The tangent stress-strain law is employed in the evaluation of the tangent stiffness matrix. If the stress-strain relationship used in the evaluation of the stiffness matrix is not a true tangent relation, then the next calculated displacement (and hence strain) increment will, in general, not give as accurate an approximation to the solution sought as possible (and this will decrease the rate of convergence of the equilibrium iterations; see Section 8.4).
A crucial requirement is that the stresses for the new state be accurately calculated. If we denote the calculated total strains as t^{+\Delta t}\mathbf{e} , then our requirement is to obtain the stresses t^{+\Delta t}\sigma accurately. We should note that, in general, any error introduced in the evaluation of t^{+\Delta t}\sigma is an error that cannot be compensated for later in the solution by some corrective iterative scheme. Instead, errors present in the stresses and plastic strains, and the state variables for time t + \Delta t will, in general, irreversibly deteriorate the subsequent response prediction.
Let us point out here that of course in an equilibrium iteration, the tangent constitutive relation in (6.10) needs to be calculated corresponding to the current state, and also, the stresses to be calculated are t+\Delta t \sigma^{(i-1)} for the given strains t+\Delta t e^{(i-1)} (see Section 8.4).
In the following presentation, we first consider the basic computational scheme for the case of elastoplasticity, and then we briefly consider the cases of creep and viscoplasticity.
Elastoplasticity
We assume that the elastoplastic response is governed by the classical incremental theory of plasticity based on the Prandtl-Reuss equations (see, for example, L. E. Malvern [A], R. Hill [B], A. Mendelson [A], and M. Życzkowski [A]).
Using the additive decomposition of strains, de_{ij} = de_{ij}^{E} + de_{ij}^{P} , the stress increment is given by the basic relationship
d \sigma_ {i j} = C _ {i j r s} ^ {E} \left(d e _ {r s} - d e _ {r s} ^ {P}\right) \tag {6.209}
where the C_{ijrs}^{E} are the components of the elastic constitutive tensor and the de_{ij} , de_{ij}^{E} , de_{ij}^{P} are the components of the total strain increment, the elastic strain increment, and the plastic strain increment, respectively. This incremental relationship holds throughout the inelastic response. To calculate the plastic strains, we use three properties to characterize the material behavior:
A yield function, which gives the yield condition that specifies the state of multiaxial stress corresponding to start of plastic flow
A flow rule, which relates the plastic strain increments to the current stresses and the stress increments
A hardening rule, which specifies how the yield function is modified during plastic flow.
The yield function has the general form at time t,
^ {t} f _ {y} \left(^ {t} \sigma_ {i j}, ^ {t} e _ {i j} ^ {P}, \dots\right) \tag {6.210}
where “ . . . ” denotes state variables that depend on the material characterization. The instantaneous material response is elastic if
^ \prime f _ {y} < 0 \tag {6.211}
and elastic or plastic depending on the loading condition if
f _ {y} = 0 \tag {6.212}
whereas f_y > 0 is inadmissible. Hence the relation (6.212) represents the yield condition, which must hold throughout the plastic response.
Assuming that for the material the associated flow rule is applicable during plastic response we use the function f_y in the flow rule to obtain the plastic strain increments
d e _ {i j} ^ {P} = d \lambda \frac {\partial {} ^ {t} f _ {y}}{\partial {} ^ {t} \sigma_ {i j}} \tag {6.213}
where d\lambda is a scalar to be determined [and depending on the state variables additional quantities appear in (6.213)]. The hardening rule—which also depends on the particular
material model used—changes the state variables in f_{y} as a consequence of plastic flow and therefore changes the yield condition during the response.
Let us consider von Mises plasticity with isotropic hardening, and a general three-dimensional stress state. In the following we present a simple solution procedure that is widely used and referred to as the radial return method (see M. L. Wilkins [A] and R. D. Krieg and D. B. Krieg [A]). Our objective is to present the procedure such that it can be seen to be the basis of a general approach toward the solution of inelastic stress conditions.
Since in von Mises plasticity the plastic volumetric strains are zero (as we shall see later), it is effective to write the general stress-strain relationship at time t + \Delta t in the form [see also (4.125) to (4.133)]
{ } ^ { t + \Delta t } \mathbf { S } = \frac { E } { 1 + \nu } ( { } ^ { t + \Delta t } \mathbf { e } ^ { t } - { } ^ { t + \Delta t } \mathbf { e } ^ { P } ) \tag {6.214}
{ } ^ { t + \Delta t } \sigma _ { m } = \frac { E } { 1 - 2 \nu } { } ^ { t + \Delta t } e _ { m } \tag {6.215}
where t+\Delta t S is the deviatoric stress tensor with components ^{13}
{ } ^ { t + \Delta t } S _ { i j } = { } ^ { t + \Delta t } \sigma _ { i j } - { } ^ { t + \Delta t } \sigma _ { m } \delta _ { i j } \tag {6.216}
^{t+\Delta t}\sigma_{m} is the mean stress
{ } ^ { t + \Delta t } \sigma _ { m } = \frac { { } ^ { t + \Delta t } \sigma _ { i i } } { 3 } \tag {6.217}
^{t+\Delta t}e' is the deviatoric strain tensor with components
{ } ^ { t + \Delta t } e _ { i j } ^ { \prime } = { } ^ { t + \Delta t } e _ { i j } - { } ^ { t + \Delta t } e _ { m } \delta _ { i j } \tag {6.218}
t + \Delta t e_{m} is the mean strain
{ } ^ { t + \Delta t } e _ { m } = \frac { { } ^ { t + \Delta t } e _ { i i } } { 3 } \tag {6.219}
and ^{t+\Delta t}\mathbf{e}^{P} is the plastic strain tensor with components ^{t+\Delta t}e_{ij}^{P} . Note that in (6.214) we have 3 \times 3 matrices on each side of the equal sign, with the components of the tensors. We should recall that the total strain at time t + \Delta t is known (see Table 6.8), and hence, ^{t+\Delta t}\mathbf{e}^{t} and ^{t+\Delta t}e_{m} are known.
The relations (6.214) and (6.215) represent the integrated forms of (6.209), where we note that the deviatoric stresses depend on the plastic strains, which are in general highly dependent on the stress history, and that the mean stress is independent of the plastic strains (because the plastic mean strain is zero; i.e., the plastic deformation is isochoric). Hence, t^{+\Delta t}\sigma_{m} is given directly by (6.215), and our task is now to calculate t^{+\Delta t}e^{P} and t^{+\Delta t}S .
Since we assume that the complete stress and strain conditions are known at time t , we can write (6.214) in the form
{ } ^ { t + \Delta t } \mathbf { S } = \frac { E } { 1 + \nu } ( { } ^ { t + \Delta t } \mathbf { e } ^ { \prime \prime } - \Delta \mathbf { e } ^ { P } ) \tag {6.220}
where t^{+\Delta t}e^{\prime\prime}=t^{+\Delta t}e^{\prime}-t e^{P} (6.221)
is a known quantity.
Hence, the task of integrating the constitutive relations is now reduced to determining in (6.220) the stresses {}^{t+\Delta t}S and the plastic incremental strains \Delta e^{P} subject to the yield condition, flow rule, and hardening rule.
In von Mises plasticity, the yield condition is at time t + \Delta t ,
{ } ^ { t + \Delta t } f _ { y } ^ { o M } = \frac { 1 } { 2 } { } ^ { t + \Delta t } \mathbf { S } \cdot { } ^ { t + \Delta t } \mathbf { S } - \frac { 1 } { 3 } ( { } ^ { t + \Delta t } \sigma _ { y } ) ^ { 2 } = 0 \tag {6.222}
where t^{+}\Delta t \sigma_{y} is the yield stress at time t + \Delta t . This stress is a function of the effective plastic strain, which defines the hardening of the material,
t + \Delta t\sigma_{y} = f_{e}(t + \Delta t\overline{e}^{P}) (6.223) with t + \Delta t\overline{e}^{P} = \int_{0}^{t + \Delta t}\sqrt{\frac{2}{3}de^{P}\cdot de^{P}} (6.224)
Figure 6.9 shows a yield curve schematically.
line
| Effective plastic strain | Effective stress and yield stress | | ------------------------ | ---------------------------------- | | t+Δtσ̄ = t+Δtσy | τf_y = 0 | | tσ̄ = tσy | τf_y < 0 | | Δe̅^P | τf_y < 0 | | t_ē^P | τf_y < 0 | | t+Δtē^P | τf_y < 0 |Figure 6.9 General generic yield curve. Yielding is assumed at times t and t + \Delta t . Let, at time \tau , \tau f_{y} = \tau \overline{\sigma} - \tau \sigma_{y} ; then the response must satisfy \tau f_{y} \leq 0 , \dot{\overline{e}}^{P} \geq 0 , \dot{\overline{e}}^{P} \tau f_{y} = 0 .
The flow rule gives for the finite step
\Delta \mathbf {e} ^ {P} = \lambda \frac {\partial^ {t + \Delta t} f _ {y} ^ {o M}}{\partial^ {t + \Delta t} \mathbf {G}} = \lambda^ {t + \Delta t} \mathbf {S} \tag {6.225}
Geometrically, this equation means that \Delta e^{P} is in the direction of ^{t+\Delta t}S , and the equation shows that the volumetric plastic strains are zero (the trace of ^{t+\Delta t}S is zero).
To determine \lambda we take the dot products of the tensor components on each side of (6.225) to obtain ^{15}
\Delta \overline {{e}} ^ {P} = \frac {2}{3} \lambda^ {t + \Delta t} \overline {{\sigma}} \tag {6.226}
where \Delta\overline{e}^{P} is the increment in the effective plastic strain and t+\Delta t\overline{\sigma} is the effective stress
\Delta \overline {{e}} ^ {P} = \sqrt {\frac {2}{3} \Delta \mathbf {e} ^ {P} \cdot \Delta \mathbf {e} ^ {P}}; \quad t + \Delta t \overline {{\sigma}} = \sqrt {\frac {3}{2} t + \Delta t \mathbf {S} \cdot t + \Delta t \mathbf {S}} \tag {6.227}
^{14} Note that here and in the following we use the notation that for second-order tensors a and b we have a \cdot b = a_{ij}b_{ij} [sum over all i, j ; see (2.79)].
^{15} Note that geometrically \lambda simply gives the difference in the “lengths” of \Delta e^{P} and ^{t+\Delta t}S .
The relations in (6.227) correspond to the plastic strain and the stress in a uniaxial stress condition (see Exercise 6.73), and using (6.225) we have that in general three-dimensional analysis the plastic work \sigma_{ij}de_{ij}^{P} is given by \overline{\sigma}d\overline{e}^{P} . From (6.226)
\lambda = \frac {3}{2} \frac {\Delta \bar {e} ^ {P}}{r + \Delta r \bar {\sigma}} \tag {6.228}
However, using the requirement that during yielding the yield stress is equal to the effective stress (because the yield condition is satisfied), we can relate \Delta \overline{e}^P to the effective stress at time t + \Delta t (and other known variables). This relation is obtained from the effective stress-effective plastic strain curve shown schematically in Fig. 6.9.
We next substitute from (6.225) into (6.220) to obtain
{ } ^ { t + \Delta t } \mathbf { S } = \frac { 1 } { a _ { E } + \lambda } { } ^ { t + \Delta t } \mathbf { e } ^ { n } \tag {6.229}
and so t^{+\Delta t} \mathbf{e}'' is also in the direction of t^{+\Delta t} \mathbf{S} . To evaluate the scalar multiplier in (6.229) we can again simply take the dot product of the tensor components on each side of (6.229). We obtain
a ^ {2 t + \Delta t} \overline {{{\sigma}}} ^ {2} - d ^ {2} = 0 \tag {6.230}
where a = a_{E} + \lambda ;\quad a_{E} = \frac{1 + \nu}{E} (6.231)
d ^ {2} = \frac {3}{2} ^ {t + \Delta t} \mathrm{e} ^ {\prime \prime} \cdot {} ^ {t + \Delta t} \mathrm{e} ^ {\prime \prime} \tag {6.232}
We note that the coefficient d is constant, whereas the coefficient a varies with \lambda , and hence with ^{t + \Delta t}\overline{\sigma} .
Let us define the function
f (\overline {{{\sigma}}} ^ {*}) = a ^ {2} (\overline {{{\sigma}}} ^ {*}) ^ {2} - d ^ {2} \tag {6.233}
and call f(\overline{\sigma}^{*}) the effective stress function. Considering (6.230) we recognize that f = 0 at \overline{\sigma}^{*} = ^{t + \Delta t}\overline{\sigma} . Hence, solving for the zero value of the effective stress function provides the solution for ^{t + \Delta t}\overline{\sigma} and \lambda , and hence, the solution for the current stress state ^{t + \Delta t}\mathbf{S} [see (6.229)] and the incremental plastic strains \Delta \mathbf{e}^{P} [see (6.225)]. Since the key step in the solution—here and in more complex inelastic analysis—is the calculation of the zero of a function f(\overline{\sigma}^{*}) , the complete solution algorithm has been termed the effective-stress-function (ESF) algorithm (see K. J. Bathe, A. B. Chaudhary, E. N. Dvorkin, and M. Kojić [A], where the method is introduced for more complex thermoelastoplastic and creep solutions).
The attractiveness of the ESF algorithm lies in the general applicability of the method. The effective stress function can be quite complicated (because of a complicated effective stress–effective plastic strain relationship) or because of thermal and creep effects (as considered later). The zero of the function is in general efficiently found by a numerical bisection scheme, which can be made very stable because the function contains only one unknown, the effective stress.
However, if the material is modeled as bilinear (see Fig. 6.10), the effective stress–effective strain relationship is a straight line with slope E_{P} and we can solve directly for the effective stress at time t + \Delta t . Namely,
\Delta \overline {{e}} ^ {P} = \frac {^ {t + \Delta t} \overline {{\sigma}} - ^ {t} \sigma_ {y}}{E _ {P}} \tag {6.234}
line
| Total strain | Stress |
|---|---|
| e^P | E |
| e^E | E_T |
text_image
Stress σyv EP Plastic strain eP
Figure 6.10 Bilinear elastoplastic material model
with
E _ {P} = \frac {E E _ {T}}{E - E _ {T}} \tag {6.235}
where E_{T} is the tangent modulus, and hence,
^ {t + \Delta t} \overline {{{\sigma}}} = \frac {2 E _ {P} d + 3 ^ {t} \sigma_ {y}}{2 E _ {P} a _ {E} + 3} \tag {6.236}
In case of perfect plasticity E_P = 0 and r^{+ \Delta t} \overline{\sigma} = \sigma_{yv} , so that from (6.230),
\lambda = \frac {d}{\sigma_ {y v}} - a _ {E} \tag {6.237}
The elastoplasticity computations are of course started by first testing whether the material is yielding from time t to time t + \Delta t . Yielding is assumed to be present throughout the time interval if t^{+\Delta t}\overline{\sigma}^{E} > {}^{t}\sigma_{y} , where {}^{t+\Delta t}\overline{\sigma}^{E} is the elastic stress solution
{ } ^ { t + \Delta t } \overline { { { \sigma } } } ^ { E } = \frac { d } { a _ { E } } \tag {6.238}
The above solution process can be interpreted to consist of, first, an elastic prediction of stress (in which \Delta e^{p} is assumed to be zero) and then, if this stress prediction lies outside the yield surface corresponding to time t, a stress correction. Figure 6.11 illustrates geometrically the solution process.
text_image
S33 tfvM = 0 t+ΔtfvM = 0 tS A C B t+Δtn S11 t+Δts S22
Figure 6.11 Geometric representation of stress solution in elastoplastic analysis (using principal stress directions)
Point A corresponds to the stress solution at time t . Point B corresponds to the elastic stress prediction, which from (6.220) is
{ } ^ { t + \Delta t } \mathbf { S } ^ { E } = \frac { E } { 1 + \nu } ( { } ^ { t + \Delta t } \mathbf { e } ^ { \prime \prime } ) \tag {6.239}
The stress correction corresponds to the vector BC and is
^ {t + \Delta t} \mathbf {S} ^ {c} = \frac {E}{1 + \nu} (- \Delta \mathbf {e} ^ {P}) \tag {6.240}
where \Delta \mathbf{e}^P is given by (6.225).
Since ^{t+\Delta t}\mathbf{S} and ^{t+\Delta t}\mathbf{S}^{c} are in the same direction and ^{t+\Delta t}\mathbf{S} = ^{t+\Delta t}\mathbf{S}^{E} + ^{t+\Delta t}\mathbf{S}^{c} , we note that the normal ^{t+\Delta t}\mathbf{n} is in the direction of ^{t+\Delta t}\mathbf{S}^{E} (and ^{t+\Delta t}\mathbf{S} , \Delta\mathbf{e}^{P} , and ^{t+\Delta t}\mathbf{e}'' ) and is given by
{ } ^ { t + \Delta t } \mathbf { n } = \frac { { } ^ { t + \Delta t } \mathbf { e } ^ { \prime \prime } } { \| { } ^ { t + \Delta t } \mathbf { e } ^ { \prime \prime } \| _ { 2 } } . \tag {6.241}
Since the stress correction along this vector returns the stress state to the yield surface radially, the method has been called the radial return method (see R. D. Krieg and D. B. Krieg [A]).
This interpretation of the solution scheme gives some indication of the numerical error in the solution process. We see that the yield condition and hardening law are accurately satisfied at the discrete solution times, but the magnitude of the plastic strains is obtained in (6.225) by an estimate based on the Euler backward integration for the flow rule. If the direction of the deviatoric stress vector does not change (hence \tau^{n} is constant for all \tau ), the solution is exact, but any change in the direction of this vector will introduce numerical integration errors. Some studies on solution errors are presented by R. D. Krieg and D. B. Krieg [A], H. L. Schreyer, R. F. Kulak, and J. M. Kramer [A], M. Ortiz and E. P. Popov [A], and N. S. Lee and K. J. Bathe [B]. These investigations indicate the high accuracy that is achieved in using the solution algorithm.
As we have pointed out, the second important requirement of the computational scheme is the evaluation of an accurate tangent stress-strain relation for the stiffness matrix. This tangent stress-strain matrix must be consistent with the numerical stress integration scheme used in order to have the full convergence benefits of the Newton-Raphson iteration (see Section 8.4.1). Let us consider that we require the tangent material matrix corresponding to time t + \Delta t . We assume the schematic yield curve in Fig. 6.9 is applicable.
For ease of writing let us define the stress and strain vectors
\begin{array}{l} { } ^ { t + \Delta t } \boldsymbol { \sigma } ^ { T } = \left[ { } ^ { t + \Delta t } \sigma _ { 1 1 } \quad { } ^ { t + \Delta t } \sigma _ { 2 2 } \quad { } ^ { t + \Delta t } \sigma _ { 3 3 } \quad { } ^ { t + \Delta t } \sigma _ { 1 2 } \quad { } ^ { t + \Delta t } \sigma _ { 2 3 } \quad { } ^ { t + \Delta t } \sigma _ { 3 1 } \right] \\ = \left[ ^ {t + \Delta t} \sigma_ {1} \quad {} ^ {t + \Delta t} \sigma_ {2} \quad {} ^ {t + \Delta t} \sigma_ {3} \quad {} ^ {t + \Delta t} \sigma_ {4} \quad {} ^ {t + \Delta t} \sigma_ {5} \quad {} ^ {t + \Delta t} \sigma_ {6} \right] \tag {6.242} \\ \end{array}
\begin{array}{l} { } ^ { t + \Delta t } \mathbf { e } ^ { T } = \left[ { } ^ { t + \Delta t } e _ { 1 1 } \quad { } ^ { t + \Delta t } e _ { 2 2 } \quad { } ^ { t + \Delta t } e _ { 3 3 } \quad { } ^ { t + \Delta t } \gamma _ { 1 2 } \quad { } ^ { t + \Delta t } \gamma _ { 2 3 } \quad { } ^ { t + \Delta t } \gamma _ { 3 1 } \right] \\ = \left[ ^ {t + \Delta t} e _ {1} ^ {t + \Delta t} e _ {2} ^ {t + \Delta t} e _ {3} ^ {t + \Delta t} e _ {4} ^ {t + \Delta t} e _ {5} ^ {t + \Delta t} e _ {6} ^ {t + \Delta t} \right] \tag {6.243} \\ \end{array}
where we use the engineering shear strains t^{+\Delta t}\gamma_{ij}=t^{+\Delta t}e_{ij}+t^{+\Delta t}e_{ji} . Then the tangent stress-strain matrix at time t+\Delta t is given by the derivative of the stresses with respect to the strains, evaluated at time t+\Delta t , which we write as
\mathbf {C} ^ {E P} = \frac {\partial^ {t + \Delta t} \mathbf {0}}{\partial^ {t + \Delta t} \mathbf {e}} \tag {6.244}
We now need to establish the derivative in (6.244) using the assumptions of the above stress integration, and in particular the fact that the effective stress alone (or effective plastic strain alone) determines the current stress state. Hence, using the stress integration assumption and that t^{+\Delta t}\sigma is a function of t^{+\Delta t}e and t^{+\Delta t}\overline{\sigma}
\mathbf {C} ^ {E P} = \left. \frac {\partial^ {t + \Delta t} \boldsymbol {\sigma}}{\partial^ {t + \Delta t} \mathbf {e}} \right| _ {t + \Delta t _ {\overline {{\sigma}}} = \text { const. }} + \frac {\partial^ {t + \Delta t} \boldsymbol {\sigma}}{\partial^ {t + \Delta t} \overline {{\sigma}}} \frac {\partial^ {t + \Delta t} \overline {{\sigma}}}{\partial^ {t + \Delta t} \mathbf {e}} \tag {6.245}
or instead of the effective stress we can use the effective strain.
The derivation is achieved by careful differentiations. Let us present the first steps using for all stresses and strains the conventions in (6.242) and (6.243).
With the earlier established relations, we have
\left.\begin{array}{l}^{t + \Delta t}\sigma_{i} = ^{t + \Delta t}S_{i} + ^{t + \Delta t}\sigma_{m};\qquad i = 1,2,3\\^{t + \Delta t}\sigma_{i} = ^{t + \Delta t}S_{i};\qquad i = 4,5,6\end{array}\right\} \tag{6.246} and
Using (6.244), we obtain for the elements of the tangent stress-strain matrix,
\left. \begin{array}{c c} C _ {i j} ^ {E P} = C _ {i j} ^ {\prime} + \frac {E}{3 (1 - 2 \nu)}; & 1 \leq i, j \leq 3 \\ C _ {i j} ^ {E P} = C _ {i j} ^ {\prime}; & \text { otherwise } \end{array} \right\} \tag {6.247}
where C_{ij}^{\prime} = \frac{\partial^{t + \Delta t}S_{i}}{\partial^{t + \Delta t}e_{j}} (6.248)
However, ^{t + \Delta t}S_{i} is given by (6.229), in which \lambda is given by (6.228) and a_{E} is a constant. Therefore, (6.248) will ultimately involve a differentiation of \lambda (and hence, of \Delta \overline{e}^{P} and ^{t + \Delta t}\overline{\sigma} ) with respect to the strain components. Using the chain rule, we obtain
\frac {\partial^ {t + \Delta t} S _ {i}}{\partial^ {t + \Delta t} e _ {j}} = \frac {\partial^ {t + \Delta t} S _ {i}}{\partial^ {t + \Delta t} e _ {k} ^ {n}} \frac {\partial^ {t + \Delta t} e _ {k} ^ {n}}{\partial^ {t + \Delta t} e _ {l} ^ {t}} \frac {\partial^ {t + \Delta t} e _ {l} ^ {t}}{\partial^ {t + \Delta t} e _ {j}} \tag {6.249}
where from (6.221),
\frac {\partial^ {t + \Delta t} e _ {k} ^ {\prime \prime}}{\partial^ {t + \Delta t} e _ {l} ^ {\prime}} = \delta_ {k l} \tag {6.250}
and from (6.218),
\left. \begin{array}{l} \left[ \frac {\partial^ {t + \Delta t} e _ {i} ^ {\prime}}{\partial^ {t + \Delta t} e _ {j}} \right] = \frac {1}{3} \left[ \begin{array}{c c c} 2 & - 1 & - 1 \\ - 1 & 2 & - 1 \\ - 1 & - 1 & 2 \end{array} \right]; \quad 1 \leq l, j \leq 3 \\ \frac {\partial^ {t + \Delta t} e _ {i} ^ {\prime}}{\partial^ {t + \Delta t} e _ {j}} = \frac {1}{2} \delta_ {i j}; \quad \text { otherwise } \end{array} \right\} \tag {6.251}
Also, using (6.229), we obtain
\frac {\partial^ {t + \Delta t} S _ {i}}{\partial^ {t + \Delta t} e _ {k} ^ {n}} = \frac {1}{a _ {E} + \lambda} \delta_ {i k} - \frac {1}{(a _ {E} + \lambda) ^ {2}} \frac {\partial \lambda}{\partial^ {t + \Delta t} e _ {k} ^ {n}} ^ {t + \Delta t} e _ {i} ^ {n} \tag {6.252}