User:ZdenekPBazant

Objective Stress Rates in Finite Strain of Inelastic Solids

In many practical problems of solid mechanics, it is sufficient to characterize material deformation by the small (or linearized) strain tensor $e_{ij}={\frac {1}{2}}(u_{i,j}+u_{j,i})$ where $u_{i}$ are the components of the displacements of continuum points, the subscripts refer to Cartesian coordinates $x_{i}$ $(i=1,2,3)$, and the subscripts preceded by a comma denote partial derivatives (e.g., $u_{i,j}=\partial u_{i}/\partial x_{j}$). But there are also many problems where the finiteness of strain must be taken into account. These are of two kinds: 1) Large nonlinear elastic deformations possessing a potential energy, $\Pi$ (exhibited, e.g., by rubber), in which the stress tensor components are obtained as the partial derivatives of $\Pi$ with respect to the finite strain tensor components; and 2) inelastic deformations possessing no potential, in which the stress-strain relation is defined incrementally. In the former kind, the total strain formulation described in another article on Finite Strain is appropriate, and need not be described here.

In the latter kind, which is what is reviewed here, an incremental (or rate) formulation is necessary and must be used in every load or time step of a finite element computer program using updated Lagrangian procedure. The absence of a potential raises intricate questions due to the freedom in the choice of finite strain measure and characterization of the stress rate.

For a sufficiently small loading step (or increment), one may use the deformation rate tensor (or velocity strain) ${\dot {e}}_{ij}={\frac {1}{2}}(v_{i,j}+v_{j,i})$ or increment $\Delta e_{ij}={\dot {e}}_{ij}\Delta t$ representing the linearized strain increment from the initial (stressed and deformed) state in the step. Here the superior dot represents the time derivative $\partial /\partial t$, $\Delta$ denotes a small increment over the step, $t$ = time, and $v_{i}={\dot {u}}_{i}$ = material point velocity or displacement rate. However, it would not be objective to use the time derivative of the Cauchy (or true) stress $S_{ij}$. This stress, which describes the forces on a small material element imagined to be cut out from the material as currently deformed (Fig. 1), is not objective because it varies with rigid body rotations of the material. The material points must be characterized by their initial coordinates $x_{i}$ (called Lagrangian) because different material particles are contained in the element that is cut out (at the same location) before and after the incremental deformation.

Consequently, it is necessary to introduce the so-called objective stress rate ${\hat {S}}_{ij}$, or the corresponding increment $\Delta S_{ij}={\hat {S}}_{ij}\Delta t$. The objectivity is necessary for ${\hat {S}}_{ij}$ to be functionally related to the element deformation. It means that that ${\hat {S}}_{ij}$ must be invariant with respect to coordinate transformations (particularly rotations) and must characterize the state of the same material element as it deforms.

The objective stress rate can be derived in two ways: 1) By tensorial coordinate transformations, which is the standard way in finite element textbooks variationally, from strain energy density in the material expressed in terms of the strain tensor (which is objective by definition)). Although the former way may be instructive and provide geometric insight, it is more complicated and does not readily reveal whether the stress and strain tensors are work-conjugate, i.e., whether they give the correct expression for the second-order energy increment. In the variational energy approach, which is explained here, the objectivity and work-conjugacy are automatic.

Consider incremental finite strain tensors $\epsilon _{ij}$ relative to the initial (stressed) state at the beginning of the load step, using the Lagrangian frame of reference with material points of initial coordinates $x_{i}$ ($i=1,2,3$). Many different definitions of tensors $\epsilon _{ij}$ are admissible, provided that they all satisfy the conditions that: 1) $\epsilon _{ij}$ vanishes for all rigid-body motions; 2) the dependence of $\epsilon _{ij}$ on strain gradient tensor $u_{i,j}$ is continuous, continuously differentiable and monotonic. For convenience, it is also desired that $\epsilon _{ij}$ reduces to $e_{ij}$ as the norm $|u_{i,j}|\to 0$.

A broad class of equally admissible finite strain tensors, which comprises virtually all those ever used to solve mechanics problems, is represented by the Doyle-Ericksen tensors. In tensor notation, they are written as $\epsilon ^{(m)}=(U^{m}-I)/m$, where $m$ is a real parameter, $I$ = unit tensor and $U$ = right-stretch tensor. The second-order approximation of these tensors is

$\epsilon _{ij}^{(m)}=e_{ij}+{\frac {1}{2}}u_{k,i}u_{k,j}-{\frac {1}{2}}(2-m)e_{ki}e_{kj},$ (1)

The case $m=2$ gives the Green-Lagrangian strain tensor, $m=1$ gives the Biot strain tensor, $m=0$ gives the Hencky (logarithmic) strain tensor (for which $\epsilon ^{(0)}=\lim _{m\to 0}\epsilon ^{(m)}=\ln U$), and $m=-2$ gives the Almansi-Lagrangian strain tensor. Note also that the tensors $\epsilon _{ij}^{(n)}=\left(U^{n}-U^{-n}\right)/2n$, which do not belong to the foregoing class, have, for all $n$, the same 2nd-order approximation as Eq. 1 for $m=0$.

Relation of Objective Stress Rates to Finite Strain Tensors

The work $\delta W$ done at small deformations of a material element of unit initial volume, starting from an initial state under initial Cauchy (or true) stress $S_{ij}^{0}$, can be expressed (in Lagrangian coordinates of the initial state) in two ways:

$\delta W=(S_{ij}^{0}+\sigma _{ij}^{(m)})\delta \epsilon _{ij}^{(m)}~~~~{\mbox{or}}~~~~\delta W=(S_{ij}^{0}+\tau _{ij})\delta u_{i,j}$ (2)

where $\delta \epsilon _{ij}^{(m)}$ is an arbitrary strain variation, and $\tau _{ij}$ is a small Lagrangian (or first Piola-Kirchhoff) stress increment referred to the deformed configuration at the beginning of the stress increment; $\tau _{ij}$ is nonsymmetric and represents the forces acting at the faces of the incrementally deformed material element in the directions of initial Lagrangian coordinates $x_{i}$ (Fig. 1).

Furthermore, $\sigma _{ij}^{(m)}$ is a small stress increment during the load step which is symmetric and objective $\sigma _{ij}^{(m)}$ may be regarded as an increment of the second Piola-Kirchhoff stress referred to the stressed initial state). The first-order part of the work, $S_{ij}^{0}\delta u_{i,j}$, figures in the virtual work equation of static (or dynamic) equilibrium of the structure, in which it is cancelled out by the work of loads. The second-order part of the work, $\sigma _{ij}^{(m)}\delta \epsilon _{ij}^{(m)}$, decides stability, as well as bifurcation of the load path in the stress space. The objectivity of $\sigma _{ij}^{(m)}$ is ensured by independence of $\sigma _{ij}^{(m)}$ from rigid body motions (objectivity) and by the correctness of $\sigma _{ij}^{(m)}\delta \epsilon _{ij}^{(m)}$ as a second-order energy expression;

Now set the two expressions in Eq. (2) equal, and substitute $S_{ij}\delta u_{i,j}=S_{ij}\delta e_{ij}=S_{ij}{\dot {e}}_{ij}\Delta t$ (by virtue of symmetry of $S_{ij}$), $\sigma _{ij}^{(m)}\delta \epsilon _{ij}^{(m)}\approx \sigma _{ij}^{(m)}{\dot {e}}_{ij}\Delta t$ (which suffices for second-order work accuracy in $u_{i,j}$), $S_{ij}\delta \epsilon _{ij}^{(m)}=S_{pq}(\partial \epsilon _{pq}^{(m)}/\partial u_{i,j})v_{i,j}\Delta t$ and $\sigma _{ij}^{(m)}={\hat {S}}_{ij}^{(m)}\Delta t$ (where $v_{i,j}\Delta t=\delta u_{i,j}$, and $v_{i}={\dot {u}}_{i}$). Imposing the variational condition that the resulting equation must be valid for any strain gradient $\delta u_{i,j}$, one gets the following general expression for the objective stress rate associated with $m$:

${\hat {S}}_{ij}^{(m)}={\dot {T}}_{ij}-S_{pq}\ {\frac {\partial }{\partial t}}{\frac {\partial (\epsilon _{pq}^{(m)}-e_{pq})}{\partial u_{i,j}}}$ (3)

where

${\dot {T}}_{ij}={\dot {S}}_{ij}-S_{ik}v_{j,k}+S_{ij}v_{k,k}$ (4)

Here ${\dot {S}}_{ij}=\partial S_{ij}/\partial t$ = material rate of Cauchy stress (i.e., the rate in Lagrangian coordinates of the initial stressed state), and $T_{ij}=S_{ij}^{0}+\tau _{ij}$ (the Lagrangian stress, aka the first Piola-Kirchhoff stress), ${\dot {T}}_{ij}=\partial T_{ij}/\partial t=\partial \tau _{ij}/\partial t$.

Eq. (3) is the key equation from the energy viewpoint. There are infinitely many objective stress rates. But a rate for which there exists no legitimate finite strain tensor $\epsilon _{ij}$ associated according to Eq. (3) is energetically inconsistent, i.e., its use violates energy balance or the first law of thermodynamics.

Evaluating Eq. (3) for general $m$ and for $m=2$, one gets a general expression for the objective stress rate:

${\hat {S}}_{ij}^{(m)}={\hat {S}}_{ij}+{\frac {1}{2}}(2-m)(S_{ik}{\dot {e}}_{kj}+S_{jk}{\dot {e}}_{ki})$ (5)

where ${\hat {S}}_{ij}={\hat {S}}_{ij}^{(2)}$ = objective stress rate associated with the Green-Lagrangian strain ($m=2$), as the reference case. In particular,

${\hat {S}}_{ij}^{(2)}={\dot {S}}_{ij}-S_{kj}v_{i,k}-S_{ki}v_{j,k}+S_{ij}v_{k,k}{\mbox{ (Truesdell stress rate)}}$ (6)

${\hat {S}}_{ij}^{(0)}={\dot {S}}_{ij}-S_{kj}{\dot {\omega }}_{ik}+S_{ik}{\dot {\omega }}_{kj}+S_{ij}v_{k,k}{\mbox{ (Jaumann rate of Kirchhoff stress)}}$ (7)

where ${\dot {\omega }}_{jk}={\frac {1}{2}}(v_{i,j}-v_{j,i})$ = material rotation rate (spin tensor). Furthermore, $m=1$ gives the Biot stress rate, and $m=-2$ gives the Lie derivative of Kirchhoff stress.

Another rate, used in most commercial codes, is

${\hat {S}}_{ij}^{J}={\dot {S}}_{ij}-S_{kj}{\dot {\omega }}_{ik}+S_{ik}{\dot {\omega }}_{kj}{\mbox{ (Jaumann, or corotational, rate of Cauchy)}}$ (8)

It differs from Eq. (7) by missing the term $S_{ij}v_{k,k}$, in which $v_{k,k}={\dot {e}}_{kk}$ = rate of relative volume change of material. This causes that ${\hat {S}}_{ij}^{J}$ is not associated by work with any admissible finite strain measure. It follows that $\delta ^{2}W=({\hat {S}}_{ij}^{J}\Delta t)\delta \epsilon _{ij}^{(0)}$ is generally not the correct expression for the second-order work of stress increments, which is a violation of energy balance. ${\hat {S}}_{ij}^{J}$ happens to be the rate that is (by 2012) used in most commercial programs, although this is usually not a serious problem since the term $S_{ij}v_{k,k}$ is negligibly small for many materials and zero for incompressible materials.

The Cotter-Rivlin rate, which corresponds to $m=-2$, is the same as the Lie derivative, except that it again misses the volumetric work term $S_{ij}v_{k,k}$, and thus is not work-conjugate to any finite strain tensor. Neither are the Oldroyd rate and the Green-Naghdi rate.

Tangential Stiffness Moduli and Their Transformations to Achieve Energy Consistency

The tangential stress-strain relation has generally the form

$S_{ij}^{(m)}=C_{ijkl}^{(m)}{\dot {e}}_{kl}$ (9)

where $C_{ijkl}^{(m)}$ are the tangential moduli (components of a 4th-order tensor) associated with strain tensor $\epsilon _{ij}^{(m)}$. They are different for different choices of $m$, and are related as follows:

$\left[C_{ijkl}^{(m)}-C_{ijkl}^{(2)}-{\frac {1}{4}}(2-m)(S_{ik}\delta _{jl}+S_{jk}\delta _{il}+S_{il}\delta _{jk}+S_{jl}\delta _{ik})\right]v_{k,l}=0$ (10)

Noting that $C_{ijkl}^{(2)}v_{k,l}={\hat {S}}_{ij}$ and $v_{i,j}+v_{j,i}=2{\dot {e}}_{ij}$, and that $C_{ijkl}^{(m)}v_{k,l}=S_{ij}^{(m)}$, one can verify that Eq. (10) is identical to Eq. (9) with (5). Now, from the fact that Eq. (10) must hold true for any velocity gradient $v_{k,l}$, it follows that:

$C_{ijkl}^{(m)}=C_{ijkl}^{(2)}+(2-m)[S_{ik}\delta _{jl}]_{sym},[S_{ik}\delta _{jl}]_{sym}={\frac {1}{4}}(S_{ik}\delta _{jl}+S_{jk}\delta _{il}+S_{il}\delta _{jk}+S_{jl}\delta _{ik})$ (11)

where $C_{ijkl}^{(2)}$ are the tangential moduli associated with the Green-Lagrangian strain $(m=2)$, taken as a reference, $S_{ij}$ = current Cauchy stress, and $\delta _{ij}$ = Kronecker delta (unit tensor).

Eq. (11) can be used to convert a black-box commercial finite element program from one objective stress rate to another. To this end, the transformation in Eq. (11) must be used in the user's material subroutine in every load step at each integration point of each finite element. Since the loading step is always finite, the best convergence at step refinement occurs when $S_{ij}$ in Eq. (11) is taken approximately as the average of the initial and final values in the step, which requires iterations of each step. To convert an eigenvalue subroutine, and initial guess of the eigenvalue must be made and then, in each subsequent iteration, the $_{ijkl}^{(m)}$-values must be updated according to the new stresses obtained the previous eigenvalue.

Since $S_{ij}{\dot {e}}_{kk}=(S_{ij}\delta _{kl})\delta e_{kl}$, the transformation

$C_{ijkl}^{\mbox{conj}}=C_{ijkl}^{\mbox{nonconj}}+S_{ij}\delta _{kl}$ (12)

can further correct for the absence of the term $S_{ij}v_{k,k}$ (the term $S_{ij}\delta _{km}$ does not allow interchanging subscripts $ij$ and $kl$, which is what breaks the major symmetry of the tangential moduli tensor). This transformation (whose special case, without consideration of work-conjugacy, is implied in the relation to Truesdell rate derived in may be used in the user's material subroutine to convert a commercial program from a non-conjugate to a work-conjugate objective stress rate. It ensures that the results satisfy the energy balance (or the first law of thermodynamics), i.e., that ${\hat {S}}_{ij}\Delta t\delta \epsilon _{ij}$ be the correct expression for the second-order work per unit volume during the loading step. If this transformation is implemented in the user's material subroutine in every load step at each integration point of each finite element, and if also $C_{ijkl}^{(2)}$ is expressed in terms of $C_{ijkl}^{(0)}$ according to the transformation in Eq. (11) for $m=0$, any commercial program using the Jaumann rate of Cauchy stress can be converted to simulate a program using the Truesdell rate, a rate that is work-conjugate with the Green-Lagrangian strain ($m=2$).

Often $S_{ij}$ are so much smaller than $C_{ijkl}$ that the differences among different objective stress rates cause only negligible errors which do not matter in most engineering applications. There are nevertheless some important cases where this is not so. In, an example is given of an indentation of a thick sandwich plate with a polymeric (vinyl) foam core, in which the non-conjugacy error of a commercial program exceeds 30% of the reaction force predicted by commercial software. Generally, similar errors must be expected for metallic and ceramic foams, honeycomb, snow, loess, silt, clay, loose sand, organic soils, tuff, pumice, many insulation materials (including space shuttle tiles), lightweight concretes, light wood, carton, osteoporotic bone and various biologic tissues.

Large strain often develops when the material behavior becomes nonlinear, due to plasticity of damage. Then the primary cause of stress dependence of the tangential moduli is the physical behavior of material. What Eq. (11) means is that the nonlinear dependence of $C_{ijkl}$ on the stress must be different for different objective stress rates. Yet none of them is more convenient than another---except if there exists one stress rate, one $m$, for which the moduli can be considered constant. This exception occurs for soft-in-shear structures such as columns, plates, shells and lattices of high shear deformability.

The Special Case of Structures Very Soft in Shear

Examples of such structures are highly orthotropic fiber-polymer composites, homogenized sandwich plates or shells with foam or honeycomb cores, homogenized lattice columns, and layered bearings. Under axial or in-plane compressive loads, these structures typically reach the critical (or bifurcation) load of buckling while all the material is still at small strain, often a strain much smaller than the limit of linearity. So one would expect that constant elastic moduli, as measured in laboratory tests, would be appropriate. Not necessarily, though.

According to Eq. (11), the elastic moduli can, at small strain, be considered constant for only one $m$-value, i.e., for one objective strain rate only. But for which one?---For the Truesdell rate ($m=2$), if the the structure is highly orthotropic and compressed in the strong direction. Only for this rate, the work of in-plane forces on lateral deflections of a soft-in-shear plate can be captured correctl while the tangential moduli are kept constant. On the other hand, when the highly orthotropic body is compressed transversely to the strong direction, energy consistency requires the Lie derivative of Kirchhoff stress, corresponding to $m=-2$.

To explain it, consider the example in Fig. which shows a stocky soft-in-shear sandwich column compressed axially and centrically by force $P$, and imagined stages of the deformations of a short element of the column during buckling from initial position (a). Stage(b) shows a small rotation $\psi$ of the element without shear (equal to the slope $w'={\mbox{d}}w/{\mbox{d}}x$ of the deflection curve $w(x)$). Stage (c) shows subsequent small shearing by angle $\gamma$ while keeping the skin rotation constant ($\gamma =w'={\mbox{d}}w/{\mbox{d}}x$ of deflection curve $w(x)$). Stage (c) shows the the final shear deformation with small angle $\gamma$ (Fig. 2d). The work of $P$ due to rotation $\psi$ is $P(1-\cos w'){\mbox{d}}x\simeq P{w'}^{2}/2$ (with second-order accuracy in $w'$), and the strain energy of the deformed element is $bhG\gamma ^{2}/2$, where $G$ = elastic shear modulus of the sandwich core, $h$ = its thickness, and $b$ = thickness of the skins ($b\ll h$). So, when the buckling begins, the change of potential energy of the column per unit length is

$\delta ^{2}W=bhG\gamma ^{2}/2-P{w'}^{2}/2$ (13)

based on first principles. But $\delta ^{2}W$ can also be calculated from $C_{ijkl}$ using Eq. (11) and making the usual simplifying assumption that the cross section of core remains plane and the skin thickness is negligible. If consideration is restricted to homogeneous strain fields in the core and skins, the result is

$\delta ^{2}W=\left[bhG+P(2-m)/4\right]\gamma ^{2}/2-P{w'}^{2}/2$ (14)

Since this equation must be identical to Eq. (13), it is necessary that $m=2$, which indicates the Truesdell rate, associated with the Green-Lagrangian finite strain. Generally, this is the rate that is required for axially or in-plane loaded beams, plates and shells. When the energy of axial strains in the skins during the buckling is included, Eq. (13) leads to the classical Engesser formula for shear buckling of columns, which has been extensively verified by experiments.

When $m=-2$ is used in Eq. (14), one gets for shear buckling of columns the Haringx formula. This formula faithfully describes the buckling of layered elastomeric bearings used for bridge supports or for seismic isolation, in which horizontal steel plates are alternated with rubber layers. A similar argument as in Fig. 2, in which the rotation of column element is followed by the shearing of the horizontal rubber layer, shows that $m$ must be equal to $-2$.

For the critical load of buckling of short bulky columns, the difference between these two formulas, associated with $m=$ or $m=-2$, can exceed 100% even if the material is in the small-strain linear range. This discrepancy shows that choosing the correct objective stress rate can be extremely important.

In thin plates and shells in which the transverse shear strains are negligible compared to the axial or in-plane strains, the in-plane rotations and in-plane strains are negligible compared to the out-of-plane rotations. In that case the distinction among all finite strains, and thus also among all objective stress rates, disappears, and the incompressibility does not matter. All the rates are then equivalent, and one can use, e.g., Eq. (8). The same applies to slender columns without shear.

Another case where the differences between the objective stress rates for different $m$ are important occurs for the damage constitutive laws in which some of the tangential moduli $C_{ijkl}^{(m)}$ gradually decrease during loading and, on approach to peak load, become of the same order of magnitude as some stress components. Then the peak load indicated by a check for singularity of the matrix of tangential moduli can be rather different for different $m$. In the constitutive models for quasibrittle materials, capturing the damage due to distributed parallel splitting cracks, the material can become incrementally strongly orthotropic, with a very low shear-to-normal stiffness ratio. Big differences due to the choice of $m$, similar to those mentioned above for shear-buckling of columns and plates, may then be encountered.

The main reason for preferring the Hencky (logarithmic) strain has been the additivity of subsequent strains and the tension-compression symmetry, features that are convenient for generalizing small-strain constitutive laws to large strains. However, it is no problem to use one strain measure in developing a constitutive law and then make a conversion to another strain measure for the purpose of structural analysis.

References

• ACI Committee 209 (1972) "Prediction of creep, shrinkage and temperature effects in concrete structures" ACI-SP27, Designing for Effects of Creep, Shrinkage and Temperature}, Detroit, pp. 51--93 (reaproved 2008)
• ACI Committee 209 (2008). Guide for Modeling and Calculating Shrinkage and Creep in Hardened Concrete ACI Report 209.2R-08, Farmington Hills.
• Brooks, J.J. (2005). "30-year creep and shrinkage of concrete." Magazine of Concrete Research, 57(9), 545--556. Paris, France.
• CEB-FIP Model Code 1990. Model Code for Concrete Structures. Thomas Telford Services Ltd., London, Great Britain; also published by Comité euro-international du béton (CEB), Bulletins d'Information No. 213 and 214, Lausanne, Switzerland.
• fib Model Code 2011. "Fédération internationale de béton (fib). Lausanne.
• Harboe, E.M., et al. (1958). "A comparison of the instantaneous and the sustained modulus of elasticity of concrete", Concr. Lab. Rep. No. C-354, Division of Engineering Laboratories, US Dept. of the Interior, Bureau of Reclamation, Denver, Colorado.
• Jirásek, M., and Bažant, Z.P. (2001). Inelastic analysis of structures, J. Wiley, London (chapters 27, 28).
• RILEM (1988a) Committee TC 69, Chapters 2 and 3 in Mathematical Modeling of Creep and Shrinkage of Concrete, Z.P. Bažant, ed., J. Wiley, Chichester and New York, 1988, 57--215.
• Troxell, G.E., Raphael, J.E. and Davis, R.W. (1958). "Long-time creep and shrinkage tests of plain and reinforced concrete" Proc. ASTM 58} pp. 1101-1120.
• Vítek, J.L., "Long-Term Deflections of Large Prestressed Concrete Bridges, 20AC CEB Bulletin d20AC%u2122Information No. 23520AC%u201CServiceability Models20AC%u201C Behaviour and Modelling in Serviceability Limit States Including Repeated and Sustained Load, CEB, Lausanne, 1997, pp. 215-227 and 245-265.
• Wittmann, F.H. (1982). "Creep and shrinkage mechanisms." Creep and shrinkage of concrete structures, Z.P. Bažant and F.H. Wittmann, eds., J. Wiley, London 129--161.
• Bažant, Z.P., and Yu, Q. (2012). "Excessive long-time deflections of prestressed box girders." ASCE J. of Structural Engrg. 138 (6), 676-686, 687-696.