Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

55, 3, pp. 977-990, Warsaw 2017
DOI: 10.15632/jtam-pl.55.3.977

VOLUMETRIC LOCKING SUPPRESSION METHOD FOR NEARLY

INCOMPRESSIBLE NONLINEAR ELASTIC MULTI-LAYER BEAMS

USING ANCF ELEMENTS

Grzegorz Orzechowski, Janusz Frączek

Warsaw University of Technology, Institute of Aeronautics and Applied Mechanics, Warsaw, Poland

e-mail: gorzech@meil.pw.edu.pl

The analysis and solution of many modern flexible multibody dynamic problems require
formulations that are able to effectively model bodies with nonlinear materials undergoing
large displacements and deformations. The absolute nodal coordinate formulation (ANCF)
in connectionwith a continuum-based approach is one way to deal with these systems. The
main objective of this work is to extend an existent approach for the modelling of slen-
der structures within the ANCF frameworkwith nonlinear, nearly incompressiblematerials
using the volumetric energy penalty technique. Themain part of the study is devoted to the
evaluation of multi-layer beam models and simplifications in the locking suppression me-
thod based onF-bar projection. The results present significantly better agreementwith the
reference solution formulti-layer structures built with the standardANCFbeam element as
compared with the earlier implementation.

Keywords: multibody dynamics, ANCF, incompressibility, locking phenomena, multi-layer
beams

1. Introduction

Thedynamic analysis of bodies that undergo large deformations and are builtwith complex and
nonlinearmaterials is a vital part of themodern computer-aided design andmodelling techniqu-
es. Therefore, such features should be included in a reliable manner in the advancedmultibody
system (MBS) simulation software. In flexible multibody dynamics, the most frequently used
method is the floating frame of reference formulation (Shabana, 1997b) that is usually limited to
linear-elastic deformations. The geometrical andmaterial nonlinearities can be included within
the finite element analysis (FEA) (Bathe, 1996), however, the FEA is not perfectly compatible
with theMBS (Wasfy and Noor, 2003).
The absolute nodal coordinate formulation (ANCF) proposed by Shabana (1997a) can be

efficiently usedwithin theflexiblemultibodydynamics.Theuniquecharacteristics of thismethod
allow straightforward modelling of beam and plate elements using nonlinear material models.
ANCF employs the slope coordinates rather than rotations to describe local orientation, which
enables, among other things, representation of complicated shapes using just a few elements.
Flexible ANCF bodies can exactly represent rigid body modes, including large rotations, and
model large body deformations. Additionally, the ANCF beam elements may employ general
constitutive formulations (in addition to the classical beam theories) for a variety of nonlinear
material models, including incompressible ones. All these features cause that the ANCF is well
suited for thedynamicanalysis of highlyflexiblebeamstructuresusingnonlinearmaterialmodels
within theMBS framework (Shabana, 2008).
Incompressible rubber-likematerials areused inmanyengineeringand industrial applications

like defence, automotive, safety andothers.Consequently, reliable and effective application of the
incompressible nonlinear materials in many biomechanical and engineeringmodels is one of the



978 G. Orzechowski, J. Frączek

key goals. However, commercial FEApackages restrict theworkwith incompressiblematerials to
shell and solid elements (ANSYS® Academic Research, 2010), also when slender structures are
considered that are otherwise modelled with beam elements. One can overcome this limitation
in the mentioned kind of common applications by applying fully parameterized ANCF beam
elements (Shabana, 2008; Orzechowski and Frączek, 2015).

Most investigations devoted to the ANCF framework assume a linear-elastic material mo-
del. Compressible and incompressible hyperelastic and isotropic material models were firstly
used within the ANCF by Maqueda and Shabana (2007). Furthermore, Maqueda et al. (2010)
presented the application of the ANCF beams with incompressible materials to rubber chains
systems. Moreover, the validation of the ANCF model based on the experiment that captured
motion of the rubber-like beam was presented in (Jung et al., 2011). Nonetheless, none of the
aboveworks have stressed the importance of using the locking alleviation techniques. It is worth
to point out thatmany issues that are actively researched in the ANCF field have already been
studied for nonlinear finite elements.

The main objective of this paper is to recall and extend the volumetric locking elimina-
tion techniques for nonlinear, hyperelastic, nearly incompressiblematerial models applied to the
ANCF beams. The previous paper by Orzechowski and Frączek (2015) showed an importance
of volumetric locking elimination techniques in typical applications, however, only higher-order
elements, like those presented in Orzechowski and Shabana (2016), which were also numerically
expensive due to the high number of coordinates and integration points, provided a reasonable
results. Therefore, exemplary techniques that may reduce the computational cost are introdu-
ced and validated with several numerical examples. Strictly speaking, the implementation of
multi-layer beammodels with appropriate continuity between layers make it possible to use the
lower-order ANCF beam element, while the use of lower-order projection basis with the F-bar
projection technique simplifies this locking suppression formulation. Therefore, in the current
study, the standard three-dimensional fully parameterized element is used (Shabana and Yako-
ub, 2001; Yakoub and Shabana, 2001) together with two simple incompressiblematerial models:
one-parameter Neo-Hookean and two-parameter Mooney-Rivlin (Shabana, 2008). Incompressi-
bility of thematerials is ensured by the penaltymethod, which is chosen due to its simple form
and common use, and two methods of the locking suppression are applied. Due to introduced
multi-layer structures, the results of numerical tests are in significantly better agreement with
the reference for models built with the standard ANCF beam element as compared with the
earlier implementation.

2. Kinematics and dynamics of deformable bodies

The nodal coordinates of the ANCF elements are prescribed with respect to the global refe-
rence frame and they include translational and slope coordinates. Consequently, no rotational
coordinates are used to identify the element orientation. Thus, the independent rotation field
interpolation is not required and only the displacement field is interpolated (Sugiyama et al.,
2006).

In this investigation, the fully parameterized ANCF beam element with twenty-four nodal
coordinates is used (Shabana and Yakoub, 2001; Yakoub and Shabana, 2001). Twelve nodal

parameters of this beam are ei
T
=
[
ri
T
ri,x
T
ri,y
T
ri,z
T
]
, where ei is the vector of nodal co-

ordinates of the i-th node, ri is the vector of the i-th node global position, while ri,k = ∂r
i/∂k

for k = x,y,z are vectors of the slope coordinates of the i-th node. The beam element used in
the study consists of two nodes, therefore, the vector of the nodal coordinates for a single-beam

element is given by eT =
[
eA
T
eB
T
]
, where A and B indicate nodes at the beam ends. The

position of an arbitrary point on the ANCF element can be obtained as follows



Volumetric locking suppression method for nearly... 979

r(x, t) =S(x)e(t) (2.1)

where x = [x y z]T, S(x) = [s1I s2I · · · s8I] is the element shape function matrix, I is a
3×3 identity matrix, and

s1 =1−3ξ
2+2ξ3 s3 = l(η− ξη) s5 =3ξ

2
−2ξ3 s7 = lξη

s2 = l(ξ−2ξ
2+ ξ3) s4 = l(ζ− ξζ) s6 = l(−ξ

2+ ξ3) s8 = lξζ
(2.2)

where l is length of the element in the undeformed state while ξ = x/l, η = y/l and ζ = z/l
are element dimensionless coordinates. It can be shown that the shape function matrix S can
describe arbitrary rigid bodymotion (Yakoub and Shabana, 2001).
The mass matrix M of the element can be written as M =

∫
V ρS

TS dV , where ρ and V
are, respectively, density and volume of the element. In theANCF, themassmatrix is constant.
In addition, the forces resulting from differentiation of kinetic energy, like Coriolis, tangential,
centrifugal and others, are equal to zero. Therefore, the only nonzero quantities in the system
equations of motion are the vectors of the elastic and external forces.
To derive the vector of the external forces, which comprise, for example, the gravitational

forces, the principle of virtual work can be used in the form δWe =F
T
e δr = F

T
eSδe = Qe

Tδe,
where δWe is the virtual work of the external force Fe and Qe is the vector of the genera-
lized external forces. For example, the nodal force vector due to gravity can be obtained as
QTeg =

∫
V F
T
egS dV , whereF

T
eg = [0 −mg 0] is the gravity force vector acting along the vertical

axis,m is total mass of an element and g is the gravitational constant.
The position vector gradient of the fully parameterized ANCF elementmay be expressed by

F = ∂r/∂X where r is given by Eq. (2.1) and X = [X Y Z]T. Using directly the expression
of the tensor F, one can evaluate the value of the strain energy for an element. In the case of
the linear-elastic materialmodel, the strain energy can bewritten asUs =

1
2

∫
V ε
TEε dV , where

ε is the strain vector associated with the Green-Lagrange strain tensor and E is the matrix of
elastic coefficients (Sopanen andMikkola, 2003). The vector of the elastic forces for an element
can be defined using the strain energyUs as follows

Qs =
(∂Us
∂e

)T
(2.3)

The present study is mainly devoted to isotropic, hyperelastic, nonlinear, and nearly incom-
pressible material models, and the proper value of the strain energy density function for these
materials is presented in the next Section of the paper.
The mass matrices and vectors of external and elastic forces of the elements follow the

standard finite element assembly procedure for each flexible body. In the case of the ANCF,
usually all the position and slope coordinates are shared between the elements. Finally, one can
write dynamic equations of motion of the constrained flexible multibody system in the general
form (Shabana, 2013)

Më+Qs+Φ
T
eλ=Qe Φ=0 (2.4)

where M, Qs and Qe are, respectively, the mass matrix, the vector of elastic forces and the
vector of external forces of the system, ë is the acceleration vector of the system, Φ represents
the vector of constraint equations (Sugiyama et al., 2003), Φe = ∂Φ/∂e is the Jacobian ma-
trix of constraints and λ is the vector of Lagrange multipliers. Equations of motion (2.4) form
a set of differential-algebraic equations with the differential index equal to 3. Finding the so-
lution to these equations is usually a more demanding task than for the solution to ordinary
differential equations (Brenan et al., 1996). Moreover, differential-algebraic equations require
special numerical techniques, as denoted by Hairer and Wanner (1996). For a review of the



980 G. Orzechowski, J. Frączek

methods used to solve Eq. (2.4), see e.g. Garćıa de Jalón and Bayo (1994). However, the most
common methods for solving differential-algebraic equations are the direct integration with,
e.g., implicit Runge-Kutta schemes (Hairer and Wanner, 1996), integration of the transformed
system with a lower index and stabilization (Gear et al., 1985) or the generalized coordinate
partitioning scheme (Wehage and Haug, 1982). In the present work, a Fortran-based research
code is used and the implicit Runge-Kutta Radau IIA scheme is utilized. It is worth noting
that if all the constrained equations were linear and time independent, system (2.4) would come
down to a set of ordinary differential equations M̂¨̂e+Q̂s = Q̂e where the vectors andmatrices
with hat are obtained after linear transformation due to constraint elimination (Garcia-Vallejo
et al., 2003).

2.1. Multi-layer beam models

Figure1presentsa two-layer beammodel.Twobeams,I and IIwith local coordinate systems
xIyIzI andxIIyIIzII (for clarity local z axes are omitted) are connected across their height. Each

Fig. 1. Two-layer beammodel

beamhas length l, widthw andheighth/2, therefore, the two-beammodel has height equal toh.
The nodes aremarked as black dots and they are shared between the elements in each layer, i.e.,
they follow the standard assembly procedure. The white dots (denoted asL1 andL2) represent
the layer connectivity points at which one can impose linear constraint equations between two
adjacent elements.X1X2X3 is the global reference frame, while the position vectors r

L2
I and r

L2
II

points to theL2 using theparameters, respectively, of elementhe I and II.The layer connectivity
constraints can be enforced at the position and slope level, thus the required continuitymight be

achieved. For exampleΦT =
[
(rL1I −r

L1
II )
T (rL2I −r

L2
II )
T
]
enforces the continuity at theposition

level only, while constraints ΦT =
[
(rL1I −r

L1
II )
T (rL1I,x−r

L1
II,x)

T (rL2I −r
L2
II )
T (rL2I,x−r

L2
II,x)

T
]
,

which are used in the numerical examples Section, impose additional constraints on the slope
along the local x axis. The total constraint vector should consist of constraints for each layer
connectivity point. This approach might be used to build a beam structures with more than
two layers across the height or to create a model with layers in both transversal directions.
In addition, the use of multi-layer structures further enables the modelling of complex systems
like vehicle tires (Patel et al., 2016). As will be shown later in the paper, the multi-layer beam
model may ensure a better convergence when the nearly incompressible materials are used. It
is worth noting that the approach for modelling a multi-layer beam used by Patel et al. (2016)
and introduced in (Liu et al., 2011) is based on subdomains with different material properties
created within a single element and cannot be easily adopted for hyperelastic and nonlinear
material models.



Volumetric locking suppression method for nearly... 981

3. Nearly incompressible polynomial material models

The hyperelastic material models, which are shown in this Section, are in the well-known form
of the Mooney-Rivlin models (Bonet and Wood, 1997; Orzechowski and Frączek, 2015). The-
se models are commonly used to represent incompressible rubber-like materials. To fulfil the
incompressibility condition, the penalty technique is employed, due to its efficiency and sim-
plicity. Therefore, any member of the Mooney-Rivlin models family can be adopted. Herein,
two simplest incompressiblemodels are used andmentioned, namely theNeo-Hookean and two-
parameter Mooney-Rivlin.

The strain-energy density function can be written for isotropic materials as a function of
the invariants of the deviatoric part of the right Cauchy-Green deformation tensor Cr =F

TF,
defined as I1 = J

−2/3I1 and I2 = J
−4/3I2, whereJ =det(F) and the invariants of the tensorCr

itself are I1 = tr(Cr) and I2 =
1
2
[tr(Cr)

2
− tr(C2r)]. In addition, the constraint J = 1 must

be ensured through the body in order to account for the material incompressibility. The strain
energy density function in the form of theMooney-Rivlin models is the following

Us =
K∑

i+j=1

µij(I1−3)
i(I2−3)

j (3.1)

where µij are material coefficients, usually determined from an experiment.K may be an arbi-
trarily large number, but in practice values ofK > 2 are rarely used. In the present paper, only
material models withK =1 are considered.

The material models from Eq. (3.1) implicitly assume that the incompressibility is ensured
by setting J equal to one. This condition is fulfilled by the penalty method (Maqueda and
Shabana, 2007;Orzechowski andFrączek, 2015). In this technique, thevolumetric energypenalty
functionUp is added to the expression of the strain energy function. This term can be expressed
as

Up =
1

2
k(J−1)2 (3.2)

where k is the penalty coefficient that represent the bulk modulus, a real material property
(Bonet and Wood, 1997). In practice, k should be selected sufficiently large to assure incom-
pressibility, but also not too large to avoid numerical complications. The use of energy function
from Eq. (3.2) with a finite coefficient k causes that the material can be considered as nearly
incompressible only.

Finally, one can combine Eqs (3.1) and (3.2) as

Usic =Us+Up (3.3)

To obtain the vector of the elastic forcesQs, the above expression should be integrated over the
flexible bodyvolumeand inserted intoEq. (2.3). Below, twomodels based on that representation
are shown.

3.0.1. Incompressible Neo-Hookean material

The incompressible Neo-Hookean is the simplest member of the Mooney-Rivlin models fa-
mily, which depends on only one elastic coefficient µ10. Therefore, the expression for the strain
energy function can bewritten asUnhs =µ10(I1−3), andµ10 is initially equal to one-half of the
shear modulus.



982 G. Orzechowski, J. Frączek

3.0.2. Incompressible Mooney-Rivlin material

A two-parameter incompressibleMooney-Rivlin material is another widely used and simple
material model which is obtained by assuming that two elastic parameters, µ10 and µ01, are not
equal to zero. Therefore, the strain energy density function takes the formUmrs =µ10(I1−3)+
µ01(I2−3). It canbe shown that in the case of small strains,Young’smodulus isE=6(µ10+µ01)
and the shear modulus is equal to µ=2(µ10+µ01) (Bathe, 1996).

4. Locking elimination techniques for nearly incompressible materials

The locking phenomena can be noticed in the case of many ANCF elements both for linear-
-elastic (Gerstmayr andShabana, 2006) andnonlinear (Orzechowski andFrączek, 2015)material
models. Its occurrence often causes an erroneously stiff bending characteristic, and is especially
noticeable for approaches that directly employ the continuummechanics.Moreover, a far greater
impact of the locking is observed in the case of an incompressible material than in the case of
compressible models. The paper presents shortly two methods that can be used for locking
suppression for incompressible material models. Both methods were introduced in the previous
work (Orzechowski and Frączek, 2015). In addition, simplifications of the projection space of
the F-bar method are introduced.

4.1. Selective reduced integration

This locking alleviation technique is commonly used to prevent Poisson’s locking in many
FEA elements (Zienkiewicz and Taylor, 2005) as well as in continuum-based ANCF elements
with a linear material model (Gerstmayr et al., 2008). In this method, the integral of the strain
energy function is split into two parts that are treated differently. In the first part, which is
fully integrated in the total element volume, one does not consider the Poisson effect, while
in the second part, which is integrated only along the beam centerline or plate midplane, i.e.
uses a reduced integration scheme, one takes the Poisson effect into account. Adequately, the
expression of the strain energy density function given by Eq. (3.1) can also be split. The first
part of theMooney-Rivlin material model, denoted asUs, can be fully integrated as it does not
consider the volumetric effect. However, the volumetric energy penalty function Up should be
considered only at the beam axis, as it accounts for the volumetric behaviour. Therefore, the
following formula is used

Usrisic =

∫

V

Usd V +A

∫

l

Up dl (4.1)

where l is length of the element, A denotes cross-section area and the index sri designates the
selective reduced integration technique.

4.2. F-bar projection method

The F-bar projection method involves product decomposition of the position vector gra-
dient into volumetric and deviatoric parts (Bonet andWood, 1997; Elguedj et al., 2008) and is
especially convenient to use when the split of the energy density function into deviatoric and
volumetric parts is not straightforward and the use of the selective reduced integration might
be troublesome. This method is a generalization of the strain projection B-bar technique to
finite-strain analysis which is considered as an extension of the selective and reduced integration
approaches (Hughes, 1987).



Volumetric locking suppression method for nearly... 983

The gradient tensor F may be split as F=FdilFdev whereFdev = J−1/3F is the deviatoric
(volume preserving) part and Fdil = J1/3I is the volumetric-dilatational part, and I is the
identity matrix. One can notice that det(Fdev) = 1 and det(Fdil) = J. Now, the tensor F can

bemodified by the use of the modified dilatational part asF
dil
= J1/3I, where

J1/3 =π(J1/3) (4.2)

and π is the linear projection operator presented in details below. Consequently, one can write

F=F
dil
F
dev =

J1/3

J1/3
F (4.3)

The modified tensor F can be directly employed to calculate the energy density function given
by Eq. (3.1). However, because Us is volume preserving, only the penalty function Up can be
affected.

4.2.1. Projection operator π

Studies performed for the isogeometric analysis beam by Elguedj et al. (2008) and further
carried out in the ANCF framework by Orzechowski and Frączek (2015) show that the L2

projection of strains is a good candidate for the projection space.Whilst the associated space on
which this projection is performed, it should have a constant value in the transversal directions.
Therefore, in the previouswork (Orzechowski andFrączek, 2015), the following lower-order basis
was employed

S̃4 =
[
1−3ξ2+2ξ3 l(ξ−2ξ2+ ξ3) 3ξ2−2ξ3 l(−ξ2+ ξ3)

]
(4.4)

where S̃4 denotes a lower order cubic basis with four components. However, even a lower order
basis may be introduced as constant, linear or quadratic (with one, two and three components,
respectively) in the longitudinal direction

S̃1 = [1] S̃2 =
[
1− ξ ξ

]

S̃3 =
[
2ξ2−3ξ+1 4ξ−4ξ2 2ξ2− ξ

] (4.5)

Next, Eq. (4.2) may be written in the new space (using any basis from Eqs (4.4) and (4.5)) as

J1/3 = S̃J̃1/3 (Elguedj et al., 2008; Hughes, 1987) where

J̃
1/3 =



∫

V

S̃
T
S̃ dV



−1∫

V

S̃
TJ1/3 dV (4.6)

The presented procedure corresponds toL2 projection of J1/3 into the S̃ basis. Next, the newly
calculated value of J1/3 may be substituted into Eq. (3.2) to obtain a modified volumetric
penalty function as

Up =
1

2
k
[
(J1/3)3−1

]2
(4.7)

Finally, the strain energydensity function for theF-bar strainprojectionmethod is expressed
as follows

UF−barsic =

∫

V

Us dV +A

∫

l

Up dl (4.8)

whereUp can be calculated only at the element centerline, as the value of J1/3 depends only on
the longitudinal coordinate.



984 G. Orzechowski, J. Frączek

5. Numerical examples

Exemplarynumerical calculations are carried outwith the fully parameterized three-dimensional
ANCF beam element with twelve nodal coordinates at each of its two nodes. In order to ef-
fectively model bodies with nearly incompressible materials, the techniques of alleviating the
volumetric locking presented in Section 4 are applied. In addition, to assemblymulti-layer beam
structures, the procedure described in Section 2.1 is used. In this study, three simple models
of highly deformable clamped beams and physical pendulums are shown. All examples use the
standard, spatial, two-node ANCF beam element with twenty-four coordinates (Shabana and
Yakoub, 2001; Yakoub and Shabana, 2001).

5.1. Physical pendulum

A dynamical analysis of a flexible beam attached at one end to the ground by a spherical
joint and falling under the gravity forces is carried out for the purpose of numerical verification.
A similar pendulumwas examined by Maqueda and Shabana (2007). In the undeformed state,
the beam has 1m in length, a square cross-section of dimension 20mm, and a material den-
sity of 7200kg/m3. In this example, the elastic coefficient for the incompressible Neo-Hookean
model is µ10 = 1MPa, while in the case of the two-parameter Mooney-Rivlin material, the
values of its coefficients are µ10 = 0.8MPa and µ01 = 0.2MPa. The material properties allow
large deformations of the body. To ensure incompressibility, the penalty term is assumed to be
k=103MPa. The pendulummodel is shown in Fig. 2. To increase body deformations, the base
of the pendulum is subjected to prescribedmotion. The constraint equations for this model can

bewritten asΦT =
[
rN1 +0.02sin(2πt) r

N
2 r

N
3

]
, where rNi for i=1,2,3 denotes the component

of the position vector of the nodeN (in meters), and t is time expressed in seconds.

Fig. 2. Very flexible physical pendulum, three-layermodel

The pendulumbody is built of four beam elements along its length. Figures 3a and 3b show
displacements of the pendulum free end tip resulting from performed dynamical simulations
for three types of models: without the locking suppression method, with the F-bar projection
and the three-layer model with the F-bar projection. Figure 3a presents results for the Neo-
Hooken material, while in Fig. 3b results for the Mooney-Rivlin material are shown. Despite
the differences inmaterials, the results in both figures are very similar. Themodels without the
locking suppressionmethodshowvery small deformationsas theybehavealmost like a rigidbody,
despite very low elastic coefficients and large pendulum density. On the other hand, all models
that employ the locking elimination by the F-bar projection show reasonable deformations.
The lack of the over-stiff response in the incompressible models during bending indicates that
the influence of the volumetric locking has been suppressed (Orzechowski and Frączek, 2015).
Therefore, it can be concluded that the influence of the volumetric locking on the results of
hyperelastic nearly incompressiblematerial modelsmight be enormous, and that themulti-layer
beam model produces acceptable results, comparable with those of the standard beam model.
The results for the selective reduced integration exhibit very similar behaviour.

Figures 4a and 4b showhowwell the incompressibility condition is preserved by the analysed
models by plotting the value of the determinant of the deformation gradient tensor J =det(F).
Figure 4a presents the value of J as a function of time at the pointA of the body (see Fig. 2).



Volumetric locking suppression method for nearly... 985

Fig. 3. Vertical position of the beam free end tip for (a) incompressible Neo-Hookean and
(b) incompressibleMooney-Rivlinmaterial: ( ) without the locking suppressionmethod, ( ) with

the F-bar projection, ( ) with the F-bar projection, three-layermodel

Fig. 4. Value of J =det(F) (a) in time and (b) along beam thickness for incompressible Neo-Hookean
material: ( ) without the locking suppressionmethod, ( ) with the F-bar projection, ( ) with

the F-bar projection, three-layermodel

Likewise, Fig. 4b presents the value of J along the line connecting pointsA andB of the body
(shown in Fig. 2) for a specific time instant (at t = 0.2s). It can be clearly seen that the
incompressibility condition for thewhole body volume is preserved only for the formulation that
does not use any locking alleviation technique.

For beams using the F-bar strain projection method, the value of J is changing noticeably,
however, in the case of the three-layer model this violation is significantly smaller than for the
one-layer beam as the region of constraint violation is bounded to the beam boundaries. Since
the lack of continuity of theJ value between the beam layers is crucial for preventing the locking
behaviour (continuous J would result in its constant value equal to one), it must be ensured
that the continuity does not occur on the gradient perpendicular to the layer boundary.

5.2. Cantilever beam

In this example, static simulation of the cantilever beam, similar to that shown in Figure 5,
is employed. The beam is clamped at one end with the help of a linear constraint equation in



986 G. Orzechowski, J. Frączek

the form of Φ= eN. The gravitational force is the only external force that acts on the system.
The material and geometrical properties of the flexible body are the same as in the previous
example, except for the size of the cross section that is increased to 40mm for both height and
width.For a verification purpose, theANCFmodelswith one, two and three layers are examined
and the results of simulations are compared with the solution of an analogous model obtained
with a commercial FEA package (ANSYS® Academic Research, 2010).

Fig. 5. Very flexible cantilever beam under gravity forces

Themainpurposeof this example is to assess the impact of thenewmodelling techniques, like
theuse of themulti-layer beamsand the loworderprojectionbasis, on the obtained solution.The
results are shown for both presented locking elimination techniques, theF-bar strain projection
and the selective reduced integration aswell as for bothmaterialmodels, the incompressibleNeo-
Hookean and two-parameter Mooney-Rivlin. For the F-bar method, four different projection
bases are employed as has been proposed in Section 4.2. The results are compared against the
FEApackage solution.For this reason,weapplied averyfinemesh in theFEApackage consisting
of 1600 higher-order solid elements (SOLID186).

In Table 1, the results of static ANCF analysis are presented for the models with one, two
and three layers across thickness, together with referenceFEA results. The results show for both
locking elimination procedures, the SRI andF-bar, the convergence to nearly the same solution
for a given number of layers and irrespective of the employed projection basis.

Table1.Deformationof the free endtipof the clampedbeamformodelswith48elements in each
layer (NoL – number of layers, ICNH – incompressible Neo-Hookean, ICMR – incompressible
Mooney-Rivlin, SRI – selective reduced integration). Results for one-layer models for SRI and
cubic F-bar locking suppressionmethods are taken from (Orzechowski and Frączek, 2015)

Method NoL
ICNH analysis ICMR analysis
X [m] Y [m] X [m] Y [m]

FEA – −0.828 −0.940 −0.828 −0.941

1 −0.894 −0.965 −0.893 −0.964
ANCF SRI 2 −0.838 −0.941 −0.838 −0.941

3 −0.829 −0.938 −0.830 −0.938

ANCF F-bar constant
projection basis

1 −0.894 −0.965 −0.893 −0.965
2 −0.840 −0.943 −0.840 −0.943
3 −0.832 −0.940 −0.832 −0.941

ANCF F-bar linear
projection basis

1 −0.894 −0.965 −0.893 −0.965
2 −0.837 −0.941 −0.838 −0.941
3 −0.829 −0.938 −0.830 −0.939

ANCF F-bar quadratic
projection basis

1 −0.892 −0.963 −0.891 −0.962
2 −0.837 −0.941 −0.838 −0.941
3 −0.829 −0.938 −0.830 −0.938

ANCF F-bar cubic
projection basis

1 −0.891 −0.962 −0.890 −0.962
2 −0.837 −0.941 −0.838 −0.941
3 −0.829 −0.938 −0.830 −0.938



Volumetric locking suppression method for nearly... 987

The results presented in Table 1 show that the convergent solutions given byANCF are no-
ticeably different when different numbers of layers are used. For the one-layer model, theANCF
results are significantly different from the reference. The noticeable discrepancy between the
reference FEA results and the ANCFmodels is mainly due to violation of the incompressibility
assumption at the beam surface. Nevertheless, this behaviour is consistent with the characteri-
stics of the finite element analysis, where incompressibility is ensured only on the average (Adler
et al., 2014). The standard approach to deal with such a problem is to employ a larger number
of elements across thickness, and the results inTable 1 show that the beammodelswith two and
three layers conformmuch better to the reference. Therefore, when the incompressible material
model employs the locking elimination technique based on reduced integration, it is crucial to
apply more elements across the transversal direction in order to obtain reasonable results. In
addition, in the case of the F-bar method, the simplest constant projection basis is able to
reproduce a proper solution that is only slightly different than the reference one, while offering
the smallest computational complexity.

5.3. Cantilever rubber-like beam

The last example is the dynamic analysis of the clampedbeammodel, wherein the constraint
equations are identical as those in thepreviousSection.The system is similar to themodel shown
by Jung et al. (2011). The body undergoes the gravitational force, and is made of a rubber-like
material having density of 2150kg/m3,Kirchhoff’smodulus of 1.91MPaand the penalty termof
103MPa. The beamhas a rectangular cross section of 7mmwidth and 5mmheight, and 0.35m
in length. Themodels use thirty fully parameterizedANCFbeamelements for each layer, whilst
one- and three-layer models are compared. Such a large number of elements is needed to obtain
a converged solution for both locking elimination procedures. The results of ANCF simulations
are compared with the results obtained with the FEA package. The beam elements from the
classic FEA package cannot be used with nonlinear incompressible material models. Therefore,
SOLID185 (ANSYS® Academic Research, 2010) elements with reduced integration are used.
To obtain a convergent solution, the use of 336 solids is sufficient. In this example, the only
considered material model is the incompressible Neo-Hookean.

Fig. 6. Beam end tip displacements, (a) one-layermodel, (b) three-layermodel: ( ) reference FEA,
( ) ANCF with selective reduced integration, ( ) ANCF full integration withF

Figure 6a shows the results for twoANCFmodelswith different locking suppressionmethods
and for the reference FEA solution (Orzechowski and Frączek, 2015). The results are recalled
here for easier comparison with the three-layer model. Both ANCF models provide almost the
same results, and for the one-layer model a reasonably good agreement with the reference can
be observed.However, the results inFig. 6b show a large improvement for the three-layermodel,
irrespective of the employed projection basis. The results for the SRI present a similar advance.



988 G. Orzechowski, J. Frączek

One can conclude that the alleviation of the volumetric locking effect and introduction of the
multi-layer structures leads to a substantial improvement as comparedwith the one-layer model
and ensures the quality of the results obtained using higher-order beam elements with forty-
-two nodal parameters as shownbyOrzechowski andFrączek (2015). In addition, in comparison
with the FEA solidmodel, the number of parameters involved is, in the given application, more
than four times less in the ANCF one-layer beam model, and two times less in the case of the
three-layer model.

6. Conclusions

TheANCF is used successfully in the analysis of flexible bodies undergoing large deformations.
The current study presents that nonlinear and nearly incompressible material models can be
included in theANCF in a straightforwardway. The implementation of nonlinear incompressible
materials and locking elimination techniqueswithin theANCF frameworkpresented so far in the
literature exhibit inaccurate behaviour in bending-dominated examples for a standard twenty-
-four parameter three-dimensional beam element.

In order to carry out accurate and efficient model simulations with incompressibility, two
locking elimination techniques are presented – the F-bar strain projection method and the
selective reduced integration. Both methods are implemented in the ANCF framework and
tested withmulti-layer beam structures. In addition, for theF-barmethod, different projection
spaces are investigated. Numerical tests show that the locking influence for an incompressible
materialmodels is enormous. The results of ANCF simulations are comparedwith a commercial
FEA package and a very good agreement is found, especially for three-layer models. In the case
of the F-bar method, the order of the projection space has aminor influence on the quality of
the solution.

In the classical FEA, solid elements must be used when the materials with a nonlinear
characteristic areused tomodel slender structures. In contrast, in theANCF, fullyparameterized
beam elements can be employed in such systems, the result of which are models withmuch less
parameters than for theFEA.However, the application ofANCFbeamsmodelledwithnonlinear
materials to complex cross-sectional shapes might require further investigations (Orzechowski,
2012; Orzechowski and Shabana, 2016).

Adesirabledirection for future research is thedevelopment of alternativemethods thatwould
enable efficient locking elimination. In addition, it is indicated in the literature that, applying
the systemswhich use incompressiblematerial models, onemay obtain solutions that are highly
imprecise in stresses, although a proper solution exist for displacements (Bathe, 1996). Hence,
the distribution of stresses in such models should also be taken into account as an essential
research topic.

Acknowledgements

This work has been supported under grantNo. DEC-2012/07/B/ST8/03993by the National Science

Centre.

References

1. Adler J.H., Dorfmann L., Han D., MacLachlan S., Paetsch C., 2014, Mathematical
and computational models of incompressible materials subject to shear, IMA Journal of Applied
Mathematics, 79, 5, 889-914

2. ANSYS®Academic Research, 2010,ANSYS Mechanical APDL Documentation, ANSYS Inc.,
release 13.0



Volumetric locking suppression method for nearly... 989

3. Bathe K.J., 1996,Finite Element Procedures, Prentice Hall, New Jersey

4. Bonet J., Wood R.D., 1997, Nonlinear Continuum Mechanics for Finite Element Analysis,
Cambridge University Press, Cambridge; NewYork, NY, USA

5. Brenan K.E., Campbell S.L., Petzold L.R., 1996, Numerical Solution of Initial-Value Pro-
blems in Differential-Algebraic Equations, SIAM, Philadelphia

6. ElguedjT.,BazilevsY.,CaloV.M.,HughesT.J., 2008,B-barandF-barprojectionmethods
for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS
elements,Computer Methods in Applied Mechanics and Engineering, 197, 33-40, 2732-2762

7. Garćıa de Jalón J., Bayo E., 1994,Kinematic and Dynamic Simulation of Multibody Systems.
The Real-Time Challenge, Springer, NewYork

8. Garcia-Vallejo D., Escalona J.L., Mayo J., Doḿinguez J., 2003, Describing rigid-flexible
multibody systems using absolute coordinates,Nonlinear Dynamics, 34, 1-2, 75-94

9. Gear C.W., Leimkuhler B., Gupta G.K., 1985, Automatic integration of Euler-Lagrange
equations with constraints, Journal of Computational and Applied Mathematics, 12-13, 77-90

10. Gerstmayr J., Matikainen M.K., MikkolaA.M., 2008,A geometrically exact beam element
based on the absolute nodal coordinate formulation,Multibody System Dynamics, 20, 4, 359-384

11. Gerstmayr J., ShabanaA.A., 2006,Analysis of thin beams and cables using the absolute nodal
co-ordinate formulation,Nonlinear Dynamics, 45, 1, 109-130

12. Hairer E.,WannerG., 1996,Solving Oridinary Differential Equations II. Stiff and Differential-
-Algebraic Problems, Springer, Berlin, 2nd edition

13. Hughes T.J.R., 1987, The Finite Element Method. Linear Static and Dynamic Finite Element
Analysis, Prentice Hall, New Jersey

14. Jung S., Park T., Chung W., 2011, Dynamic analysis of rubber-like material using absolute
nodal coordinate formulation based on the non-linear constitutive law, Nonlinear Dynamics, 63,
1, 149-157

15. LiuC.,TianQ.,HuH., 2011,Dynamics of a large scale rigid-flexiblemultibody systemcomposed
of composite laminated plates,Multibody System Dynamics, 26, 3, 283-305

16. Maqueda L.G., Mohamed A.N.A., Shabana A.A., 2010, Use of general nonlinear material
models in beam problems: Application to belts and rubber chains, Journal of Computational and
Nonlinear Dynamics, 5, 2, 021003, 10 pages

17. Maqueda L.G., Shabana A.A., 2007, Poissonmodes and general nonlinear constitutive models
in the large displacement analysis of beams,Multibody System Dynamics, 18, 3, 375-396

18. Orzechowski G., 2012, Analysis of beam elements of circular cross section using the absolute
nodal coordinate formulation,Archive of Mechanical Engineering,LIX, 3, 283-296

19. Orzechowski G., Frączek J., 2015, Nearly incompressible nonlinear material models in the
large deformation analysis of beams using ANCF,Nonlinear Dynamics, 82, 1, 451-464

20. Orzechowski G., Shabana A.A., 2016, Analysis of warping deformation modes using higher
order ANCF beam element, Journal of Sound and Vibration, 363, 428-445

21. Patel M., Orzechowski G., Tian Q., Shabana A.A., 2016, A new multibody system ap-
proach for tire modeling using ANCF finite elements,Proceedings of the Institution of Mechanical
Engineers, Part K: Journal of Multi-body Dynamics, 230, 1, 69-84

22. Shabana A.A., 1997a, Definition of the slopes and the finite element absolute nodal coordinate
formulation,Multibody System Dynamics, 1, 3, 339-348

23. Shabana A.A., 1997b, Flexible multibody dynamics: review of past and recent developments,
Multibody System Dynamics, 1, 2, 189-222

24. Shabana A.A., 2008, Computational Continuum Mechanics, Cambridge University Press, Cam-
bridge, first edition



990 G. Orzechowski, J. Frączek

25. Shabana A.A., 2013,Dynamics of Multibody Systems, Cambridge University Press, 4th edition

26. Shabana A.A., Yakoub R.Y., 2001, Three dimensional absolute nodal coordinate formulation
for beam elements: theory, Journal of Mechanical Design, 123, 4, 606-613

27. Sopanen J.T., Mikkola A.M., 2003, Description of elastic forces in absolute nodal coordinate
formulation,Nonlinear Dynamics, 34, 1, 53-74

28. Sugiyama H., Escalona J.L., Shabana A.A., 2003, Formulation of three-dimensional joint
constraints using the absolute nodal coordinates,Nonlinear Dynamics, 31, 2, 167-195

29. Sugiyama H., Gerstmayr J., Shabana A.A., 2006, Deformation modes in the finite element
absolute nodal coordinate formulation, Journal of Sound and Vibration, 298, 4-5, 1129-1149

30. WasfyT.M.,NoorA.K., 2003,Computational strategies for flexiblemultibody systems,Applied
Mechanics Reviews, 56, 6, 553-613

31. Wehage R.A., Haug E.J., 1982,Generalized coordinate partitioning for dimension reduction in
analysis of constrained dynamic systems, Journal of Mechanical Design, 104, 1, 247-255

32. YakoubR.Y., ShabanaA.A., 2001,Threedimensional absolutenodal coordinate formulation for
beam elements: implementation and applications, Journal of Mechanical Design, 123, 4, 614-621

33. Zienkiewicz O.C., Taylor R.L., 2005, The Finite Element Method for Solid and Structural
Mechanics, Butterworth-Heinemann, Oxford, 6th edition

Manuscript received September 19, 2016; accepted for print April 5, 2017