Work conjugacy error in commercial finite-element codes: its magnitude and how to compensate for it

Zdeněk P. Bažant, Mahendra Gattu, Jan Vorel


Most commercial finite-element programs use the Jaumann (or co-rotational) rate of Cauchy stress in their incremental (Riks) updated Lagrangian loading procedure. This rate was shown long ago not to be work-conjugate with the Hencky (logarithmic) finite strain tensor used in these programs, nor with any other finite strain tensor. The lack of work-conjugacy has been either overlooked or believed to cause only negligible errors. Presented are examples of indentation of a naval-type sandwich plate with a polymeric foam core, in which the error can reach 28.8 per cent in the load and 15.3 per cent in the work of load (relative to uncorrected results). Generally, similar errors must be expected for all highly compressible materials, such as metallic and ceramic foams, honeycomb, loess, silt, organic soils, pumice, tuff, osteoporotic bone, light wood, carton and various biological tissues. It is shown that a previously derived equation relating the tangential moduli tensors associated with the Jaumann rates of Cauchy and Kirchhoff stresses can be used in the user’s material subroutine of a black-box commercial program to cancel the error due to the lack of work-conjugacy and make the program perform exactly as if the Jaumann rate of Kirchhoff stress, which is work-conjugate, were used.

1. Problem

The finite-element analysis of inelastic solids subjected to finite strain is generally conducted by an incremental (or Riks) updated Lagrangian procedure in which the material constitutive properties are in each small loading step characterized in a rate (or tangential) form as follows: Embedded Image1.1Here the subscripts refer to Cartesian coordinates xi (i=1,2,3); Embedded Image rate of the small (linearized) strain Embedded Image, Cijkl=tangential moduli, Embedded Image stress rate and repetitions of subscripts imply summation.

There exist many equally valid objective stress rates associated by work with various choices of the incremental finite strain tensor. Most commercial programs use the Hencky (or logarithmic) finite strain tensor Embedded Image, for which the associated (or work-conjugate) objective stress rate is the Jaumann rate of Kirchhoff stress, expressed as Embedded Image1.2where Sij is the current Cauchy (or true) stress, Embedded Image is the material rotation rate (or spin tensor), vi is the material point velocity and vk,k represents the rate of change of volume.

However, most commercial programs use, instead of Embedded Image, the Jaumann (or co-rotational) rate of Cauchy (or true) stress (Bažant & Cedolin 1991; Belytschko et al. 2000), Embedded Image1.3This rate misses the last, volumetric, term in equation (1.2). It has been shown (Bažant 1971; Bažant & Cedolin 1991) that Embedded Image is work-conjugate to the Hencky finite strain tensor Embedded Image, but Embedded Image is not, and in fact there exists no finite strain tensor to which it would be work-conjugate. This means that the error due to the missing volumetric term Sijvk,k violates energy balance, i.e. the first law of thermodynamics (the reason for superscript (0) is that the Hencky finite strain tensor is a special case for m→0 of the Doyle–Ericksen finite strain tensors with parameter m).

The error in using the Jaumann rate of Cauchy stress was pointed out in the literature long ago, are but has either been ignored or thought to cause only negligible errors. Indeed, the error is generally less than 0.1 per cent for metals and other materials that are inelastically incompressible. It now appears, though, that this is not true for highly compressible inelastic materials; for example, rigid foams (polymeric, metallic and ceramic), honeycomb, certain soils (loess, silt, under-consolidated and organic soils), some rocks (pumice and tuff), osteoporotic bone, light wood, carton and various biological tissues. The above-mentioned problem afflicts the most widely used commercial software—ABAQUS, ANSYS, LS-DYNA and some versions of NASTRAN.

2. Remedy

The best remedy for the software makers would be to change the black-box coding of these programs. This is the only way when built-in material models are used. However, most commercial programs feature a user’s material subroutine (in ABAQUS called UMAT), in which the user can define the constitutive law by specifying Cijkl as a function of the state variables. In that case, a remedy can easily be made by the user.

Substituting Embedded Image into equation (1.1) and noting that Embedded Image, one gets Embedded Image. This may be rearranged as Embedded Image. Consequently, the change of the objective stress rate from equations (1.3) to (1.2) is equivalent to the following transformation of the tangential moduli (Bažant & Cedolin 1991, eqn. 11.4.6): Embedded Image2.1This conjugacy transformation, which restores major symmetry of the fourth-order tensor (i.e. the symmetry at interchange of subscripts ij and kl), ensures that the results satisfy the energy balance, or the first law of thermodynamics, i.e. that Embedded Image is the correct expression for the second-order work per unit volume. A more detailed and general derivation of this transformation was presented in eqn 11.4.6 (on p. 727) of Bažant & Cedolin (1991). A special case of it was implied in a transformation presented in Hibbitt et al. (1970) (with no attention to work-conjugacy). However, it was not realized that this transformation could avoid huge errors in some situations.

The reason that the conjugacy transformation in equation (2.1) does not break the major symmetry of the tangential moduli of material, and the symmetry of the structural stiffness matrix, is that this transformation is specified in UMAT simultaneously with the given constitutive law. The combination of both is thus handled by ABAQUS as if it were the constitutive law. A symmetric stiffness matrix is then generated for this combined constitutive law.

The user’s material subroutine is written under the assumption that Embedded Image is specified in this subroutine. But what the user correctly needs to use is Embedded Image. Therefore, transforming the tangential moduli specified by the constitutive law according to equation (2.1) will cause the commercial program to deliver the correct results, preserving energy balance. This transformation must be implemented in the user’s material subroutine in every load step at each integration point of each finite element.

3. An example of indentation of a foam-core sandwich for a light ship hull

This problem is of interest for prospective new designs of very light and fuel-efficient large ships in which the hull is clad by sandwich panels that have a foam core approximately 150 mm thick and that are supported on tubular steel ribs. One important design consideration is the impact of floating objects or ice floes onto the sandwich.

Two sets of case studies are conducted. In the first case study, a sandwich panel of given geometry and material properties intended for practical application is analysed to determine the work of indentation and the relationship of load and indentation depth. In the second case study, a sandwich panel of the same core thickness but different material properties is studied to document a possible magnitude of error.

Case study 1. The sandwich shown in figure 1 is indented in a two-dimensional plane-strain mode by a rigid rectangular object of width b=160 mm. The bottom skin is assumed to be fixed on a rigid base and to extend far enough so that the boundary condition at the rims do not matter. The object is assumed to be in perfect contact with the top skin.

Figure 1.

(a) Undeformed configuration of a sandwich plate with finite-element mesh; (b) deformed configuration with deformed mesh; (c) contours of equal volumetric strain.

The core consists of a vinyl foam of mass density 200 kg m−3, marketed as Divinycell 200, typically used for marine sandwich structures. The core thickness is tc=150 mm and the skin thickness is ts=9.6 mm. The skin is elastic and quasi-isotropic, with Es=46 GPa and νs=0.3.

The constitutive properties of the foam are simplified as follows: von Mises (or J2) plasticity in shear, with shear yield stress τ0=62 MPa, and the volumetric compression defined by the tri-linear diagram plotted in figure 2. The middle, nearly flat, portion of this diagram corresponds to collapse of foam cells (modelled in Brocca et al. (2001)). The final stiffening from initial bulk modulus K to bulk modulus γK is caused by the closing of collapsed cells with the opposite walls pressed into contact. Young’s elastic modulus of the foam core is Ec=230 MPa, and the elastic Poisson ratio is νc=0.353. The parameters σa=−2 MPa, σb=−2.1 MPa, ϵa=σa/3K, ϵb=σb/3K, α=(σaσb)/(ϵaϵb)=0.02 and γ=1.5.

Figure 2.

(a) Volumetric stress–strain relation of the core; (b) deviatoric (von Mises) yield surface of the core. (Online version in colour.)

Because the objective is to reveal the differences between the predictions of the corrected and uncorrected code, rather than to obtain a precise prediction of the indentation depth, the straight line simplifications of the stress–strain relation should not seriously affect the comparisons. Nevertheless, the stress–strain diagram in figure 2 is similar to those identified from foam crushing experiments (Shuaeib & Soden 1997; Deshpande & Fleck 1999; Tagarielli et al. 2005). The linear response to volumetric expansion is also a simplification, but it is irrelevant because no finite elements experienced such expansion. Another simplification is the decoupling of the deviatoric plastic strain from the volumetric compression, but previous experiments give no indication of any significant coupling of this kind.

The sandwich is discretized by the mesh shown in figure 1a, which contains 6422 four-node plane-strain elements (of the type CPE4). At the end of loading the mesh gets heavily distorted, which may, of course, introduce a non-negligible error. However, such an error cannot not matter for the comparison of the corrected and uncorrected forces P at the same u (figures 3 and 4), because the mesh distortion for the same u is about the same.

Figure 3.

Load–displacement curves for case study 1. (Online version in colour.)

Figure 4.

(a) Curve of secant stiffness versus displacement (for E=500). (b) Load–displacement curves obtained in case study 2. Solid lines, ABAQUS corrected; dashed lines, ABAQUS. (Online version in colour.)

The analysis is carried out by applying incremental displacements Δu=1.2 mm and 0.12 mm using the full Newton iterative scheme implemented in ABAQUS (2009). The total reaction force (or load P) is computed with ABAQUS. Two cases are run: (i) with and (ii) without the transformation of tangential moduli according to equation (2.1).

The resulting load deflection curves P(u) are plotted in figure 3. The difference in the maximum load reached is 17.67 per cent of the uncorrected ABAQUS value. The load versus displacement (Pu) curves are visually indistinguishable for Δu=1.2 mm and 0.12 mm, as shown in figure 3. This illustrates that the forward Euler solution has converged for these load increments.

The area under the load-deflection curve P(u), representing the work of load, is Embedded Image. The area between the curve computed for the uncorrected moduli (the default case) and the curve for the corrected moduli is the error in energy. This error is 12.07 per cent of the uncorrected ABAQUS value of work.

Case study 2. The second case study deals with a sandwich panel whose skin is elastic and orthotropic, made of unidirectional carbon-epoxy laminate AS-4/3501-6. Its elastic properties are characterized by E1=E3=146.6 GPa, E2=72.4 GPa, G12=G32=7.6 GPa, ν12=ν32=0.28, ν21=ν23=0.02, G13=56.39 GPa and ν13=0.3, where subscripts 1 and 2 label the longitudinal and transverse directions. The thickness of the core, tc=150 mm, and impacting object width, b=160 mm, are the same as mentioned before. The skin thickness is ts=6 mm.

The sandwich is discretized by a similar mesh as mentioned before, having 6422 four-node plane-strain elements (of the type CPE4). For Young’s modulus of the core, three values are considered: Ec=500 MPa, 230 MPa and 120 MPa, with the corresponding yield stress values τ0=184 MPa, 90 MPa and 49 MPa, respectively. Furthermore, νc=0, σa=−2 MPa, σb=−2.1 MPa, ϵa=σa/3K, ϵb=σb/3K, α=(σaσb)/(ϵaϵb)=0.02, γ=1.5. The volumetric stress–strain relation of the core is shown in figure 5. The resulting load-deflection curves P(u) are plotted in figure 4. Numerically, the results are summarized in table 1.

Figure 5.

Volumetric stress–strain relations for the core with various E-values, considered in case study 2. (Online version in colour.)

View this table:
Table 1.

Difference between the results of uncorrected and corrected ABAQUS, relative to the uncorrected values.

It may now be noted that, without the correction in equation (2.1), the maximum relative error ΔP, error in reaction P, and the maximum relative error in energy, ΔW, error, as computed by ABAQUS, are Embedded Image3.1relative to the uncorrected ABAQUS results. These errors occur for E=500 MPa at the end of computation, at which the core is compressed to 80 per cent of its original thickness. If the displacement were incremented further, still higher percentage errors would be obtained.

These errors are certainly severe enough not to be ignored. Similar corrections must be made in other commercial programs in which the incremental analysis is based on the Jaumann rate of Cauchy stress. They include LS-DYNA, ANSYS and some versions of NASTRAN (on the other hand, the commercial code ATENA and the open-source code OOFEM use the Truesdell rate, and thus require no correction).

There is also the issue of numerical error due to finiteness of the load steps. Ideally, Sij in the transformation equation (2.1) is the current Cauchy (true) stress tensor. But for a finite load step Δu, the value of this tensor is available only for the beginning of the current step (or the end of the previous step). This is the value that has been used in each load step in the computations for figure 3. However, using these delayed stress values represents the forward Euler integration method, which must introduce a certain numerical error. The accuracy may be poor if the load steps are not small enough.

For a central difference scheme, the accuracy would be better and the convergence upon step refinement would become quadratic rather than linear.

To verify the accuracy of the forward Euler method, the response of a single four-node plane-strain element (CPE4) was computed. The four nodes of this element lie at the vertices of a square of side 1 mm (figure 6a). The element is compressed by applying vertical displacement increments Δu, for which three values are used: Δu=0.01 mm, 0.0005 mm and 0.000125 mm, at no lateral displacements. The simulations are run in ABAQUS again, both with and without correction. As seen in figure 6b, the step of 0.01 mm is too large because it produces a large error. But the steps of 0.0005 mm and 0.000125 mm are sufficient because the corresponding responses in figure 6b are visually indistinguishable. This confirms convergence of the forward Euler method as load steps are refined.

Figure 6.

(a) Square four-node finite element CPE4; (b) P(u)-curves from corrected and uncorrected ABAQUS obtained for loading steps Δu=0.01 mm, 0.0005 mm and 0.000125 mm (the last two are visually indistinguishable). Solid lines, ABAQUS corrected; dashed lines, ABAQUS. (Online version in colour.)

4. Broader aspects and related issues

Sandwich shell indentations are also of interest for aircraft. They often happen accidentally during maintenance. A sandwich must be sufficiently tolerant of such accidents. Another problem where a similar kind of error takes place is the delamination fracture of the skin, which cannot occur without significant volume changes in the foam core.

The definition of the finite strain measure ϵij is a matter of choice, though within some restrictions. For each chosen tensor, there exists one and only one work-conjugate objective stress rate and tangential stiffness tensor (Bažant 1971). The work-conjugacy is properly established by variational energy analysis, which is presented in detail in Bažant & Cedolin (1991, §§11.1–11.5) (in brief, see the appendix).

The stress dependence of the tangential stiffness tensor Cijkl is different for each conjugate finite strain tensor. By introducing in UMAT a transformation of Cijkl that is linear in current Cauchy stress Sij, one can make the commercial program such as ABAQUS simulate the response for any chosen finite strain measure (see Bažant & Cedolin 1991, eqn. 11.4.5). Generally, the dependence of Cijkl on stress is a manifestation of material nonlinearity, and for each chosen finite strain measure, the dependence on stress is different.

An interesting problem that was for a long time regarded as perplexing arises in buckling of highly orthotropic plates and sandwich plates very soft in shear. They can buckle under in-plane compressive forces, while the strain is everywhere still small compared with the material strength limit and lies within the linear range of the material. This means that the elastic moduli of the material should be considered as constant. But because of the above-mentioned linear transformation Cijkl, the constancy is possible only for one finite strain tensor. Which one?

For in-plane loaded sandwich plates or soft-in-shear orthotropic plates, it is the Green–Lagrangian finite strain tensor corresponding to m=2 (see Bažant & Beghini (2005, 2006) and the examples in Ji & Waas (2009); Ji et al. (2010)). To be able to use a commercial program such as ABAQUS, the transformation in (Bažant & Cedolin 1991, eqn. 11.4.5) (or equation (A5) in the appendix) needs to be superposed on that in equation (2.1).

For transverse compression of a foam layer, a transformation to m=2 is not necessary, provided that the measurements of the nonlinear foam properties have been interpreted in figure 2 on the basis of the Hencky strain tensor.

5. Conclusions

  • (i) The lack of work-conjugacy of the Jaumann rate of Cauchy stress, which was pointed out long ago but has generally been ignored owing to the allegation that the resulting errors are negligible, is found to cause major errors in analysing structures consisting of highly compressible materials. In a practical example of indentation of a naval-type foam-core sandwich plate, the error in the force and energy required for indentation is shown to reach 28.8 per cent and 15.3 per cent, respectively.

  • (ii) The demonstrated error magnitude has important implications for the developers of commercial finite-element programs as well as their users. The objective stress rate used in these programs should be replaced by a work-conjugate rate.

  • (iii) But even if this replacement is not done, the user can make a simple correction to eliminate the non-conjugacy error. A previously derived relation between the tangential moduli tensors of the Jaumann rates of Cauchy stress and of Kirchhoff stress may be used in the user’s material subroutine of a black-box commercial program to modify the tangential moduli obtained from the constitutive law.

  • (iv) Although this correction should ideally be based on the mean Cauchy stress in the current load step, this stress value is unavailable. So the correction must be based on the stress obtained for the end of the previous step. Nevertheless, the numerical integration error caused by delaying the stress value is, for normal load step size, negligible.


Financial support under grant no. N00014-11-1-0155 from the Office of Naval Research to Northwestern University, monitored by Dr Roshdy S. Barsoum, is gratefully appreciated.

Appendix A. Review of energy-consistent objective stress rates

The objective stress rates are usually derived by tensorial coordinate transformations. However, to reveal work-conjugacy, it is better to use a variational energy approach (Bažant 1971). Consider incremental finite strain tensors ϵij relative to the initial (stressed) state at the beginning of the load step, using the initial (Lagrangian) coordinates xi (i=1,2,3) of material points. A broad class of equally admissible finite strain tensors is represented by the Doyle–Ericksen tensors whose second-order approximation is Embedded ImageA1The case m=2 gives the Green–Lagrangian strain tensor, m=1 gives the Biot strain tensor, m=0 gives the Hencky (logarithmic) strain tensor and m=−2 gives the Almansi–Lagrangian strain tensor. The work δW done at small deformations of a material element of unit initial volume, starting from an initial state under initial Cauchy (or true) stress Embedded Image, can be expressed in two equivalent ways, Embedded ImageA2where Embedded Image represents an arbitrary strain variation, and Embedded Image a small stress incre- ment—which is symmetric and objective (an incremental second Piola–Kirchhoff stress); and τij a non-symmetric small Lagrangian (or first Piola–Kirchhoff) stress increment. Since the first-order work, Embedded Image is cancelled in the virtual work equation of equilibrium by the work of loads, only the second-order work is of interest.

The two expressions in equation (A2) must be equal. Impose this equality and substitute Embedded Image (by virtue of symmetry of Sij), Embedded Image (which suffices for second-order work accuracy in ui,j), Embedded Image and Embedded Image (where vi,jΔt=δui,j, and Embedded Image). Then introduce the variational condition that the resulting equation must be valid for any δui,j. This yields (Bažant 1971; Bažant & Cedolin 1991) Embedded ImageA3where Embedded Image, Embedded Image, and Embedded Image rate of Cauchy stress. Evaluating equation (A3) for general m and for m=2, one gets a general expression for the objective stress rate (Bažant 1971; Bažant & Cedolin 1991), Embedded ImageA4where Embedded Image rate. In particular, for m=0, this yields equation (1.2) for the Jaumann rate of Kirchhoff stress, which confirms its work-conjugacy to the Hencky strain. For m=2, equation (A3) gives the Truesdell rate, for m=1 the Biot rate and for m=−2 the Lie derivative. The Jaumann (co-rotational) rate of Cauchy stress, equation (1.3), cannot be obtained from equation (A3) and thus is not work-conjugate with any finite strain tensor (the same is true for the Cotter–Rivlin, Oldroyd and Green–Naghdi objective stress rates).

When different m are considered, the tangential stress–strain relation in equation (1.1) must be rewritten as Embedded Image, where moduli Embedded Image are associated with strain tensor Embedded Image. They are different for different choices of m, and are related as follows (Bažant 1971; Bažant & Cedolin 1991): Embedded ImageA5and Embedded ImageA6Here, Embedded Image are the tangential moduli associated with the Green–Lagrangian strain (m=2), taken as a reference, Sij is the current Cauchy stress, and δij is the Kronecker delta. Using equation (A5), one can convert a black-box commercial finite-element program from one objective stress rate to another, in the same way as shown here for the conversion to work-conjugacy.


  • On leave from Czech Technical University in Prague.

  • Received March 17, 2012.
  • Accepted April 26, 2012.


View Abstract