Review of energy conservation errors in finite element softwares caused by using energy-inconsistent objective stress rates

The paper briefly summarizes the theoretical derivation of the objective stress rates that are work-conjugate to various finite strain tensors, and then briefly reviews several practical examples demonstrating large errors that can be used by energy inconsistent stress rates. It is concluded that the software makers should switch to the Truesdell objective stress rate, which is work-conjugate to Green’s Lagrangian finite strain tensor. The Jaumann rate of Cauchy stress and the Green-Naghdi rate, currently used in most software, should be abandoned since they are not work-conjugate to any finite strain tensor. The Jaumann rate of Kirchhoff stress is work-conjugate to the Hencky logarithmic strain tensor but, because of an energy inconsistency in the work of initial stresses, can lead to severe errors in the cases of high natural orthotropy or strain-induced incremental orthotropy due to material damage. If the commercial softwares are not revised, the user still can make in the user’s implicit or explicit material subroutines (such as UMAT and VUMAT in ABAQUS) a simple transformation of the incremental constitutive relation to the Truesdell rate, and the commercial software then delivers energy consistent results.     

On the Consequences of the Adoption of the Zaremba–Jaumann Objective Stress Rate in FEM Codes

This paper deals with a particular issue of computational mechanics in main FEM codes nowadays available, i.e. the outcomes of implementations of large strain constitutive models based on the adoption of so-called objective stress rates, in order to satisfy objectivity requirements. The point here is that of directly inquiring whether well-known incoherencies due to the adoption of the Zaremba–Jaumann objective stress rate may manifest themselves when the most used elastic and elastoplastic constitutive models are adopted. The present investigation aims at providing a comprehensive review of the theoretical aspects and at developing an informed knowledge to final users of FEM codes, in terms of exposing which constitutive models and FEM implementations may be affected by Zaremba–Jaumann objective stress rate induced incoherencies. Towards this end, local FEM simple shear tests are explored and clearly show that kinematic cases characterized by a non zero spin may be heavily affected by oscillatory incoherencies, which arise for expected cases, i.e. Cauchy stress responses, but also for other less expected cases, i.e. strain responses, whether they are total, elastic or plastic. Beyond local tests, structural simple shear tests are also performed and show as well that oscillatory incoherencies found in local simple shear tests may heavily influence the overall structural outcomes. A non-secondary target of the paper is that of reviewing the relevant scientific and technical literature about objective stress rates, by critically analyzing correlated issues and proposed solutions, considering scientific contributions spanning over a century, keeping specific attention to the treatment of the Zaremba–Jaumann objective stress rate and to the possible flaws related to its adoption.     

Geometrically nonlinear formulation for the axi-symmetric transition finite elements

This paper presents a geometrically nonlinear formulation for the axi-symmetric transition finite elements using total lagrangian approach. The basic element is formulated using properties of the axi-symmetric solids and the axi-symmetric shells. A novel feature of the formulation presented here is that the restriction on the magnitude of the rotations for the shell nodes of the transition element is eliminated. This is accomplished by retaining true nonlinear functions of nodal rotations in the definition of the element displacement field. Such transition elements are essential for geometrically nonlinear applications requiring both axi-symmetric solids and the axi-symmetric shells. They ensure proper connection of the axi-symmetric solid portion of the structure to the shell like portion of the structure. It is shown that the selection of different stress and strain components at the integration points does not effect the overall linear response of the element. However, in the geometrically nonlinear formulation, it is necessary to select appropriate components of the stresses and the strains at the integration point for accurate and converging element behavior. Numerical examples are presented to demonstrate such characteristics of the transition elements.     

A nonlinear visco-elastic constitutive model for large rotation finite element formulations

Although all known materials have internal damping that leads to energy dissipation, most existing large deformation visco-elastic finite element formulations are based on linear constitutive models or on nonlinear constitutive models that can be used in the framework of an incremental co-rotational finite element solution procedure. In this investigation, a new nonlinear objective visco-elastic constitutive model that can be implemented in non-incremental large rotation and large deformation finite element formulations is developed. This new model is based on developing a simple linear relationship between the damping forces and the rates of deformation vector gradients. The deformation vector gradients can be defined using the decomposition of the matrix of position vector gradients. In this paper, the decomposition associated with the use of the tangent frame that is equivalent to the QR decomposition is employed to define the matrix of deformation gradients that enter into the formulation of the viso-elastic constitutive model developed in this investigation. Using the relationship between the deformation gradients and the components of the Green–Lagrange strain tensor, it is shown that the damping forces depend nonlinearly on the strains and linearly on the classical strain rates. The relationship between the damping forces and strains and their rates is used to develop a new visco-elastic model that satisfies the objectivity requirements and leads to zero strain rates under an arbitrary rigid body displacement. The linear visco-elastic Kelvin–Voigt model frequently used in the literature can be obtained as a special case of the proposed nonlinear model when only two visco-elastic coefficients are used. As demonstrated in this paper, the use of two visco-elastic coefficients only leads to viscous coupling between the deformation gradients. The model developed in this investigation can be used in the framework of large deformation and large rotation non-incremental solution procedure without the need for using existing co-rotational finite element formulations. The finite element absolute nodal coordinate formulation (ANCF) that allows for straightforward implementation of general constitutive material models is used in the validation of the proposed visco-elastic model. A comparison with the linear visco-elastic model is also made in this study. The results obtained in this investigation show that there is a good agreement between the solutions obtained using the proposed nonlinear model and the linear model in the case of small deformations.     

Geometrically exact physics-based modeling and computer animation of highly flexible 1D mechanical systems

This paper presents a geometrically exact beam theory and a corresponding displacement-based finite-element model for modeling, analysis and natural-looking animation of highly flexible beam components of multibody systems undergoing huge static/dynamic rigid-elastic deformations. The beam theory fully accounts for geometric nonlinearities and initial curvatures by using Jaumann strains, concepts of local displacements and orthogonal virtual rotations, and three Euler angles to exactly describe the coordinate transformation between the undeformed and deformed configurations. To demonstrate the accuracy and capability of this nonlinear beam element, nonlinear static and dynamic analysis of two highly flexible beams are performed, including the twisting a circular ring into three small rings and the spinup of a flexible helicopter rotor blade (Graphical abstract). These numerical results reveal that the proposed nonlinear beam element is accurate and versatile for modeling, analysis and 3D rendering and animation of multibody systems with highly flexible beam components.     

Isogeometric configuration design optimization of shape memory polymer curved beam structures for extremal negative Poisson’s ratio

Using a continuum-based design sensitivity analysis (DSA) method, a configuration design optimization method is developed for curved Kirchhoff beams with shape memory polymers (SMP), from which we systematically synthesize lattice structures achieving target negative Poisson’s ratio. A SMP phenomenological constitutive model for small strains is utilized. A Jaumann strain, based on the geometrically exact beam theory, is additively decomposed into elastic, stored, and thermal parts. Non-homogeneous displacement boundary conditions are employed to impose mechanical loadings. At each equilibrium configuration, an additional nonlinear analysis is performed to calculate the Poisson’s ratio and its design sensitivity of the SMP material. The design objectives are twofold: for purely elastic materials, lattice structures are designed to achieve prescribed Poisson’s ratios under finite compressive deformations. Also, SMP-based lattice structures are synthesized to possess target Poisson’s ratios in specified temperature ranges. The analytical design sensitivity of the Poisson’s ratio is verified through comparison with finite difference sensitivity. Several configuration design optimization examples are demonstrated.     

Choice of objective rate in single parameter hypoelastic deformation cycles

Hill [Hill R. The mathematical theory of plasticity. Oxford: Clarendon Press; 1950] demonstrated that “infinitesimal-displacement theory, used in classical elastoplasticity, may no longer be valid in elastic–plastic analysis because the convected terms in the rate of change of the stress acting on material particle may then not be negligible” [Lee EH. Some anomalies in the structure of elastic–plastic theory at finite strain. In: Carroll MM, Hayes M. editors. Nonlinear effects in fluids and solids, New York: Plenum Press; 1996. p. 227–49]. From this we may deduce that

• •the elastic deformation part may have considerable influence on the total deformation, even when it is relatively small, i.e. we are confronted with the vital requirement of properly computing elastic deformations.
• •Further, finite deformation kinematics should be applied, i.e. they should take account of possibly large rotations, e.g. through a formulation in an Eulerian frame.

Xiao et al. [Xiao H, Bruhns OT, Meyers ATM. Self-consistent Eulerian rate type elasto-plasticity models based upon the logarithmic stress rate. Int J Plast 1999;15:479–520] gave the mathematical proof, that Bernstein’s consistency criterion [Bernstein B. Relation between hypo-elasticity and elasticity. Trans Soc Rheol 1960;4:23–8; Bernstein B. Hypoelasticity and elasticity. Arch Ration Mech Anal 1960;6:90–104] is fulfilled in a hypoelastic law of grade zero if, and only if, the objective logarithmic stress rate [Xiao H, Bruhns OT, Meyers A. Hypo-elasticity model based upon the logarithmic stress rate. J Elasticity 1997;47:51–68] has been applied. This proof is of complicated mathematical nature. Here, we compare several objective Eulerian stress rates of corotational and non-corotational type for the hypoelastic law cited above in closed single parameter deformation cycles. It is found that the logarithmic stress rate returns the element to its stress-free original state after the closed cycle, thus confirming the findings in Xiao et al. (1999). We show that for some other objective rates the errors are accumulating to considerable amounts after several cycles, even when the deformation in investigation is relatively small. Interestingly, for Jaumann stress rate, the error may vanish for specified deformation measures.     

A unified formulation of small-strain corotational finite elements: I. Theory

This paper presents a unified theoretical framework for the corotational (CR) formulation of finite elements in geometrically nonlinear structural analysis. The key assumptions behind CR are: (i) strains from a corotated configuration are small while (ii) the magnitude of rotations from a base configuration is not restricted. Following a historical outline the basic steps of the element independent CR formulation are presented. The element internal force and consistent tangent stiffness matrix are derived by taking variations of the internal energy with respect to nodal freedoms. It is shown that this framework permits the derivation of a set of CR variants through selective simplifications. This set includes some previously used by other investigators. The different variants are compared with respect to a set of desirable qualities, including self-equilibrium in the deformed configuration, tangent stiffness consistency, invariance, symmetrizability, and element independence. We discuss the main benefits of the CR formulation as well as its modeling limitations.     

Modeling large strain anisotropic elasto-plasticity with logarithmic strain and stress measures

In this paper we present a model and a fully implicit algorithm for large strain anisotropic elasto-plasticity with mixed hardening in which the elastic anisotropy is taken into account. The formulation is developed using hyperelasticity in terms of logarithmic strains, the multiplicative decomposition of the deformation gradient into an elastic and a plastic part, and the exponential mapping. The novelty in the computational procedure is that it retains the conceptual simplicity of the large strain isotropic elasto-plastic algorithms based on the same ingredients. The plastic correction is performed using a standard small strain procedure in which the stresses are interpreted as generalized Kirchhoff stresses and the strains as logarithmic strains, and the large strain kinematics is reduced to a geometric pre- and post-processor. The procedure is independent of the specified yield function and type of hardening used, and for isotropic elasticity, the algorithm of Eterovi? and Bathe is automatically recovered as a special case. The results of some illustrative finite element solutions are given in order to demonstrate the capabilities of the algorithm.     

On instabilities in large deformation simple shear loading

We obtain finite strain solutions for an elastic-perfectly-plastic isotropic material under simple shear loading. The effects of using the Jaumann, Green-Naghdi, and Truesdell stress rates are examined. An analytic solution is obtained using the Jaumann rate, and numerical solutions are obtained for the Green-Naghdi and Truesdell rates. The shear stress, as a function of the shear strain is nonmonotonic, and consequently ‘unstable’, for all three stress rates, which is contrary to previously published results.     

An extension of 3-D procedure to large strain analysis of shells

A numerical stress integration procedure for general 3-D large strain problems in inelasticity, based on the total formulation and the governing parameter method (GPM), is extended to shell analysis. The multiplicative decomposition of the deformation gradient is adopted with the evaluation of the deformation gradient practically in the same way as in a general 3-D material deformation. The calculated trial elastic logarithmic strains are transformed to the local shell Cartesian coordinate system and the stress integration is performed according to the GPM developed for small strain conditions. The consistent tangent matrix is calculated as in case of small strain deformation and then transformed to the global coordinate system.A specific step in the proposed procedure is the updating of the left elastic Green–Lagrangian deformation tensor. Namely, after the stresses are computed, the principal elastic strains and the principal vectors corresponding to the stresses at the end of time step are determined. In this way the shell conditions are taken into account appropriately for the next step.Some details are given for the stress integration in case of thermoplastic and creep material model.Numerical examples include bulging of plate (plastic, thermoplastic, and creep models for metal) and necking of a thin sheet. Comparison of solutions with those available in the literature, and with solutions using other type of finite elements, demonstrates applicability, efficiency and accuracy of the proposed procedure.     

Fast Corotated FEM using Operator Splitting

In this paper we present a novel operator splitting approach for corotated FEM simulations. The deformation energy of the corotated linear material model consists of two additive terms. The first term models stretching in the individual spatial directions and the second term describes resistance to volume changes. By formulating the backward Euler time integration scheme as an optimization problem, we show that the first term is invariant to rotations. This allows us to use an operator splitting approach and to solve both terms individually with different numerical methods. The stretching part is solved accurately with an optimization integrator, which can be done very efficiently because the system matrix is constant over time such that its Cholesky factorization can be precomputed. The volume term is solved approximately by using the compliant constraints method and Gauss‐Seidel iterations. Further, we introduce the analytic polar decomposition which allows us to speed up the extraction of the rotational part of the deformation gradient and to recover inverted elements. Finally, this results in an extremely fast and robust simulation method with high visual quality that outperforms standard corotated FEMs by more than two orders of magnitude and even the fast but inaccurate PBD and shape matching methods by more than one order of magnitude without having their typical drawbacks. This enables a very efficient simulation of complex scenes containing more than a million elements.     

Dynamics modeling for a rigid-flexible coupling system with nonlinear deformation field

In this paper, a moving flexible beam, which incorporates the effect of the geometrically nonlinear kinematics of deformation, is investigated. Considering the second-order coupling terms of deformation in the longitudinal and transverse deflections, the exact nonlinear strain-displacement relations for a beam element are described. The shear strains formulated by the present modeling method in this paper are zero, so it is reasonable to use geometrically nonlinear deformation fields to demonstrate and simplify a flexible beam undergoing large overall motions. Then, considering the coupling terms of deformation in two dimensions, finite element shape functions of a beam element and Lagrange’s equations are employed for deriving the coupling dynamical formulations. The complete expression of the stiffness matrix and all coupling terms are included in the formulations. A model consisting of a rotating planar flexible beam is presented. Then the frequency and dynamical response are studied, and the differences among the zero-order model, first-order coupling model and the new present model are discussed. Numerical examples demonstrate that a ‘stiffening beam’ can be obtained, when more coupling terms of deformation are added to the longitudinal and transverse deformation field. It is shown that the traditional zero-order and first-order coupling models may not provide an exact dynamic model in some cases.     

On optimal geometrically non-linear trusses

The paper is devoted to the investigation of regularities inherent to optimal geometrically non-linear trusses. The single static loading case is considered, a single structural material is used (except specially indicated cases) and buckling effects are neglected. The so-called small strains and large rotations case is investigated. Some regularities inherent to the kinematic and static variational principles for geometrically non-linear trusses are considered. Then the strain compatibility conditions resulting from the static variational principle are obtained and explored. It is shown that 1) the conditions are linear with respect to subcomponents of rod Green strains such as rotations and geometrically linear strains, 2) strains (in particular, rotations and geometrically linear strains) within rods which are not members of the so-called basic structure are fully determined by geometrically linear strains in rods of the basic structure.Extensions of some theorems (Maxwells theorem, Michells theorem, theorems on the stiffness properties of equally-stressed structures, etc.) known for geometrically linear structures are proved.Conditions assuring better or worse quality of equally-stressed geometrically non-linear truss as compared to geometrically linear ones are obtained.It is shown that in numerical optimization of geometrically non-linear trusses in the case of negligible rotations of compressed rods some updated analytical optimization algorithms (derived earlier for geometrically linear case) are monotonic. A simple numerical example confirming the features is presented.     

Hybrid strain based three node flat triangular shell elements—I. Nonlinear theory and incremental formulation

Aspects and theories of nonlinear analysis of structures, with special emphasis on structures that are discretized by the finite element method, are discussed. The updated Lagrangian formulation and the incremental Hellinger-Reissner variational principle are adopted. The independently assumed fields employed are the incremental displacements and incremental strains. Accordingly, the incremental second Piola-Kirchhoff stress and the incremental Washizu strain are selected as the incremental stress and strain measures. Various schemes for the transformation of the second Piola-Kirchhoff stress to Cauchy stress are included. Two versions of linear and nonlinear element stiffness and mass matrices are considered. These are the director and simplified versions. Variable thickness of the shell is considered so as to account for the ‘thinning effect’ due to large strain. Material nonlinearity studied in this paper is of elasto-plastic type with isotropic strain hardening. Cases in which small elastic but large plastic strain condition applies are considered and the J₂ flow theory of plasticity, in conjunction with Ilyushin's yield criterion, is employed. To simplify the derivation of (small displacement) stiffness matrix and to facilitate the derivation of explicit expressions for the element matrices, the non-layered approach has been applied.     

An implicit incrementally objective formulation for the solution of elastoplastic problems at finite strain by the F.E.M.

The mechanical formulation presented in this paper is based on an incremental updated Lagrange procedure using the principle of virtual work at the end of each load increment and an implicit incremental flow rule obtained by an approximate time integration of the objective rate constitutive equations. The approximate time integration is carried out along a particular path in the deformation and rotation space. This path ensures the incremental objectivity and minimizes the equivalent strain over the increment among all the possible paths, consequently avoiding an artificial increase of the plastic equivalent strain during the interpolation. The mechanical formulation presented leads to a fixed set of nonlinear equations, whose unknowns are the nodal displacements of the structure. A numerical algorithm based on a quasi Newton-Raphson method is then proposed to solve this system. The separation of the mechanical formulation from the resolution algorithm ensures the path independence. Numerical tests are carried out for a material obeying an isotropic with work-hardening von Mises criterion and associated flow rule. Single element tests show that this approach gives a very accurate solution even when the strain increment reaches twenty times the elastic strain up to yield. A structural test on a beam measures the influence of the incremental objectivity on the displacements, the equivalent plastic strain and the stresses.     

Static and dynamic finite deformations of cables using rate equations

A new incremental technique is presented which facilitates the numerical solution of static and dynamic finite deflection problems. By means of this technique, termed the Rate Equation Method, a set of linear differential equations in the displacement rates and stress rates are obtained from the nonlinear differential equations associated with finite deflections. These linear differential equations are solved numerically for the rate variables which in turn are integrated forward in time to yield the displacements and stresses at each time increment. By using second order finite differences to express the rates of accelerations in terms of velocities, static and dynamic problems may be solved by the identical procedure. To illustrate the details of the method, static and dynamic geometrically nonlinear deformations of cables under various loads are investigated and numerical results presented for several cases.     

Comparison between two theories of plasticity

Among numerous large strain elasto-plasticity theories, Green-Naghdi's theory and E. H. Lee's theory are distinguished and distinctive. In Green-Naghdi's theory, the Green-Lagrange strain tensor is decomposed into the elastic and the plastic parts. On the other hand, E. H. Lee started with a decomposition of the deformation gradient into a product of two parts: elastic and plastic. In the case of simple tension, the essential differences arc found between these two theories. In E. H. Lee's theory, the unloading curves are parallel on the plots of Cauchy stress vs natural strain. However, this parallel relation does not exist on the plots of Piola-Kirchhoff stress vs Green-Lagrange strain. In Green-Naghdi's theory, the results are reversed. The unloading curves are not parallel on the plots of Cauchy stress vs natural strain, but parallel on the plots of Piola-Kirchhoff stress vs Green-Lagrange strain. The significance of this finding is further discussed.