Nearly incompressible nonlinear material models in the large deformation analysis of beams using ...

14 Pages·2015·1.28 MB·English

by  Grzegorz Orzechowski

NonlinearDyn(2015)82:451–464
DOI10.1007/s11071-015-2167-1
ORIGINAL PAPER
Nearly incompressible nonlinear material models
in the large deformation analysis of beams using ANCF
GrzegorzOrzechowski · JanuszFra˛czek
Received:4February2015/Accepted:15May2015/Publishedonline:2June2015
©TheAuthor(s)2015.ThisarticleispublishedwithopenaccessatSpringerlink.com
Abstract Many modern applications of the flexible show important improvement as compared with pre-
multibodysystemsrequireformulationsthatcaneffec- viousapplicationsandgoodagreementwithreference
tivelysolveproblemsthatincludelargedisplacements results.
anddeformationshavingtheabilitytomodelnonlinear
materials. One method that allows dealing with such Keywords Flexible bodies · Absolute nodal
systemsiscontinuum-basedabsolutenodalcoordinate coordinateformulation·Incompressiblematerials·
formulation (ANCF). The objective of this study is Lockingphenomena·Selectivereducedintegration·
to formulate an efficient method of modeling nonlin- F-bar
ear nearly incompressible materials with polynomial
Mooney–Rivlinmodelsandvolumetricenergypenalty 1 Introduction
functionintheframeworkoftheANCF.Themainpart
of this paper is dedicated to the examination of sev- Inthemoderndesignofthesoftwareandapplications,
eral ANCF fully parameterized beam elements under thehighlyflexiblebodiesbuiltwithnonlinearmateri-
incompressibleregime.Moreover,twovolumetricsup- als are in common use. Consequently, the multibody
pression methods, originating in the finite element system (MBS) simulation software should be able to
analysis,areproposed:awell-knownselectivereduced include such bodies in a reliable manner. The most
integration and F-bar projection. It is also presented commonlyusedmethodinflexiblebodysimulationsis
thattheuseofthesemethodsiscrucialforperforming thefloatingframeofreferenceformulation[33].This
reliable analysis of models with bending-dominated method is usually able to include only small defor-
loads when lower-order elements are employed. The mations and linear material models. The large defor-
results of the simulations carried on with considered mationsandnonlinearitiesmightbemodeledwiththe
elementsandproposedmethodsarecomparedwiththe finiteelementanalysis(FEA)[4],buttheFEAdoesnot
resultsobtainedfromcommercialfiniteelementpack- workperfectlywiththeMBS[42].
age and existing ANCF implementation. The results Absolutenodalcoordinateformulation(ANCF)[32]
is one of the methods that can be used for efficient
B
G.Orzechowski( )·J.Fra˛czek dynamic simulation of flexible multibody systems.
InstituteofAeronauticsandAppliedMechanics,Facultyof
This formulation has some unique features allowing
PowerandAeronauticalEngineering,WarsawUniversity
formodelingnonlinearmaterialmodelswithbeamele-
ofTechnology,Nowowiejska24,00-665Warsaw,Poland
e-mail:[email protected] ments in a straightforward manner. In the ANCF, to
describe local orientation, the slope coordinates are
J.Fra˛czek
e-mail:[email protected] usedinsteadofrotations.Thisallowscomplexshapes
123
452 G.Orzechowski,J.Fra˛czek
to be represented using a small number of elements. [41];however,thistopicwillnotbefurtherexamined
Theelasticforcesofthebeamelementscanbederived inthepaper.
withnotonlytheclassicalbeamtheories(suchasthe The main objective of this paper is to examine
Bernoullibeammodel)butalsothecontinuummechan- theuseofnonlinear,hyperelastic,nearlyincompress-
icsapproachmaybeemployed.Thisallowsforuseof ible material models using the absolute nodal coor-
general constitutive equations for the nonlinear (both dinate formulation beam elements. While this objec-
compressibleandincompressible)materialmodels. tive was addressed in previous publications [23], the
The absolute nodal coordinate formulation is well main contribution of this paper is to show the impor-
suitedtomodelflexiblebodiesintheframeworkofthe tance of volumetric locking elimination in the appli-
multibodysystems.TheANCFcaneasilybeusedwith cations when bending-dominated models are consid-
the nonincremental formulation [34]. Flexible ANCF ered.Foradetailedoverview ofthelockingphenom-
bodiesmayundergolargedeformations,aswellaslarge enaoccurringwithinANCFframeworkaswellasgen-
translational and rotational displacements. Addition- eralreviewofdevelopmentsinANCF,consult[14,26].
ally, it allows for exact modeling of the rigid body Moreover, exemplary techniques that may reduce the
motions.ThischaracteristicmakestheANCFmethod lockinginfluenceontheresultsarepresentedandval-
well suited for the dynamic simulations of the slen- idated with several numerical examples. Two basic
der, flexible bodies modeled with nonlinear material incompressiblematerialmodelsarepresented,incom-
models. pressible Neo-Hookean and two-parameter Mooney–
Efficient and reliable modeling of the nonlinear, Rivlin [34]; however, presented methods may be eas-
incompressiblematerialsisanimportanttaskinmod- ilyextendedtothehigher-ordermaterialmodels.The
ernindustrialandengineeringapplications.Manyareas incompressibility of the materials is ensured in the
ofengineeringandbiomechanicsaccountforthemate- current study by the penalty method, which is cho-
rial incompressibility. For example, the vessel wall sen due to its popularity and simplicity. Furthermore,
material,heartmuscle,aswellassomebioelastomers, fourfullyparameterizedelementsareemployedasthey
suchaselastinandresilin,aretreatedasincompressible allowforastraightforwardimplementationofthecon-
materials [10]. The rubber-like materials are used in stitutive nonlinear models. The first element is stan-
manypracticalapplicationssuchasautomotive(tires, dardbeamelementwithtwenty-fournodalcoordinates
rubber seals, boots), defense (coating for equipment, introducedin[35,45]thatisreferredinthefollowing
radar absorbing materials), safety (seat belts, impact as beam element 24 (BE24). The second element is
absorbers) and many others. Incompressible materi- higher-order beam with thirty nodal coordinates that
als have to be modeled in commercial finite element suppress shear locking introduced in [13]. This ele-
packageswithshellorsolidelements[2],evenforthe mentisdenotedinthepaperasbeamelement30with
slender structures that are usually approximated with shearlockingsuppression(BE30S).Thethirdelement
beams.TheuseoftheANCFinsuchapplicationallows is higher-order element with thirty nodal coordinates
overcomingthislimitationbyusingthefullyparame- thatincludetrapezoidaldeformationmodeofthecross
terizedbeamelements[34]. sectionthatsuppressesvolumetriclockingintroduced
Most of the papers devoted to the ANCF assume in [24,25]. This element is referred in the following
linearconstitutiverelationbetweenstressesandstrains asbeam element 30withtrapezoidalmode(BE30T).
[26]. The isotropic and hyperelastic material mod- The last element is higher-order beam element with
els(includingincompressiblemodels)wereintroduced forty-twonodalcoordinatesthatincludefullquadratic
intoANCFbyMaquedaandShabana[23].Moreover, approximation in transversal direction introduced in
theapplicationsoftheincompressiblematerialstothe [36]whereitisdenotedasB2andfurtherinvestigated
rubber chains build with ANCF beam elements were in [28]. This element is denoted in the following as
presented[22].Inaddition,Jungetal.[18]presentthe beamelement42(BE42).SimilarlylikeBE30T,BE42
experimentalvalidationoftherubber-likebeammod- suppresses the volumetric locking phenomena using
eledwithANCFelements,althoughtheydonothigh- higher-orderapproximationinthecrosssection.
lighttheimportanceofthevolumetriclockingsuppres- Withinthefiniteelementmethod,awidevarietyof
sionmethods.Furthermore,theformulationsofplastic- element technologies is known which main goal is to
itywereexaminedwithintheframeworkoftheANCF obtain a general element formulation that provides a
123
Nearlyincompressiblenonlinearmaterialmodels 453
goodbendingbehavior,nolockingforincompressible issuesconnectedwithvolumetriclockingphenomena
materialsandverythinelements,efficiency,distortion is depicted in Sect. 4. This section also presents the
insensitivityandgoodcoarsemeshaccuracy[44].The proposedmethodsforsuppressingthevolumetriclock-
notable and important methods that are used to solve ing.InSect.5,severalnumericalexamplesofthemod-
theissuesassociatedwithlockingsaretheselectiveand els that incorporate incompressible material with and
reducedintegrationtechniques[21]whichareeffective withoutlockingeliminationtechniquesarepresented.
andeasytoimplement.Alternatively,afamilyofhybrid ConclusionsofthisstudyarepresentedinSection6.
andmixedfiniteelementshasbeenproposed[5].The
mostcommonlyknownmixedformulationsarebased
on the two-field variational Hellinger–Reissner prin-
ciple [30] and three-field Hu–Washizu principle [43]; 2 Kinematicsofthedeformablebodies
however, the development of the finite elements via usingtheANCF
these principles requires the approximation of addi-
tionaltensorfields,andtheresultingnumberofscalar In the absolute nodal coordinate formulation, nodal
variables might be very large. Therefore, these prin- coordinatesconsistofthetranslationalandslopescoor-
ciplesareseldomimplementeddirectly[5].Notewor- dinates,describedwithrespecttotheglobalcoordinate
thy, the shear locking suppression methods based on system. There are no rotational coordinates, such as
these variational principles were employed in elastic EulerparametersorEulerangles,usedtodescribethe
line and mid-surface approaches with beam [31] and orientationoftheelement.Duetothisfeature,onlythe
plate [7] elements. More recently, the simpler, pure displacementfieldisinterpolated,andthereisnoneed
displacement approaches could extend the selective fortheindependentrotationfieldinterpolation[40].For
reduced integration method to more complex cases example,thenodeofthestandardbeamelementBE24
emerged.Amongthemarethe B-barstrainprojection hasthefollowingcoordinates:
method [16,37] and an enhanced assumed strain for- (cid:2) (cid:3)
T
mulation [38]. An important advantage of the B-bar ei = riT ri T ri T ri T (1)
,x ,y ,z
approachovertheenhancedformulationisthatitonly
results in a modified strain–displacement matrix with where vector ei contains nodal parameters of the i
no additional degrees of freedom or parameters [5]. node, vector ri is the global position vector of the i
The enhanced assumed strain method was applied to node, while vectors ri = ∂ri/∂k for k = x,y,z are
,k
theANCFbeamelementtoalleviatelockingsbyGer- first-ordergradientsofthepositionvector(slopecoor-
stmayrandShabana[12],andmorerecently,itwasused dinates) of the i node. The vector of the nodal coor-
(cid:2) (cid:3)
withshellelementbyYamashitaetal.[46].Inthecur- dinates for single BE24 element is e = eAT eBT T
rentpaper,itisshownthatthevolumetriclockingsup-
where A and B denote element nodes. All elements
pression iscrucialfor theBE24 and BE30S elements
considered in the paper are beam elements with two
and two methods of the locking alleviation are pre-
nodes;however,theyemployhigher-ordergradientsas
sented.Thefirstoneisbasedontheselectivereduced
nodalcoordinates.Fordetaileddescriptionofthekine-
integration,anditissimilartothePoissonlockingelim-
matics of these elements, see [45] for BE24, [13] for
inationmethodusedwithlinearmaterialswithinANCF
BE30S,[25]forBE30Tand[36]forBE42.Usingthe
framework[11].ThesecondmethodisbasedontheB-
vectorofnodalcoordinatese,onecangettheposition
barstrainprojectiontechnique.Bothmethodsarepre-
ofthearbitrarypointontheelement:
sentedindetailsandverifiedwithnumericalexamples.
Thepaperisorganizedasfollows.Section2presents r = Se (2)
the basic kinematic equations of the spatial absolute
nodalcoordinateformulationbeamelements.Inaddi- whereSisthematrixthatcontainselementshapefunc-
tion, the standard form of the equations of motion is tionsanddependsonlyonthespatialcoordinates.Usu-
shown. Section 3 introduces the well-known incom- ally, shape functions are written as a function of ele-
pressible material models, the penalty method that mentdimensionlesscoordinatesS= S(ξ,η,ζ),where
ensurestheincompressibilityandmodifiedversionof ξ = x/l,η= y/landζ =z/lforldenotingthelength
thenonlinearmaterialmodels.Theformulationofthe oftheelementintheundeformedstate.
123
454 G.Orzechowski,J.Fra˛czek
Thekineticenergyofthesingleelementisdefined of the Mooney–Rivlin polynomial are presented and
bythefollowingformula: twospecialcasesofthismodelaredescribed,namely
(cid:4)
1 1 the incompressible Neo-Hookean and two-parameter
T = ρr˙Tr˙dV = e˙TMe˙ (3) incompressibleMooney–Rivlinmaterialmodels.
2 2
V
whereρandV are,respectively,densityandvolumeof
theelementand M isthemassmatrixoftheelement. 3.1 StandardMooney–Rivlinpolynomialmodels
After differentiation of Eq. (2) with respect to time,
massm(cid:4)atrixcanbewrittenasfollows: Fullyparameterizedbeamandplateelementshavethe
followingexpressionforthematrixofthepositionvec-
M = ρSTSdV (4)
torgradient:
V
ThemassmatrixofANCFbeamelementsisconstant. ∂r
F = (8)
Intheabsolutenodalcoordinateformulation,Corio- ∂x
lis,tangential,centrifugal,aswellasotherforcesresult- wherer isdefinedbyEq.(2)andx =[x y z]T.
ingfromdifferentiationofthekineticenergyareequal
Thestrainenergydensityfunctionfortheisotropic
tozero.Thus,thenonzerotermsinsystemequationsof
materials depends only on the invariants of the right
motioncomefromthevectorsoftheelasticandexter- Cauchy–GreendeformationtensorC = FTF,defined
r
nal forces. The vector of the external forces, which
as:
include the gravitational forces, can be derived using
theprincipleofvirtualworkasfollows: I1 =tr((cid:2)Cr) (cid:3) (9)
1
δWe = QeTδe (5) I2 = 2 tr(Cr)2−tr(Cr2) (10)
where δWe is the virtual work of the external forces I3 =det(Cr) (11)
and Qeisthevectorofthegeneralizedexternalforces. Ifthematerialisincompressible,theconstraint J =
Vectoroftheelasticforcesisgivenby: 1, where J = det(F), or identically I = J2 = 1,
(cid:5) (cid:6) 3
∂U T needstobeensurethroughthematerial.Therefore,the
Q = s (6)
s ∂e strain energy density function for the incompressible
materials depends only on the I and I . A general
whereU istheenergyoftheelasticdeformation(for 1 2
s
form of the strain energy function for incompressible
linearelasticmaterialmodel)orelasticenergydensity
rubber-likematerialsintheformoftheMooney–Rivlin
function (e.g., for the isotropic, hyperelastic material
modelsisfollowing:
models). Elastic energy of the linear elastic material
(cid:7)∞ (cid:7)∞
modelmaybefoundintheliterature[34].Thepresent
U = μ (I −3)i(I −3)j, μ =0 (12)
paperismainlyfocusedonthenonlinear,hyperelastic, s ij 1 2 00
i=0 j=0
isotropicandnearlyincompressiblematerialmodels.
Finally, one can write the equations of motion in where coefficients μij are material constants, usually
theindependentcoordinatesforthethree-dimensional determinedfromtheexperiment.
ANCFbeamasfollows: MaterialmodelsfromEq.(12)assumethatmaterial
isfullyincompressible,thatis, J isequaltoone.This
Me¨+ Q = Q (7)
s e condition must be fulfilled by using additional tech-
niques,suchasLagrangemultipliersorpenaltymethod
[23].Inthecurrentpaper,thepenaltymethodisused,
3 Existingnonlinearconstitutivemodels
duetosimplicityandefficiency.Inthepenaltymethod,
oftheincompressiblematerials
additionaltermsareincludedintotheexpressionofthe
strainenergyfunction:
Inthecurrentsection,apolynomialhyperelasticmate-
rialmodelisshownintheformoftheMooney–Rivlin U = 1k(J −1)2 (13)
p
models [34]. The incompressibility condition is sat- 2
isfied by the technique of penalty method. In addi- where U is volumetric energy penalty function and
p
tion, the required modifications to the original model k is penalty coefficient, which must be selected large
123
Nearlyincompressiblenonlinearmaterialmodels 455
enough to ensure the incompressibility. However, k
must not be too large because this can cause serious
numericalissues.Itshouldbepointedoutthattheuse
ofthevolumetricenergyfunctionU causesmaterial Fig.1 Veryflexibleclampedbeamundergravityforces
p
to be considered as nearly, not fully, incompressible. where μ = 0.5μ, for μ initially equal to the shear
10
In addition, penalty term k describes a real material
modulus.
property,thebulkmodulus[6].
Finally, the strain energy density function of the
3.4 IncompressibleMooney–Rivlinmaterial
incompressibleMooney–Rivlinmodelsiswrittenas:
U =U +U (14) Anothersimpleandcommonlyusedmaterialmodelis
sic s p
two-parameterincompressibleMooney–Rivlinmater-
Inthecurrentstudy,twomodelsbasedonthatrep-
ial,whichdependsontwoelasticcoefficients,μ and
resentationarepresented. 10
μ . Expression for the strain energy function of the
01
Mooney–Rivlinmaterialmaybewrittenasfollows:
3.2 ModifiedMooney–Rivlinmodels U =μ (I¯ −3)+μ (I¯ −3) (19)
mr 10 1 01 2
The Mooney–Rivlin models expressed by Eq. (12) One can notice for μ01 = 0 this model reduces to
accountsforthevolumetric(pressure)behavior,which previouslyintroducedNeo-Hookeanmaterial.Inaddi-
may cause numerical problems during computations. tion, if small strains are considered, the shear modu-
Therefore, it is necessary to separate the volumetric lusisμ = 2(μ10 +μ01)andtheYoung’smodulusis
anddeviatoriccomponentsofthedeformation[6].To E =6(μ10+μ01)[4].
uncouplethisbehavior,insteadoftheinvariantsofthe
rightCauchy–Greendeformationtensor,theinvariants 4 Lockingeliminationmethodsforincompressible
ofthedeviatoricpartofthattensorareused: materials
I¯ = J−2/3I (15)
1 1
I¯ = J−4/3I (16) It was reported in the literature [13] that the locking
2 2
phenomena affect many ANCF elements with linear
The final form of the Mooney–Rivlin models for elastic material model. The locking behavior is espe-
incompressibleandhyperelasticmaterialsusedinthat ciallynoticeablewhencontinuummechanicsapproach
paper is given by substitution of the invariants from isusedandrevealedthatflexiblebodieshaveanerro-
Eqs.(15)and(16)intoEq.(14): neouslystiffercharacteristicduringbending.However,
(cid:7)∞ (cid:7)∞ in the case of the incompressible material model, the
Usic =Us +Up = μij(I¯1−3)i(I¯2−3)j lockinginfluencecanbemuchlargerthanincaseofthe
i=0 j=0 compressiblematerials,aswillbeshowninthecurrent
1 section.
+ k(J −1)2, μ =0 (17)
00
2
Thevalueoftheelasticforces Qs maybeobtained 4.1 Problemformulation
according to Eq. (6), after integration of the above
expressionoverthevolumeoftheflexiblebody. Theuseofhyperelastic,nearlyincompressiblematerial
models [with strain energy density function given by
Eqs.(18)or (19)]mayleadtoresultsthatdonotcon-
3.3 IncompressibleNeo-Hookeanmaterial
vergetothecorrectsolution.Thisproblemisdenotedin
theliteratureasvolumetriclockingphenomena[6,34].
The simplest material model prescribed by Mooney–
Thefollowingexamplepresentsthisbehavior.
RivlinmodelsistheincompressibleNeo-Hookean.In
As an illustration of the problem, the highly flex-
thatmodel,onlythecoefficientμ doesnotequalzero,
10 ible, clamped beam under gravity forces is analyzed.
sothestrainenergyfunctionhasthefollowingform:
Thebeam,asshownonFig.1,isof1mlongandhas
U =μ (I¯ −3) (18) a square cross section of 20 mm in size. The model
nh 10 1
123
456 G.Orzechowski,J.Fra˛czek
1.2 numerically. Moreover, in the finite element analysis,
1 usually the Gaussian quadrature is used to integrate
0.8
overthecuboidelementvolume[4,34]:
m] 0.6 (cid:4)
[
nt 0.4 G(ξ,η,ζ)dV
e
em 0.2 V(cid:4) (cid:4) (cid:4)
ac 0 1 ηmax ζmax
pl = G(ξ,η,ζ)dξdηdζ
Dis--00..42 0 ηmin ζmin
(cid:7)nξ (cid:7)nη (cid:7)nζ
-0.6
≈ w w w G(ξ ,η ,ζ ) (20)
-0.8 i j k i j k
0 0.2 0.4 0.6 0.8 1 i=1 j=1k=1
Time[s]
where G is any scalar, vector or matrix function of
Fig.2 Lockinginfluenceontheincompressiblematerialmodel. spatialcoordinatesξ,η andζ,wi,wj andwk arethe
( ) Compressible Neo-Hookean, ( ) Incompressible weightsofthequadrature,ξ ,η andζ arethepoints
i j k
Neo-Hookean
ofthequadratureandnξ,nηandnζ aretheordersofthe
quadratureindirections,respectively,ξ,η andζ.The
usessixfullyparameterizedANCFthree-dimensional numberofintegrationpointsineachdirectionmightbe
beamelementsBE24andincludestwomaterialmod- different.Ingeneral,thequadratureofthefunctionGis
els:compressibleNeo-Hookean[23]andincompress- onlyanapproximationofitsintegral.However,when
ible Neo-Hookean [see Eq. (18)]. The values of the functionGisgivenasapolynomial(whichisusuallya
Lamé parameters, used for the compressible Neo- casefortheANCFelements),theordersofthequadra-
Hookean, are λ = 1.15 MPa as the first Lamé para- turenξ,nη andnζ canbeselectedinthewaythatwill
meterandμ = 0.769MPaasthesecondLamépara- guaranteetheexactevaluationoftheintegral(thatorder
meter, called the shear modulus. Selected values of dependsonthedegreesofthepolynomialinappropri-
the parameters allow for large deformations. For the atedirections).Detailedinformationabouttherequired
compressible Neo-Hookean, the material constant is numberofintegrationpoints,typesofquadraturesand
μ10 =0.333MPaandpenaltyconstantisk =109Pa. dependencies between the polynomial order and the
Thematerialdensityisassumedtobe7200kg/m3. quadratureordercanbefoundintheliterature[12].
In the Fig. 2, one can see the beam modeled In the case of reduced integration, the quadrature
with incompressible Neo-Hookean material, which is ordersareselectedaslowerthanthoserequiredforfull
describedinSect.3.3,exhibitsaverystiffresponse— integration. In the selective reduced integration, only
the tip displacements for the compressible model are some parts of the spatial function G are evaluated by
much larger. This behavior is caused by the volumet- meansofthereducedintegration.
riclocking phenomena, which affects twoconsidered Theselectivereducedintegrationtechniqueisused
elements: BE24 and BE30S as they do not include a topreventPoissonlockingintheANCFelementswitha
trapezoidalcross-sectionaldeformationmode[25]. continuummechanicsapproachandlinearelasticmate-
Thepresentedexampleindicatesthatinsimulations rialmodel[11],aswellasformanyFEAelements[47].
involving bending, the locking effect should be sup- Thismethod’sideaistosplittheintegralofthestrain
pressedortherelevantelementthatisnotaffectedby energy function into two parts. The first part, which
thelockingshouldbeconsidered.Thepaperdiscusses does not take into account the Poisson effect, is fully
twomethodsthatcanbeusedforlockingeliminationin integrated.Thesecondpart,whichconsidersthePois-
connectionwithincompressiblematerialmodels.One soneffect,usesareducedintegrationscheme.
isbasedontheselectivereducedintegration[11],and Theexpressionforthestrainenergydensityfunction
thesecondisbasedonthe F-barmethod[8]. of the incompressible Mooney–Rivlin models, given
by Eq. (17), can also be split adequately. The first
4.2 Selectivereducedintegration part,denotedasU ,doesnotaccountforthevolumet-
s
ric behavior; therefore, it is fully integrated. The sec-
TheintegrationoftheintegralspresentedinSect.2,in ond part, the volumetric energy penalty function U ,
p
particulartheelasticforcesvector,isusuallyperformed accountsforthevolumetriceffectandisonlyconsid-
123
Nearlyincompressiblenonlinearmaterialmodels 457
eredatthebeamaxis.Thestrainenergydensityvalue toincludethefactorsderivedfrom F intheU partof
s
canbeexpressedbythefollowingformula: theenergydensityfunctiongivenbyEq.(17),asthat
(cid:4) (cid:4)
partisvolumepreserving.
USRI = U dV + A U dl (21)
sic s p
V l
4.3.1 Projectionoperator
whereindexSRIdenotestheselectivereducedintegra-
tion scheme, while A andl are, respectively, the area
Theuseofthestrainprojectionmethodrequiresfurther
andthelengthoftheelement,whichindicatesthatfunc-
investigations on the choice of the projection opera-
tionU isevaluatedonlyalongbeamcenterline.This
p
torandthespaceonwhichtheprojectionwillbeper-
alsomeansthatintheformulaforquadraturegivenby
formed. In the original F-bar method, the modified
Eq.(20),theordersnη andnζ areequaltoone.
dilatationalpartofthepositionvectorgradientiseval-
uatedattheelementcenteronly[39];however,itwas
4.3 F-barmethod appliedtolow-orderelementsanddoesnotprovidesat-
isfactoryresultsforfullyparameterizedANCFbeams.
Themethodofselectivereducedintegrationissuitable AssuggestedbyElguedjetal.[8],theL2projectionof
for the problems where it is easy to split the energy strains is appropriate in the case of the higher-order
density function into deviatoric and volumetric parts. elements. Moreover, the associated space on which
In the classical finite element method, the strain pro- theprojectionisperformedshouldhavealowerorder
jectionB-barapproachistreatedasanextensionofthe and then the original element approximation given in
reduced and selective integration techniques [17]. In Sect.2.Itturnsoutthattheproperspaceforthestan-
turn,thegeneralizationofB-barmethodtofinitestrain dardfullyparameterizedspatialbeamelementshould
formulationisso-called F-barmethod,introducedby beconstantinthetransversaldirection.Therefore,one
deSouzaNetoetal.[39],whichinvolvestheproduct candefinethelower-orderbasisusingtildenotationas:
(cid:2) (cid:3)
decompositionofthepositionvectorgradientintovol- S˜ = 1−3ξ2+2ξ3,l(ξ−2ξ2+ξ3),3ξ2−2ξ3,l(−ξ2+ξ3)
umetricanddeviatoricparts[6,8].
(27)
The matrix F may be multiplicative split into its
deviatoric(volumepreserving)andvolumetric-dilata-
and the Eq. (25) may be written in the new space as
tionalparts:
[8,17]:
F = FdilFdev, Fdev = J−1/3F, Fdil = J1/3I J1/3 = S˜J˜1/3 (28)
(22)
where:
(cid:4)
where I isidentitymatrix.Individualpartsofthissplit
J˜1/3 = M˜−1 S˜TJ1/3dV (29)
havethefollowingproperties:
V
det(F)=det(Fdil)= J, det(Fdev)=1 (23) ˜
where M isthe“mass”matrixofthetildebasis:
(cid:4)
Nowwecandefineamodifiedmatrix Fintermsof M˜ = S˜TS˜dV (30)
thedeviatoricpartandthemodifieddilatationalpartas:
V
F¯ = F¯dilFdev (24) ThepresentedprocedurecorrespondstoL2projection
where F¯dil = J1/3I and J1/3maybegivenby: oftheJ1/3intotheS˜andcanbeusedwithbothelement
types affected by the volumetric locking phenomena,
J1/3 =π(J1/3) (25) namelyBE24andBE30S.
where π is the linear projection operator described
4.3.2 Implementationdetails
below.Finally,onecanwrite:
F¯ = J1/3F (26) ThemodifiedpositionvectorgradientF¯ shouldbeused
J1/3 toevaluatethestrainenergydensityfunction:
(cid:4) (cid:4)
¯
Thevalueofthe Fisdirectlyusedintheexpression
UF-bar = U dV + A U¯ dl (31)
for the volumetric penalty function. There is no need sic s p
V l
123
458 G.Orzechowski,J.Fra˛czek
whereindexF-bardenotesthestrainprojectionmethod
and the modified volumetric energy penalty function
¯
U may be calculated by substitution of the Eq. (26)
p
into(13)andusingtheproperties(23)as: Fig.3 Veryflexiblephysicalpendulum
(cid:8)(cid:9) (cid:10) (cid:11)
U¯ = 1k J1/3 3−1 2 (32)
p 0.2
2
Finally:
(cid:12)(cid:5) (cid:4) (cid:6) (cid:13) 0
U¯ = 1k S˜M˜−1 S˜TJ1/3dV 3−1 2 (33) m] -0.2
p 2 V n[
o
Aftertheevaluationofthe J˜1/3vector,thevalueof ositi -0.4
p
the J1/3 dependsonlyonthelongitudinalcoordinate; p
¯ ti -0.6
therefore, in Eq. (31), theUp is evaluated only along al
c
thebeamcenterline. erti -0.8
V
5 Numericalexamples -1
Presented methods of preventing volumetric locking -1.2
0 0.5 1 1.5 2
are applied to the fully parameterized spatial ANCF Time [s]
beam elements BE24 and BE30S. Elements BE30T
and BE42 include the trapezoidal deformation mode Fig.4 Verticalpositionofthebeamtip.( )Compressible
Neo-Hookean,( )IncompressibleNeo-HookeanwithF-bar
of the cross section that suppresses the occurrence of
projection, ( ) Incompressible Mooney–Rivlin with F-bar
thevolumetriclockingintheseelements;therefore,the projection
methods from the previous section are not applied to
them.Inthefollowingexamples,theefficiencyofthe
absolutenodalcoordinateformulationformodelingthe Thependulumatoneendisconnectedtotheground
nearlyincompressiblematerialmodelsisstudiedusing by a spherical join, as shown in Fig. 3. The base
threesimplemodelsofhighlydeformablebeams. of the pendulum is subjected to prescribed motion
u = −0.02sin(2πt) m, where u denotes the dis-
X X
5.1 Physicalpendulum placementsofthebaseinthe X direction,andt isthe
timeexpressedinseconds.Thatadditionalbasemotion
Aflexiblependulumundergravityforcesisemployed increasesbodydeformations.
for the purpose of numerical verification of the vol- Pendulum body is built with four fully parameter-
umetric elimination methods and to shed some light ized ANCF BE24 elements. It should be pointed out
onissuesassociatedwithelementsBE24andBE30S. that the behavior of the element BE30S is the same
Asimilarpendulumwasexaminedin[23].Intheunde- like the BE24; therefore, for clarity only, the BE24
formedstate,beamhasalengthof1mandauniform is presented in this section. However, the results and
squarecrosssectionofdimension20mm.Thedensity conclusionsarevalidforbothelementtypes.Theper-
of the pendulum is assumed to be 7200 kg/m2. The formed dynamical simulations’ results are shown in
materialmodelishighlyflexible,toallowlargedefor- Figs. 4 and 5. In the current example, the following
mations of the body. In this example, the modulus of materialmodelsareused:compressibleNeo-Hookean
elasticityisequalto6MPaandPoissonrationisequal [34], incompressible Neo-Hookean with and without
to 0.3 (for compressible material model). In the case F-barprojectionmethodandincompressibleMooney–
oftheMooney–Rivlinmaterialmodel,thevaluesofits Rivlinwithandwithout F-barprojection.Inaddition,
coefficients are μ = 0.8MPa and μ = 0.2MPa. topresentthelockinginfluencemoreclearly,theresults
10 01
For the incompressible materials, the penalty method forthemodelswithincompressiblematerialswithout
is used to ensure incompressibility with penalty term locking prevention methods are compared with result
assumedtobek =109Pa. oftherigidpendulum.Onceagain,theresultsarepre-
123
Nearlyincompressiblenonlinearmaterialmodels 459
0.2 1.05
1.04
0
1.03
m] -0.2
[
osition -0.4 F) [–] 11..0012
pp et (
ti -0.6 d 1
al
c
erti -0.8 0.99
V
0.98
-1
0.97
0 0.5 1 1.5 2
-1.2
0 0.5 1 1.5 2 Time [s]
Time [s]
Fig.6 Valueofthedet(F)intime.( )CompressibleNeo-
Fig. 5 Vertical position of the beam tip. ( ) Incompress- Hookean,( )IncompressibleNeo-HookeanwithF-barpro-
ibleNeo-HookeanwithF-barprojection,( )Incompressible jection,( )IncompressibleNeo-Hookean
Neo-Hookean,( )Rigid
1.05
sented only for the F-bar method for clarity and for 1.04
the numerical verification. The behavior of the selec-
1.03
tivereduced integrationisanalogous inthefollowing
examplestothe F-barprojection,andtheconclusions –] 1.02
arevalidforbothlockingeliminationtechniques. F) [ 1.01
Figures 4 and 5 show the pendulum’s free end-tip et (
displacements. All the results in Fig. 4 show large d 1
deformationsofthependulum.Theresultsforthecom- 0.99
pressibleandincompressiblemodelsareverysimilar,
0.98
and they are barely distinguishable at this resolution.
This indicates influence of the volumetric locking in 0.97
-10 -5 0 5 10
incompressiblemodelshasbeensuppressed,sinceno
Y [mm]
over-stiffbehaviormaybenoticedduringbending,and
thatmodelsgivereasonableresults. Fig.7 Valueofthedet(F)alongbeamthickness.( )Com-
The plot in Fig. 5 shows comparisons between pressible Neo-Hookean, ( ) Incompressible Neo-Hookean
withF-barprojection,( )IncompressibleNeo-Hookean
resultsoftheincompressiblemodelswithlockingpre-
ventionbymeansoftheF-barmethod,incompressible
modelswithoutanylockingeliminationtechniqueand
therigidmodel.Itshouldbepointedoutthattheresults onthebody(point AinFig.3)intime.Inturn,Fig.7
oftheincompressiblemodelswithoutF-barprojection showsthevalueofthe J acrossbeamthickness—from
areveryclosetotherigidbodysolution.Atthesame point A to B (see Fig. 3) for given time (t = 0.2 s).
time,differenceinresultsofmodelswithandwithout Onecannoticethatmaterialisfullyincompressiblein
F-barprojectionappliedissignificant.Oneshouldcon- wholevolumeonlyfortheoriginalformulationwith-
cludevolumetriclockinginfluenceonincompressible outalockingeliminationtechnique.Forthecompress-
hyperelasticmaterialmodelresultsisenormous. ible material models and the incompressible material
Figures 6 and 7 show the value of the determinant with F-bar projection, the change of the value of J
ofthedeformationgradienttensor J =det(F),which is noticeable. One can also note that for incompress-
indicatesthecompressibility,asafunctionoftimeand ible material with F-bar projection, J is distributed
position.Figure6showsthevalueof J atafixedpoint symmetricallythroughbeamthicknessasexpected(the
123
460 G.Orzechowski,J.Fra˛czek
0.2
0.15
Fig.8 Veryflexibleclampedbeamundergravityforces m]
[
n 0.1
o
siti
o
projection lower-order basis is considered only at the p 0.05
p
beamcenterline). ti
n
i 0
or
Err
-0.05
5.2 Staticsimulationoftheclampedbeam
-0.1
Inthisexample,theclampedbeamundergravityforces 0 10 20 30 40 50 60 70 80
is analyzed, as shown in Fig. 8. The flexible body Number of elements [–]
has the same material properties and dimensions like
Fig. 9 Convergence analysis—Neo-Hookean material model.
in previous example, and only the cross section size ( )BE42,( )BE30T,( )BE30SSRI,( )BE30SF-
is changed—the height and width of the beam are bar,( )BE24SRI,( )BE24F-bar
increasedto40mm.Pendulumisrigidlyattachedtothe
groundwithoutanyprescribedmotion.Theonlyexter- 0.2
nal force that is acting on the system is gravitational
force.Inthisexample,theresultsoftheANCFmodels 0.15
arecomparedwiththesimulationperformedwithcom- m]
[
mercialfiniteelementpackage[2].Itworthtonotethat on 0.1
thecontinuitybetweenelementsisimposedattheposi- siti
o
tionandfirst-ordergradientlevel(thatisC1continuity p 0.05
p
is imposed). Moreover, the clamped condition is also ti
n
i 0
imposedatthosecoordinates.Foramoredetaileddis- or
cussiononthecontinuitycondition,see,e.g.,[15,28]. Err
-0.05
Themainpurposeofthisexampleistocomparethe
resultsofdifferentANCFbeamelementsandlocking
-0.1
elimination techniques presented in this paper. Simu- 0 10 20 30 40 50 60 70 80
lationisperformedwithbothincompressiblematerial Number of elements [–]
models: Neo-Hookean and Mooney–Rivlin. In addi-
Fig. 10 Convergence analysis—Mooney–Rivlin material
tion,theresultsforbothlockingeliminationmethods,
model. ( ) BE42, ( ) BE30T, ( ) BE30S SRI, ( )
theselectivereducedintegrationandF-barprojection, BE30SF-bar,( )BE24SRI,( )BE24F-bar
arepresented.ThesolutionobtainedwithFEAanaly-
sis is treated as the reference. Therefore, in the FEA
package, a very fine mesh of 1600 higher-order solid Resultsforbothmaterialmodelshavethesamechar-
elements(SOLID186)isemployed. acteristics. Only models with BE42 converge to the
Figures9and10showresultsofthestaticanalysis. referencesolutionwithacceptableerror,andthemodel
Theerrorintippositionsismeasuredbysign(d )|d|, with40BE42elementspredictsthesolutioncorrectly.
z
where |•| denotes quadratic norm, d is a z compo- However, the convergence rate of this element to the
z
nent of the vector d, while d = u − u where u steady response is the slowest, which may be caused
r A r
isreferencetipdisplacementand u istheANCFtip bytheshearlocking[28].ElementBE30Tconvergesto
A
displacement.SolutionsforelementsBE24andBE30S wrongover-stiffsolutionwithnoticeableerror,despite
withoutlockingeliminationtechniquesarenotshown that this element includes a trapezoidal displacement
astheyproduceanenormouserror—irrespectiveofthe mode.Thisisprobablycausebythefactthatwhenvery
element type, material model or number of elements, largedeformationsoccurduetobending,thecrosssec-
theerrorintippositionhasavalueofabout1.19. tiontakestheformoftrapezoidbutwithcurvedsides,
123