Corotational formulation for 3d solids. An analysis of geometricaly nonlinear foam deformation.
Abstract
This paper presents theory for the Lagrange corotational (CR) formulation of finite elements in the geometrically nonlinear analysis of 3D structures. In this paper strains are assumed to be small while the magnitude of rotations from the reference configuration is not restricted. A new best fit rotator and consistent spin filter are derived. Lagrange CR formulation is applied with Hybrid Trefftz Stress elements, although presented methodology can be applied to arbitrary problem formulation and discretization technique, f.e. finite volume methods and lattice models, discreet element methods. Efficiency of CR formulation can be utilized in postbuckling stability analysis, damage and fracture mechanics, modelling of dynamic fragmentation of bodies made from quasibrittle materials, solid fluid interactions and analysis of poststressed structures, discreet body dynamics.
keywords:
Lagrange and corotated formulation, large rotations, foam, geometrically nonlinear, Trefftz elements, and Corresponding author.
1 Introduction
Nowadays the Lagrangian kinematics description for Finite Element Method of geometric nonlinear structures is utilized mainly in two solution techniques, i.e. Total Lagrangian Formulation and Updated Lagrangian Formulation. In this paper we investigate a less common solution technique for geometrically nonlinear problems, i.e. the Lagrange Corotational (CR) formulation.
Key features of the Lagrange CR formulation are efficiency and robustness, i.e. computation of tangent matrices is is lowpriced, the method is easy to parallelize and can be used with various discretization techniques. Large class of problems with finite rotations but small strains, e.g. postbuckling stability analysis, damage and fracture mechanics, modelling of dynamic fragmentation of bodies made from quasibrittle materials, solid fluid interaction, analysis of poststressed structures or discreet body dynamics can be solved efficiently with the use of Corotational formulation.
As noted above, the Lagrange CR formulation can be applied together with discretization techniques not suitable for Total Lagrangian Formulation or Updated Lagrangian Formulation. In this paper the CR formulation is applied with Hybrid Trefftz Stress (HTS) finite elements. Application of the CR formulation together with HTS elements results in a powerful synergy, which is exemplified on a case study problem presented in the paper. Although the presented methodology is applied to problems of solid mechanics, the presentede method can be applied to image recovery, where extraction of rigid body motion from image is an issue.
This paper follows earlier contributions [4, 5, 6, 7]. Full literature review can be fond in [7]. Despite the large body of literature on corotational formulations, formulation of a best fit rotation functional remains an open reasearch question. In this paper we introduce a new general bestfit rotator functional. Additionally a new SpinFilter operator, mapping increment of degrees of freedom into increment of rotation vector, consistent with new bestfit rotator functional is determined.
In the first section approximation functions for the HTS elements are briefly described. Next, kinematics of rotating deformable motion is described and the bestfirt rotator functional is formulated. Section 4 is supplemented with linearized equations, used to determine the rotation operator. In Section 5 a SpinLiver linear operator transforming an increment of rotation angle into an increment of displacement degrees of freedom, and SpinFilter are formulated. With projection matrices at hand, a linearization equations of HTS elements are shown in Section 4. In Section 7 two numerical examples are presented. The paper is concluded in Section 8.
2 Approximation functions
Lagrange CR formulation can be applied to a large class of discretization methods. In this paper we apply it to a particular instance of the Hybrid Trefftz Stress elements. In authors view combination of the HTS element with the CR formulation is a good example of the robustness of the corotational approach: we apply the CR formulation to a discontinuous displacement approximation, exclusively determined only on element boundary. Moreover, unlike for the classical displacement Finite Element Method, the employed displacement approximation does represent a partition of unity.
2.1 Approximation of displacements on element faces
In the spirit of hybrid FE formulations, displacement field is approximated on element boundary, independently on each element face. In the presented approach we used simple monomials given in the base coordinate system of element faces, see Eq. (2). Displacements of face at point (expressed in face coordinates) read
(1) 
Matrix of displacement approximation functions has simple form
(2) 
where is a coordinate vector in the reference system of element face.
Such approximation technique allows for an easy extension of the approximation space by arbitrary functions, which are not necessarily a partition of unity. Moreover, since degrees of freedom are associated with elements faces, for threedimensional tetrahedrons (utilized here), an element face can have no more than six neighbours (adjacent through elements on both sides). As a result, this formulation is suitable for parallel computing where intercommunication between processors affects efficiency.
2.2 Approximation of stresses in element volume
Detailed description of HTS elements can be found in [2, 3]. In this paper we limit ourselves to basic equations. In this section approximation of stresses is briefly described. The Cauchy stress field within an element is approximated directly as:
(3) 
where is a matrix of field approximation functions and is the unknown vector of generalised stress degrees of freedom. We note that is defined on element face in Lagrange (corotated) coordinates. In eq. (3), the stress approximation field is chosen in order to automatically satisfy the equilibrium condition (LABEL:eq1):
(4) 
where is a differential operator related to the equilibrium equations. According to the Cauchy equilibrium condition, the corresponding stress induced traction field on the element boundary is
(5) 
Furthermore, assuming small strains and elastic body, the displacement field within the domain is directly associated with this stress approximation expressed as:
(6) 
where is a displacement approximation matrix, such that , where is a compliance matrix.
The stress field and the associated displacement field are derived from appropriate Trefftz functions. For example, the stress and displacement field could be given as polynomial functions derived from the Airy’s stress function. Alternatively, the Trefftz functions can be extended to functions which can take into account singularities or are suitable for a particular boundary value problem. Moreover, the approximation of stress is not a priori dependent on a particular constitutive law.
3 Kinematics of CR frame
For completeness some basic equations and theory of rotations are presented.
Consider a body under rotational motion
(7) 
where is x is the rotation operator (also called rotor), is a selected material point and is the spatial image of . Since is the rotation operator, distances between points and orientation are preserved, i.e. and .
In mathematical terms the set of all orthogonal matrices with positive determinant forms the special orthogonal group , which is a manifold. Since rotations are points of the manifold , at each point in a tangent Euclidean space collects rotation velocities .
Consider now a ninedimensional space of all matrices and a surface defined by the constraint . The surface normal is orthogonal to vectors or in other words
(8) 
The operator bears the name of the spatial angular velocity
(10) 
For the ordinary differential equation (10) has a known solution
(11) 
where time is expressed explicitly by and Exp stands for the matrix exponential. Although in general is not constant in time, the above formula serves as a useful device for implementing incremental updates of .
3.1 Variations of rotation operator
With the brief theory at hand, we introduce variations in order to construct tangent operators for the HTS formulation. From (9) we can see that and hence the variation of can be expressed in the spatial form
(12) 
where .
3.2 Variation of displacements
For simplicity and without loss of generality we analyse motion without the rigid body translations. New position of a particle is given by
(13) 
where is the position vector in the reference configuration. With that at hand we define displacement by
(14) 
Displacements are additively decomposed into the deformational and the rotational part, which yields
(15) 
where and are the deformational and the rotational parts of displacements, respectively.
In order to derive tangent operators, a useful variation of deformational displacements (15) is given by
(16) 
where variation of the axial rotation vector is used. We note, that is a pseudovector, i.e. a 3dimensional representation of a antisymmetric matrix. The operator takes form
(17) 
for an arbitrary vector .
4 Best fit rotor
Computation of the rotation operator for a given displacement field is at the hart of the corotational formulation. In our view, an ideal corotational formulation should address the following objectives:

Versatility. It must work for static and dynamic problems, with arbitrary discretization of displacements.

Rigid bodies. For rigid body motion, it should give the same results as a bodyattached frame. This simplifies coupling of FEMCR methods with multibody dynamics codes [7].

Finite stretch and pure sheer. If a element is subjected to a uniform finite stretch or pure sheer, no spurious rotation should appear.
In order to find the rotation operator , a new functional is proposed here. Stationary points of this functional determine rotation, which complies with the above objectives.
We note, that for a rotation free deformation, the antisymmetric part of the displacement gradient disappears. With that observation at hand, assuming that the rotation is constant in the element domain, the antisymmetric part of the volume averaged deformational displacement is used to construct the best fit rotator.
(18) 
Noting that the antisymmetric part of the displacement gradient has three nonzero linearly independent components, the antisymmetric part of the displacement gradient can be uniquely expressed by a pseudovector . Utilizing that and applying integration by parts, a pseudovector can be expressed by the boundary integral
(19) 
where is the spatial normal vector field on the element boundary in the current or reference configuration (since the motion is roughly rigid we can assume the Jacobian to be nearly one). A length of the pseudovector is used to formulate the best fit functional. Utilizing that , the functional takes form
(20) 
where is the referential normal vector field on the element boundary in the current or reference configuration. We note, that the functional expressed in the above form can be applied to problems where displacements are only known on elements boundaries (and can be discontinuous).
For given displacements , a new position vector in the current configuration is easy to determine. For a fixed , the condition of stationarity of represents the set of nonlinear equations for finding
(21) 
where (12) has been utilized. We then requre that
(22) 
and solve the above sytem by means of the NewtonRaphson iterations
(23) 
(24) 
In order to approximate (23) we linearise
(25) 
and noting that we represent (23) by up to first order terms of . That is
(26) 
where is the Cartesian base vector, and represents the th column of .
In this section a new bestfit functional, based on boundary integrals is presented. Stationary point of this functional lead to a system of nonlinear algebraic equations, form which current rotation operator can be calculated by means of Newton iterations. It should be noted, that the above procedure allows to determine the rotation quite independently of the motion history. However, in order to avoid troubles with multiple minima and correctly trace the history of rotation it is appropriate to initialize the iterative procedure with the most recent value of rotation.
5 Projection matrices
In order to construct tangent operators for the consistent Newton scheme projection matrices filtering a variation of the deformational displacements from the variation of the total displacements have to be derived.
(27) 
Projection matrix , SpinLiver and SpinFilter depend on the approximation method. In this paper we construct these matrices for the particular approximation of HTS finite elements. However, the methodology is general and can be applied to a large class approximation methods.
5.0.1 Spin Liver
A SpinLiver matrix represents a linear operator transforming a variation of the rotation axis into a variation of linear displacements. SpinLiver depends only on constant and linear approximation functions. Higher order terms in displacements approximation basis are neglected in the formulation of SpinLiver. It can be noted, that for any displacement function given by higherorder terms the motion is not stress free, hence higherorder therms cannot be a part of the rigid body motion.
Exploiting features of the HTS displacement approximation, see Eq. (2), the SpinLiver can be defined for each face independently, which yields
(28) 
where is a face geometrical center. The SpinLiver matrix itself, for our particular implementation of the HTS element has simple form
(29) 
where and . Definition of the SpinLiver for other types of elements can be found in [7].
5.0.2 Spin Filter
Matrix represents a linear operator which transforms a variation of linear displacements into a variation of rotation axis. In other words SpinFilter accounts for a change of magnitude and direction of the rotation axis due to the change of linear displacements.
In order to get quadratic convergence, SpinFilter matrix has to be consistent with functional , defined in Section 4. For a given current configuration, i.e. rotation , and linearized displacements, functional takes form
(30) 
Stationary point of this functional yields
(31) 
Equation of the stationary point is given by the system of three linear equations
(32) 
Substituting yields
(33) 
As a result, a matrix representing the SpinFilter
(34) 
is expressed in the form
(35) 
The novelty of this approach in comparison to those found in other papers related to the CR formulation, is that the SpinFilter is derived form best fit functional and can be therefore applied to a large class of approximation methods.
6 HTS formulation
Although the corotational formulation can be derived independently form a problem formulation by consistent differentiation of projection matrices [7], in this paper we will formulate tangent operators and residuals directly from weak forms. In authors opining this allows for better understanding of the problem, without loss of generality.
The corotational Hybrid Trefftz Stress element is built using Lagrange coordinates, i.e. integrals are over reference volumes and faces. Stresses are expressed in the corotated (current) configuration and displacements are expressed in the reference coordinate system. At each iteration corotated (Cauchy) stress is pulledback to the reference coordinate system (or conversly, the displacement is pushedforward to the corotated reference system).
6.1 Integral form
Static boundary conditions are enforced in strong integral form, i.e. we weight static boundary condition by function and integrate over element boundary.
(36) 
Function is given in the reference coordinate system and stresses are pulledback to the reference coordinate system. This results in geometrical nonlinearities, i.e. the integral becomes a function of the unknown rotor .
6.2 Weak form
The compatibly conditions within a body domain are imposed in a weak sense, where by the differentiation by parts Dirichlet boundary conditions are included
(37) 
(38) 
(39) 
(40) 
(41) 
A weighting function is defined in the corotated coordinate system and is transformed to the reference coordinate system. We note, that the term under the volume integral is a scalar, so it not depend on the coordinate system. As a result it is constant for all steps analysis and is computed only once.
6.3 Linearization
6.4 Matrix form
The linearized equations are reformulated, with stresses and displacements substituted by suitable products approximation functions and increments of vectors storing stress and displacement degrees of freedom. SpinFilter and SpinLiver are utilized in order to express linearized equations exclusively in terms unknown vectors, which yields
(48) 
(49) 
The above linearized equations can are expressed conveniently in the matrix form
(50) 
where
(51) 
(52) 
(53) 
(54) 
6.5 Partial Solution
The size of the problem to be solved can be reduced by static condensation, taking advantage of the element level approximations of stress fields. On element level, solution for stress degrees of freedom is expressed by
(55) 
Substituting in Eq. (48) yields
(56) 
This leads to a modified system of equations
(57) 
In the above, the stress degrees of freedom are eliminated conveniently from the global system of equations, leaving only the displacement degrees of freedom to be determined. This does not affect the final result in any way and simply represents a solution technique that reduces the number of equations that have to be solved simultaneously. Following the solution of the displacement degrees of freedom, the stress degrees of freedom can be recovered on an element by element basis, a process that can be easily parallelized.
It is noted, that independent displacement approximation on elements faces, as presented here, lead to larger number of unknowns in comparison to classical FEM. Nevertheless, the matrix sparsity, which has a significant impact on the efficiency of direct and iterative solvers, is minimised in the presented approach due to the association of the displacement degrees of freedom with faces rather than vertices and due to the fact that any given face can only have two neighbouring elements in either 2D or 3D. For the same reason the communication between parallel processors is significantly reduced. As shown in [br], this approach results in very good scalability.
7 Numerical examples
To illustrate features of the CRHTS approach, two numerical examples are presented. The first one is an academic problem investigating model’s accuracy. The second example further illustrates model’s accuracy for large scale 3D representative model of complex micro structure.
Examples are analysed with the use Portable Extensible Toolkit for Scientific Computation (PETSc) [8] and A MeshOriented datABase (MOAB) [9]. Finite elements and CR formulation were implemented in C and C++.
7.1 Pure bending in 3D
In the first numerical example the numerical results are compared with an analytical solution. A beam (0.2x0.2x5.0) with square crosssection is stressed under pure bending.
In order to get analytical solution, a special case is analysed, i.e. and . A plane curvature is given explicitly as
(58) 
and the bending moment according to Bernoulli theory is given by
(59) 
Strain energy is calculated from
(60) 
With the analytical solution at hand, a 3D numerical model is analysed. A discretized beam consisting 1160 tetrahedrals as shown in Fig. 5. A HTS element with quadratic displacements discretization and 3rd order polynomial approximation of stress field within element is used. At both ends the kinematic boundary conditions are applied in order to force pure bending
(61) 
where the length of the beam . Despite the coarse mesh, numerical solution is in excellent agreement with the analytical one, cf. Fig. 6. For beam forms a circle, which is in agreement with the analytical solution. At each time step quadratic convergence is observed.
This test has been reproduced for several bar placements in 3D, in order to verify the implementation and the formulation of bestfit rotator. We note that the Cauchy stress field distribution is in an exact much which the analytical solutions for all cases.
At each increment an updated rotor is computed for each element independently. At all times the procedure presented in section 4 required few iterations to achieve the machine precision. Moreover, profiling the runtimes showed that the time related to the rotor computations was negligible.
7.2 Analysis of foam
In order to present one of the many potential applications of the Lagrange CR formulation we investigate deformation of a hypothetical foam. The mesh was obtained from a highresolution 3D scan, included as a case study in the Simpleware software [10]. The mesh comprised 95956 tetrahedrons and 40725 nodes. The constitutive parameters of the model were chosen arbitrarily, i.e. and . At the bottom plate a fixed boundary conditions were applied. At the top plate tractions in direction were applied. On the front, back, left and right side the displacement in the normal direction to the side wall was blocked.
Displacement and stresses filed were approximated using HybridTrefftz elements. Two cases were considered with linear and quadratic approximation of displacements. For the first case, i.e. linear functions approximating displacements, polynomial Trefftz functions of order two were used to approximate stresses. For the second case, polynomial Trefftz functions of order three were used to approximate the stress field. This led to the model with million DoFs and million DoFs, for the case one and case two, respectively. We note, that we not observe shear locking phenomenon for the HTS elements [11], despite tetrahedrons being used to approximate the geometry.
Large rotations led to local buckling in the microstructure which had an impact on the convergence of the Newton method. In order to improve efficiency, a set of nonlinear equations was supplemented with the arclength control. We introduced additional unknown , such that
(62) 
and additional control equation
(63) 
is an increment of the rotation angle , which yields
(64) 
where elements of matrix are determined at the beginning of a load increment. and are scaling parameter and load factor respectively. is radius of a hypersphere. We note that the change of the rotation angle is a source of nonlinearities and a direct control of the rotation angle positively affects the stability of the Newton method. Since controlling geometrical instabilities in a structure like foam is an open problem, the above solution has been adopted only as a demonstration, while further improvements could well be introduced.
Analysis of foam deformation under compression and tension showed expected foam response, i.e. overall response softening and stiffening for compression and tension, respectively. Deformation, rotation angle and stress are shown in Fig. 10 and Fig. 11. We can note that for compression we observed interpenetration and the model should be extended to take into account selfcontact. For compression, a transition to bending dominated behavior was observed. For tension, a straightening of foam’s branches was observed. In both cases local buckling was a source of numerical difficulties.
Comparison between linear and quadratic approximation of displacements, i.e. case 1 and case 2, showed a close much between the numerical responses. Good quasiquadratic convergence for case 1 and case 2 was achieved until the final point of analysis. However poor mesh quality and local instabilities were the source of problems with the convergence for extreme (nonphysical) loads. Those problems need further investigation and in authors opinion can be overcome by a better local arclength control and line searches.
In figures showing deformation, displacements are plotted using Trefftz displacement function , see Eq. 6, supplemented with rigid body rotation . Trefftz displacement field is adjoint to stresses calculated to fulfill equilibrium, hence displacement vector in Eq. 6 do not carry any information about translations and rotations. By visual inspection of the figures in this example and in Fig. 6 we can verify an excellent consistency of the calculated rotor with the overall deformation.
8 Conclusions
This paper presents a geometrically nonlinear Lagrange corotational formulation applied to Hybrid Stress Trefftz Elements. All terms in the internal force and tangent matrices are accounted for.
A new bestfit rotator for continuum elements is developed. The bestfit rotator functional is constructed by utilizing a boundary integral over a deformable element/body. Stationary point of the functional lead to the set of three nonlinear equations, solved by means of the Newton method. Formulation of bestfit rotator, load vector and tangent matrices was verified by two examples: pure bending test and large deformation of foam.
Findings from foam analysis can be used to model auxetic and nonauxetic polymeric foams [12], remodelling and growth of tubercular bone [13] and fracture of ceramic/quasibrittle foams.
In authors opinion exciting research opportunities arise form this contribution. The CRHTS method can be applied to many problems involving quasibrittle materials and large deformations. We note, that for fracturing materials very often large rotations are observed for fractured parts. Moreover the current approach can be applied to many dynamic problems involving large rotations and small deformations. For models involving in addition large strains, it is possible to use the CR formulation together with Total Lagrangian or Updated Lagrangian formulations, where the CR frame could be attached to rotating parts of the model or individual elements.
Another interesting research opportunity is application of the Lagrange CR formulation to very large problems, solved with Krylov iterative solver without explicit formation of the tangent matrix. Matrix free methodology is attractive for memory demanding problems. We note, that in case of the presented methodology, all integrals can be computed only at the beginning of an analysis, and the tangent matrix can be calculated with the use of projection matrices only. As a result, the matrix free methodology can be used efficiently with nonlinear problems and higherorder elements without the need for integration at each iterative solver step.
References
 [1]
 [2] J.A. Teixeira de Freitas, Formulation of elastostatic hybridTrefftz stress elements, Comput. Methods Appl. Mech. Engrg. 153 (1998), pp. 127â151.
 [3] Lukasz Kaczmarczyk, Chris J. Pearce, A corotational hybridTrefftz stress formulation for modelling cohesive cracks,Comput. Methods Appl. Mech. Engrg, Volume 198, Issues 1516, 15 March 2009, Pages 12981310.
 [4] B. NourOmid, C.C. Rankin, Finite rotation analysis and consistent linearization using projectors, Comput. Methods Appl. Mech. Engrg. 93 (1991) 353â384.
 [5] M.A. Crisfield, A consistent corotational formulation for nonlinear threedimensional beam element, Comput. Methods Appl. Mech. Engrg. 81 (1990) 131â150.
 [6] M.A. Crisfield, G.F. Moita, A unified corotational for solids, shells and beams, Int. J. Solids Struc. 33 (1996) 2969â2992.
 [7] C.A. Felippa, B. Haugen, A unified formulation of smallstrain corotational finite elements: I. Theory, Comput. Methods Appl. Mech. Engrg. 194 (2005) 2285â2335
 [8] Satish Balay and Kris Buschelman and William D. Gropp and Dinesh Kaushik and Matthew G. Knepley and Lois Curfman McInnes and Barry F. Smith and Hong Zhang, PETSc Web page, http://www.mcs.anl.gov/petsc, 2009.
 [9] Timothy J. Tautges, MOABSD: integrated structured and unstructured mesh representation, Engineering with Computers, 20(3) 2004 286293
 [10] Simpleware Ltd Web page, http://www.simpleware.com
 [11] J Jirousek, A Wroblewski, QH Qin, XQ H, A family of quadrilateral hybridTrefftz pelements for thick plate analysis, Comput. Methods Appl. Mech. Engrg. 127 (1995) 315344.
 [12] Samuel A. McDonald, Ghislain DedreuilMonet, Yong Tao Yao, Andrew Alderson, and Philip J. Withers, In situ 3D Xray microtomography study basic solid state physics comparing auxetic and nonauxetic polymeric foams under tension, Phys. Status Solidi B, 1â7 (2010) / DOI 10.1002/pssb.201083975.
 [13] L. Kaczmarczyk and C.J. Pearce, Efficient Numerical Analysis of Bone Remodelling, Journal of the Mechanical Behavior of Biomedical, 2010.