Corotational formulation for 3d solids. An analysis of geometricaly nonlinear foam deformation.
This paper presents theory for the Lagrange co-rotational (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 post-buckling stability analysis, damage and fracture mechanics, modelling of dynamic fragmentation of bodies made from quasi-brittle materials, solid fluid interactions and analysis of post-stressed structures, discreet body dynamics.
keywords:Lagrange and co-rotated formulation, large rotations, foam, geometrically nonlinear, Trefftz elements
, and Corresponding author.
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 Co-rotational (CR) formulation.
Key features of the Lagrange CR formulation are efficiency and robustness, i.e. computation of tangent matrices is is low-priced, 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. post-buckling stability analysis, damage and fracture mechanics, modelling of dynamic fragmentation of bodies made from quasi-brittle materials, solid fluid interaction, analysis of post-stressed structures or discreet body dynamics can be solved efficiently with the use of Co-rotational 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 . Despite the large body of literature on co-rotational formulations, formulation of a best fit rotation functional remains an open reasearch question. In this paper we introduce a new general best-fit rotator functional. Additionally a new Spin-Filter operator, mapping increment of degrees of freedom into increment of rotation vector, consistent with new best-fit 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 best-firt rotator functional is formulated. Section 4 is supplemented with linearized equations, used to determine the rotation operator. In Section 5 a Spin-Liver linear operator transforming an increment of rotation angle into an increment of displacement degrees of freedom, and Spin-Filter 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 co-rotational 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
Matrix of displacement approximation functions has simple form
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 three-dimensional 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:
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 (co-rotated) coordinates. In eq. (3), the stress approximation field is chosen in order to automatically satisfy the equilibrium condition (LABEL:eq1):
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
Furthermore, assuming small strains and elastic body, the displacement field within the domain is directly associated with this stress approximation expressed as:
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
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 nine-dimensional space of all matrices and a surface defined by the constraint . The surface normal is orthogonal to vectors or in other words
Let us define as
Since equation (8) is satisfied for all and we note that is anti-symmetric, .
The operator bears the name of the spatial angular velocity
For the ordinary differential equation (10) has a known solution
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
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
where is the position vector in the reference configuration. With that at hand we define displacement by
Displacements are additively decomposed into the deformational and the rotational part, which yields
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
where variation of the axial rotation vector is used. We note, that is a pseudo-vector, i.e. a 3-dimensional representation of a antisymmetric matrix. The operator takes form
for an arbitrary vector .
4 Best fit rotor
Computation of the rotation operator for a given displacement field is at the hart of the co-rotational formulation. In our view, an ideal co-rotational 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 body-attached frame. This simplifies coupling of FEM-CR methods with multibody dynamics codes .
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.
Noting that the antisymmetric part of the displacement gradient has three non-zero linearly independent components, the antisymmetric part of the displacement gradient can be uniquely expressed by a pseudo-vector . Utilizing that and applying integration by parts, a pseudo-vector can be expressed by the boundary integral
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 pseudo-vector is used to formulate the best fit functional. Utilizing that , the functional takes form
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
where (12) has been utilized. We then requre that
and solve the above sytem by means of the Newton-Raphson iterations
In order to approximate (23) we linearise
and noting that we represent (23) by up to first order terms of . That is
where is the Cartesian base vector, and represents the th column of .
In this section a new best-fit functional, based on boundary integrals is presented. Stationary point of this functional lead to a system of non-linear 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.
Projection matrix , Spin-Liver and Spin-Filter 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 Spin-Liver matrix represents a linear operator transforming a variation of the rotation axis into a variation of linear displacements. Spin-Liver depends only on constant and linear approximation functions. Higher order terms in displacements approximation basis are neglected in the formulation of Spin-Liver. It can be noted, that for any displacement function given by higher-order terms the motion is not stress free, hence higher-order therms cannot be a part of the rigid body motion.
Exploiting features of the HTS displacement approximation, see Eq. (2), the Spin-Liver can be defined for each face independently, which yields
where is a face geometrical center. The Spin-Liver matrix itself, for our particular implementation of the HTS element has simple form
where and . Definition of the Spin-Liver for other types of elements can be found in .
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 Spin-Filter 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, Spin-Filter 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
Stationary point of this functional yields
Equation of the stationary point is given by the system of three linear equations
As a result, a matrix representing the Spin-Filter
is expressed in the form
The novelty of this approach in comparison to those found in other papers related to the CR formulation, is that the Spin-Filter is derived form best fit functional and can be therefore applied to a large class of approximation methods.
6 HTS formulation
Although the co-rotational formulation can be derived independently form a problem formulation by consistent differentiation of projection matrices , 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 co-rotational Hybrid Trefftz Stress element is built using Lagrange coordinates, i.e. integrals are over reference volumes and faces. Stresses are expressed in the co-rotated (current) configuration and displacements are expressed in the reference coordinate system. At each iteration co-rotated (Cauchy) stress is pulled-back to the reference coordinate system (or conversly, the displacement is pushed-forward to the co-rotated 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.
Function is given in the reference coordinate system and stresses are pulled-back 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
A weighting function is defined in the co-rotated 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.
The system of nonlinear algebraic equations is linearized and solved using Newton method. The rotation operator is a function of current displacements and it is computed at each iteration using the best-fit functional (20). At a current cofiguratiom, linearized Eq. (36) takes from
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. Spin-Filter and Spin-Liver are utilized in order to express linearized equations exclusively in terms unknown vectors, which yields
The above linearized equations can are expressed conveniently in the matrix form
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
Substituting in Eq. (48) yields
This leads to a modified system of equations
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 CR-HTS 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)  and A Mesh-Oriented datABase (MOAB) . 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 cross-section 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
and the bending moment according to Bernoulli theory is given by
Strain energy is calculated from
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
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 best-fit 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 high-resolution 3D scan, included as a case study in the Simpleware software . 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 Hybrid-Trefftz 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 , despite tetrahedrons being used to approximate the geometry.
Large rotations led to local buckling in the micro-structure 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 arc-length control. We introduced additional unknown , such that
and additional control equation
is an increment of the rotation angle , which yields
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 hyper-sphere. 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 self-contact. 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 quasi-quadratic 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 (non-physical) loads. Those problems need further investigation and in authors opinion can be overcome by a better local arc-length 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.
This paper presents a geometrically nonlinear Lagrange co-rotational formulation applied to Hybrid Stress Trefftz Elements. All terms in the internal force and tangent matrices are accounted for.
A new best-fit rotator for continuum elements is developed. The best-fit 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 best-fit rotator, load vector and tangent matrices was verified by two examples: pure bending test and large deformation of foam.
In authors opinion exciting research opportunities arise form this contribution. The CR-HTS method can be applied to many problems involving quasi-brittle 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 higher-order elements without the need for integration at each iterative solver step.
-  J.A. Teixeira de Freitas, Formulation of elastostatic hybrid-Trefftz stress elements, Comput. Methods Appl. Mech. Engrg. 153 (1998), pp. 127â151.
-  Lukasz Kaczmarczyk, Chris J. Pearce, A corotational hybrid-Trefftz stress formulation for modelling cohesive cracks,Comput. Methods Appl. Mech. Engrg, Volume 198, Issues 15-16, 15 March 2009, Pages 1298-1310.
-  B. Nour-Omid, C.C. Rankin, Finite rotation analysis and consistent linearization using projectors, Comput. Methods Appl. Mech. Engrg. 93 (1991) 353â384.
-  M.A. Crisfield, A consistent corotational formulation for nonlinear three-dimensional beam element, Comput. Methods Appl. Mech. Engrg. 81 (1990) 131â150.
-  M.A. Crisfield, G.F. Moita, A unified co-rotational for solids, shells and beams, Int. J. Solids Struc. 33 (1996) 2969â2992.
-  C.A. Felippa, B. Haugen, A unified formulation of small-strain corotational finite elements: I. Theory, Comput. Methods Appl. Mech. Engrg. 194 (2005) 2285â2335
-  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.
-  Timothy J. Tautges, MOAB-SD: integrated structured and unstructured mesh representation, Engineering with Computers, 20(3) 2004 286-293
-  Simpleware Ltd Web page, http://www.simpleware.com
-  J Jirousek, A Wroblewski, QH Qin, XQ H, A family of quadrilateral hybrid-Trefftz p-elements for thick plate analysis, Comput. Methods Appl. Mech. Engrg. 127 (1995) 315-344.
-  Samuel A. McDonald, Ghislain Dedreuil-Monet, Yong Tao Yao, Andrew Alderson, and Philip J. Withers, In situ 3D X-ray microtomography study basic solid state physics comparing auxetic and non-auxetic polymeric foams under tension, Phys. Status Solidi B, 1â7 (2010) / DOI 10.1002/pssb.201083975.
-  L. Kaczmarczyk and C.J. Pearce, Efficient Numerical Analysis of Bone Remodelling, Journal of the Mechanical Behavior of Biomedical, 2010.