Corotational formulation for 3d solids. An analysis of geometricaly nonlinear foam deformation.

# Corotational formulation for 3d solids. An analysis of geometricaly nonlinear foam deformation.

Łukasz Kaczmarczyk Tomasz Koziara Chris J. Pearce Department of Civil Engineering, University of Glasgow, Rankine Building
###### Abstract

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.

## 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 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 [7]. 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

 u(~Xf)=UΓ(~Xf)q. (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 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:

 σ=Sv(Xf)v (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 (co-rotated) coordinates. In eq. (3), the stress approximation field is chosen in order to automatically satisfy the equilibrium condition (LABEL:eq1):

 LTSv=0. (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

 t=nσ=nSvv (5)

Furthermore, assuming small strains and elastic body, the displacement field within the domain is directly associated with this stress approximation expressed as:

 u=Uvv. (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

 x(X)=RX (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 nine-dimensional space of all matrices and a surface defined by the constraint . The surface normal is orthogonal to vectors or in other words

 ˙RRT+R˙RT=0. (8)

Let us define as

 ω=˙RRT. (9)

Since equation (8) is satisfied for all and we note that is anti-symmetric, .

The operator bears the name of the spatial angular velocity

 ˙x=ωx. (10)

For the ordinary differential equation (10) has a known solution

 x(t)=Exp(tω)R(0)X(0). (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

 δR=δΦR (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

 xr(X)=RX (13)

where is the position vector in the reference configuration. With that at hand we define displacement by

 ur(X)=xr(X)−X=(R−I)X (14)

Displacements are additively decomposed into the deformational and the rotational part, which yields

 ud(X)=u(X)+(I−R)X=u(X)−ur(X) (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

 δud(X)=δu−δRX=δu−Spin[δ→Φ]R(X−a)=δu−Spin[δ→Φ]xr=δu+Spin[xr]δ→Φ (16)

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

 Spin[x]=⎡⎢⎣0−x3x2x30−x1−x2x10⎤⎥⎦ (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 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 [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.

 ∫Ω⎛⎜⎝[∂ud∂X]−[∂ud∂X]T⎞⎟⎠dΩ=0 (18)

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

 →h=∫Γn×uddΓ=∫ΓSpin[n]uddΓ (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 pseudo-vector is used to formulate the best fit functional. Utilizing that , the functional takes form

 J(R)=12[→h]T[→h]=12[∫ΓSpin[n]uddΓ]T[∫ΓSpin[n]uddΓ]=12[∫ΓSpin[RN](u−ur)dΓ]T[∫ΓSpin[RN](u−ur)dΓ]=12[∫ΓSpin[N]RTxdΓ]T[∫ΓSpin[N]RTxdΓ] (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

 H=[∫Γ% Spin[N]RTxdΓ]T[∫ΓSpin[N]RTSpin[x]dΓ]=0 (22)

and solve the above sytem by means of the Newton-Raphson iterations

 ∂H∂ΦΔΦ=−H (23)
 Rk+1=Exp(ΔΦ)Rk. (24)

In order to approximate (23) we linearise

 Exp(ΔΦ)=I+sin∥Δ→Φ∥∥Δ→Φ∥ΔΦ+1−cos∥→ΔΦ∥∥Δ→Φ∥2ΔΦ2≈I+ΔΦ (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 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.

 δqd=Pδq=(I−SG)δq (27)

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

 δur(X)=UΓ(~Xf)Sf(Xfgc)δ→Ω (28)

where is a face geometrical center. The Spin-Liver matrix itself, for our particular implementation of the HTS element has simple form

 S(Xfgc)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0xr3−xr2−xr30xr1xr2−xr100∂xr3/∂~X1−∂xr2/∂~X1−∂xr3/∂~X10∂xr1/∂~X1∂xr2/∂~X1−∂xr1/∂~X100∂xr3/∂~X2−∂xr2/∂~X2−∂xr3/∂~X20∂xr1/∂~X2∂xr2/∂~X2−∂xr1/∂~X20000…⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (29)

where and . Definition of the Spin-Liver 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 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

 J(δ→Φ)=12[∫ΓSpin[N]RT(δu+Spin[xr]δ→Φ)dΓ]T[∫ΓSpin[N]RT(δu+% Spin[xr]δ→Φ)dΓ] (30)

Stationary point of this functional yields

 ∂J(δ→Φ)∂δ→Φ:δ→Φ=[∫ΓSpin[N]R% TSpin[xr]dΓ]T[∫ΓSpin[N]RT(δu+Spin[xr]δ→Φ)dΓ]=0. (31)

Equation of the stationary point is given by the system of three linear equations

 [∫ΓSpin[N]RTSpin[xr]dΓ]T[∫ΓSpin[N]RTSpin[xr]dΓ]δ→Φ=−[∫ΓSpin[N]R% TSpin[xr]dΓ]T[∫ΓSpin[N]RTδudΓ]. (32)

Substituting yields

 [∫ΓSpin[N]RTSpin[xr]dΓ]T[∫ΓSpin[N]RTSpin[xr]dΓ]δ→Φ=−[∫ΓSpin[N]R% TSpin[xr]dΓ]T[∫ΓSpin[N]RTUΓdΓ]δq (33)

As a result, a matrix representing the Spin-Filter

 δ→Φ=Gδq, (34)

is expressed in the form

 G=−{[∫ΓSpin[N]RTSpin[xr]dΓ]T[∫ΓSpin[N]RTSpin[xr]dΓ]}−1[∫ΓSpin[N]RT% Spin[xr]dΓ]T[∫ΓSpin[N]RTUΓdΓ]. (35)

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 [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 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.

 (36)

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

 Rε=∫Ωtr[w2ε]dΩ−∫Γ% u(w2n)T(uΓ−c−ur)dΓ−∫Γσ(w2n)T¯¯¯¯udΓ=0 (37)
 Rε=∫Ωtr[w2ε]dΩ−∫Γ% u(w2n)TudΓdΓ−∫Γσ(w2n)T¯¯¯¯udΓ=0 (38)
 Rε=∫Ωtr[R^w2RTε]dΩ−∫Γu(R^w2^n)TudΓdΓ−∫Γσ(R^w2^n)T¯¯¯¯udΓ=0 (39)
 (40)
 Rε=∫Ωtr[^w2^ε]dΩ−∫Γu(^w2^n)TRTudΓdΓ−∫Γσ(^w2^n)TRT¯¯¯¯udΓ=0 (41)

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.

### 6.3 Linearization

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

 (42)
 (43)
 (44)

In a similar manner to Eq. (36), residual Eq. (41) is linearized and yields

 Rε+∫Ωtr[^w2δ^ε]dΩ−∫Γ(^w2^n)TRTδudΓdΓ−∫Γ(^w2^n)TδRTudΓdΓ=0 (45)
 Rε+∫Ωtr[^w2δ^ε]dΩ−∫Γ(^w2^n)TRTδudΓdΓ+∫Γ(^w2^n)TRTδΦudΓdΓ=0 (46)
 (47)

A key term in above equations is the change of the angular vector . We note depends on the current configuration and is a linear function of the variation of displacement degrees of freedom, see Eq. (34) and Eq. (35).

### 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

 Rσ+∫ΓUTΓ^n^SvdΓδv−∫ΓUTΓSpin[R^σ^n]dΓGδq=0, (48)
 Rε+∫Γ^n^SvTUvdΓδv−∫Γ^n^SvTRTUΓdΓSGδq−∫Γ^n^SvTRTSpin[udΓ]dΓGδq=0. (49)

The above linearized equations can are expressed conveniently in the matrix form

 (50)

where

 F=∫Γ^n^SvT^UvdΓ, (51)
 A=∫Γ^n^SvTRTUΓdΓ, (52)
 Q=∫Γ^n^SvTRTSpin[udΓ]dΓ, (53)
 U=∫ΓUTΓSpin[R^σ^n]dΓ. (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

 {−AF−1[AT(I−SG)+QG]−UG}δq=−Rσ−AF−1Rε. (56)

This leads to a modified system of equations

 ^Fδq=−^Rσ. (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 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) [8] and A Mesh-Oriented 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 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

 κ=|w′′|(1+w′2)3/2 (58)

and the bending moment according to Bernoulli theory is given by

 M=EIκ. (59)

Strain energy is calculated from

 E=12∫l/2−l/2EIκ2dl (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

 u=0,v=ycos12Lκ−y,w=κ−1sin12Lκ−z (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 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 [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 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 [11], 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

 Rq=^Rσ−λTλ, (62)

 Rλ=GΔqTGΔq+Ψ2λ2−s2. (63)

is an increment of the rotation angle , which yields

 Rλ=Δ→ΦTΔ→Φ+Ψ2λ2−s2, (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 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.

## 8 Conclusions

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.

Findings from foam analysis can be used to model auxetic and non-auxetic polymeric foams [12], remodelling and growth of tubercular bone [13] and fracture of ceramic/quasi-brittle foams.

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.

## References

• [1]
• [2] J.A. Teixeira de Freitas, Formulation of elastostatic hybrid-Trefftz stress elements, Comput. Methods Appl. Mech. Engrg. 153 (1998), pp. 127â151.
• [3] 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.
• [4] B. Nour-Omid, 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 three-dimensional beam element, Comput. Methods Appl. Mech. Engrg. 81 (1990) 131â150.
• [6] M.A. Crisfield, G.F. Moita, A unified co-rotational for solids, shells and beams, Int. J. Solids Struc. 33 (1996) 2969â2992.
• [7] 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
• [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, MOAB-SD: integrated structured and unstructured mesh representation, Engineering with Computers, 20(3) 2004 286-293
• [10] Simpleware Ltd Web page
• [11] 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.
• [12] 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.
• [13] L. Kaczmarczyk and C.J. Pearce, Efficient Numerical Analysis of Bone Remodelling, Journal of the Mechanical Behavior of Biomedical, 2010.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters