Composite Structures 119 (2015) 799–814

Contents lists available at ScienceDirect

Composite Structures
journal homepage: www.elsevier.com/locate/compstruct

A geometrically nonlinear FEM formulation for the analysis of fiber

reinforced laminated plates and shells
M.S.M. Sampaio ⇑, R.R. Paccola, H.B. Coda
São Carlos School of Engineering, University of São Paulo, Brazil

a r t i c l e i n f o a b s t r a c t

Article history: An alternative finite element (FE) formulation for the analysis of fiber reinforced laminated plates and
Available online 18 September 2014 shells developing small and large deformations is presented. Kinematic relations that ensure the adher-
ence of fibers nodes to the continuum medium are used to introduce long and short fibers immersed in
Keywords: any positions of the laminate. The proposed formulation does not increase the number of degrees of free-
Fiber reinforced laminated plates and shells dom of the analyzed shell and do not require the necessity of matching nodes in the discretization of
Long fibers fibers and matrix. Following this procedure, it is possible to model spatially oriented or random fibers
Random short fibers
at any volume rate simplifying the reproduction of real problems. Some understanding examples are
Geometrically nonlinear analysis
Finite element method
tested to demonstrate the good behavior of the proposed formulation and new applications that show
its potential are presented.
Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction increase the number of degrees of freedom of the problem

[25,29]. Moreover, in this case, the reinforcements present the
To well design and satisfy the growing demands of high rotation degree of freedom and bending stiffness that are not
strength-to-weight and high stiffness-to-weight relations of sev- required for fiber reinforcement analysis, the subject of the present
eral engineering materials and structures is necessary to employ formulation.
methods that allow a proper description of fiber reinforced lami- In general the above mentioned formulations exhibit one or
nated medium and structural nonlinear geometric analysis. There both of the following features, that is, increase the number of
are, at least, four different ways to model stiffened composite lam- degrees of freedom of the problem or present the very restrictive
inates available in literature: (i) models in which the laminated requirement that the fiber finite element nodes coincide with the
heterogeneous medium is replaced by an equivalent homogeneous shell finite element ones. Moreover, in some cases, these models
one [1–18]; (ii) models in which the reinforced problem is fully are restricted to geometric linear analysis.
discretized by three-dimensional solid finite elements [19]; (iii) In this sense, the aim of this paper is to develop a FE formulation
models in which the matrix is discretized by plate or shell finite for inclusion of long and short fibers in any position of laminated
elements and the stiffeners are discretized by beam elements shells developing small and large deformations without increasing
[20–28]; and (iv) models in which the fiber stiffness is assembled the number of degrees of freedom. Moreover, in this formulation
in the stiffness of the plate or shell finite element [29–31]. there is not necessity of matching fiber nodes with shell nodes in
An extensive literature review of different composite finite ele- the discretization. The triangular laminated shell finite element used
ment formulations can be seen in [16,25,29]. In short, methods that to discretize the matrix has ten nodes and seven degrees of freedom
replace the original heterogeneous medium by an equivalent per node, that are, three translations, three components of a general-
homogeneous one are useful when the aim is the global system ized vector and the linear rate of strain variation along the thickness
behavior, but make it difficult to identify the fiber–matrix contact [33,34]. The curved fibers, long or short, are introduced in any layer
stresses necessary for some failure analysis, see [16,32] for exam- of the laminate shell by means of kinematic relations to ensure its
ple. Methods that adopt the same element type to approximate adherence to the laminated shell without introducing new degrees
the stiffeners and the continuum present the disadvantage of of freedom in the resulting system of equations.
One-dimensional finite elements of any polynomial order with
three degrees of freedom per node are used to discretize fibers.
⇑ Corresponding author.
These elements are developed to consider geometrical nonlinearity
E-mail addresses: [email protected] (M.S.M. Sampaio), [email protected]
(R.R. Paccola), [email protected] (H.B. Coda).
by mean of a deformation function based on positional mapping.

http://dx.doi.org/10.1016/j.compstruct.2014.09.009
0263-8223/Ó 2014 Elsevier Ltd. All rights reserved.
800 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

The Total Lagrangian description of motion is adopted and the As the current position or the equilibrium position is the
Saint–Venant–Kirchhoff Constitutive Law that relates linearly the unknown of the problem, when a trial solution (Y 0l ) is tested, a
second Piola–Kirchhoff stress tensor and the Green–Lagrange non-null vector gk appears. In this situation gk, understood as the
strain tensor is used to model the material. The equilibrium is unbalanced force of the mechanical system, is non-linear regarding
achieved from the Principle of Minimum Potential Energy and nodal parameters and indicates that a correction in the trial posi-
the Newton–Raphson iterative procedure is used to solve the non- tion should by applied in order to achieve equilibrium. To find this
linear system of equations. correction it is necessary to expand gk, as follows:
Compared to elements based on classical beam, plate and shell   @g 
j
g j ðY l Þ ¼ g j Y 0l þ DY k þ O2j ¼ 0 ð1:4Þ
theories the present shell element uses a complete elastic constitu- @Y k ðY 0 Þ
tive relation and does not use nodal rotations as parameters. Fur- l

thermore, the good performance of the element relative to Neglecting higher order terms O2j in Eq. (1.4) results:
volumetric and shear locking phenomena is ensured by the intro-  !1
@g j   
duction of the linear rate of thickness variation [33,34]. DY k ¼   g j Y 0l ð1:5Þ

@Y k ðY 0 Þ
To present the proposed formulation the paper is organized as  
l
@g j  2 
follows: Section 2 presents the geometrically nonlinear equilib- where ¼ @W 
@Y k 
is the complete Hessian matrix of the
rium problem; in Section 3 the laminated shell finite element used ðY 0l Þ @Y k @Y j ðY 0l Þ
total potential energy.
to model laminated orthotropic plates or shells is described; Sec-
tion 4 presents the developed fiber finite element including its From Eq. (1.5) one achieves the full Newton–Raphson procedure
kinematics, internal force vector and Hessian matrix achieve- to solve nonlinear system of equations, summarized as follows:
ments; Section 5 introduces the spreading strategy that makes One chooses a trial position Y 0l and calculates the unbalanced force
 
possible a complete analysis of reinforced laminated plate or shell vector g j Y 0l according to Eq. (1.3). Using Eq. (1.5) one finds the
without increasing the number of degrees of freedom; Section 6
position correction DYk applied over the trial solution Y 0l as
presents numerical examples validating the proposed formulation
and, finally, conclusions are presented in Section 7. Y l ¼ Y 0l þ DY l . With this new position vector, one repeats the pro-
 
Index notation is used in most sections; however in some parts cedure until DYl or g j Y 0l become smaller than assumed
of the text dyadic notation is employed to facilitate understanding. tolerances.
_
The upper symbols ð Þ and ðÞ~ are used to distinguish, respectively, The first and the second derivatives of the strain energy regard-
shell and fiber variables. ing shell nodal positions, including matrix and fiber contributions,
are performed and presented in next sections.
2. The fiber reinforced laminated shell nonlinear equilibrium It is important to mention that the strain energy per unit of ini-
problem _
tial volume used in this work, for the shell element w and for the
fibers element w, ~ describes a Saint–Venant–Kirchhoff material,
The geometrically nonlinear analysis starts considering that the
that is, its constitutive Law is given by the linear relationship
change configuration process is governed by the First Law of Ther- _
modynamics. The total potential energy P generated by the sys- between the second Piola–Kirchhoff stress tensor S ij and the
_
tem, ignoring the contribution of the kinetic energy and Green–Lagrange strain tensor E ij . The expressions for
considering adiabatic process, is given by: _ _ _
~ S ij and E ij are given, respectively, by:
w; w;
P¼W þT ð1:1Þ
_ 1_ _ _
where T =  FkYk is the potential energy of external conservative w¼ E ij Cijkl E kl ð1:6Þ
_
~ is the 2
forces Fk applied in the current positions Yk and W ¼ W þW 1 ~ ~2
strain energy stored in the body including the matrix strain energy ~ ¼ E
w E ð1:7Þ
_
~ As all involved quantities are Lagrang- 2
W and fiber strain energy W. _
ian, written regarding laminated shell nodal positions, the whole
_ @w _ _
S ij ¼ _ ¼ Cijkl E kl ð1:8Þ
strain energy stored in the initial volume of the body is given by: @ E ij
Z   _ Z  _     
_ _ _ 1 _ 1 _ _
W¼ _
w Y li dV 0 þ ~ Y~ pj Y li dV
w ~0 ð1:2Þ E ij ¼ C ij  dij ¼ Aki Akj  dij ð1:9Þ
V0 V~ 0 2 2
_ _ _ _
in which w and w ~ are, _
respectively, shell and fiber strain energy per ~ C ij ; Aki and dij are, respectively, the linear elastic con-
~ E;
where Cijkl ; E;
unit of initial volume; Y li are shell nodal positions; Y ~ p are fiber nodal stitutive fourth-order tensor, the one-dimensional elastic modulus,
_ j
positions; V 0 is the initial volume of shell finite elements and V ~ 0 is the one-dimensional Green strain, the right Cauchy-Green stretch
the initial volume of fiber finite elements. tensor, the deformation gradient and the Identity tensor.
In Eq. (1.2) index l = 1, . . . , 10 and i = 1, . . . , 7 represent, respec-
tively, the shell finite element nodes and the shell degrees of free- 3. The kinematic description and the strain energy function of a
dom per node and index p = 1, . . . ,N and j = 1, 2, 3 represent, laminated shell
respectively, the fiber finite element nodes and fiber degrees of
_
freedom per node. Let x m
i be a point in the reference surface of the initial configu-
_
As the external forces are conservative, the Principle of Mini- ration X0 of the shell and g 0i a generic vector related to this point as
_
mum Total Potential Energy is applied in Eq. (1.1) and the equilib- shown in Fig. 1. In the same way, let y m i be a point
_ in the reference
rium equation can be written as: surface of the current configuration X and g 1i a generic vector
Q _
@ @W related to this point. Then, a generic point x i at the initial configu-
_
gk ¼ ¼  Fk ¼ 0 ð1:3Þ ration and a generic point y i at the current configuration can be
@Y k @Y k
@W
written, respectively, as:
where @Y is the gradient vector of the strain energy understood as
k _ _ _
the internal force of the system. From Eq. (1.3) one observes that xi ¼ xm 0
i þ gi ð1:10Þ
the solid is at equilibrium when the internal force vector is equal _ _m _1
to the applied one. yi ¼ yi þ gi ð1:11Þ
 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814 801

Fig. 1. Position vectors for kinematic description.

_ _ _
The points x m
i in the reference surface of the initial configura-
in which N0il and Gik are, respectively, nodal values of initial and cur-
_
tion are mapped from a dimensionless space X1 with coordinates rent generalized vectors and a l are nodal values of the linear rate of
ni using any order shape functions u ^ l ðn1 ; n2 Þ and the coordinates thickness variation [33,34].
_
^ l of nodes l in the initial configuration. The same is true for the
X From Eqs. (1.10)–(1.15) the mapping u0lam from auxiliary con-
i i _
current configuration. Then, we have, respectively, figuration to the initial configuration and the mapping u1lam i from
_m _ _ auxiliary configuration to the current configuration can be written,
x i ¼ /l ðn1 ; n2 ÞX li ð1:12Þ respectively, as:
_ _ _ 0 _ 1
ym l
i ¼ /l ðn1 ; n2 ÞY i ð1:13Þ _ 0lam _ _ _ _ h lam _ _
_ ui ¼ x i ¼ /l ðn1 ; n2 ÞX li þ @ lam
d þ 0
n A/l ðn1 ; n2 ÞN0il ð1:16Þ
Let d lam be the distance from the reference surface to the center 2 3
_
of a considered_lamina following direction g 0i in the initial config-
lam
uration and let h 0 be the initial thickness of this lamina
_ indicated and
_
by the superscript lam, see Fig. 2. The vectors g 0i and g 1i can be writ-
_ _ _ _
ten, respectively, as: u1lam ¼ y i ¼ /l ðn1 ; n2 ÞY li
i
0 1 20 1 0 12 3
_ _ _
_ _
@ lam h lam _ _ lam
6@_lam h 0 _ _ _ h lam
7
g 0i ¼ d þ 0
n3 A/l ðn1 ; n2 ÞN0il ð1:14Þ þ4 d þ n A þ /l ðn1 ; n2 Þa l @ d þ 0 n3 A 5
lam
2 2 3 2
20 1 0 12 3
_ _ _ _
_
6 _ h lam _ _ _ h lam 7  /k ðn1 ; n2 ÞGik ð1:17Þ
g 1i ¼ 4@ d lam þ 0
n3 A þ /l ðn1 ; n2 Þa l @ d lam þ 0
n3 A 5
2 2
_
The linear rate of thickness variation, a l , introduced in Eqs.
_ _
 /k ðn1 ; n2 ÞGik (1.15) and (1.17) avoid volumetric and shear locking from usual
homogeneous or laminated plate or shell formulations, see
ð1:15Þ
[33,34] for more details.
_
Following Fig. 1, function u0lam
i is used to find the initial gradi-
_ _
ent of the auxiliary mapping A0lam
ij while function u1lam
i is used to
_
find the current gradient of the auxiliary mapping A1lam
ij . The com-
position of these two values for each integration point gives the
_
numerical value of the deformation gradient Alam , i.e.:
_ _ _ 1
_ @ u0lam _ @ u1lam _ _
A0lam
ij ¼ i
; A1lam
ij ¼ i
; Alam ¼ A1lam  A0lam
@nj @nj
ð1:18Þ
_ _
The mapping gradient A1lam
or the derivatives of mapping u1lam
ij i
regarding the non-dimensional variables, that is, from the auxiliary
Fig. 2. Kinematic description for laminated shells. configuration to the current configuration are given as:
802 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

0 _ 12
lam
_ _ _ _ _ _ h _ _
A1lam
i1 ¼ /l;1 Y li þ /l;1 a l @ d lam þ 0 n3 A /k Gik
2
20 1 0 12 3
_ _
lam lam
6 _ h __ _ h 7_ _
þ 4@ d lam þ 0 n3 A þ /l a l @ d lam þ 0 n3 A 5/k;2 Gik ð1:19Þ
2 2
0 _ 12
_ _ _ _ _ _
@ lam h lam _ _
A1lam
i2 ¼ /l;2 Y li þ /l;2 a l d þ 0
n3 A /k Gik
2
20 1 0 12 3
_ _
lam lam
6 _ h _ _ _ h 7_ _
þ 4@ d lam þ 0 n3 A þ /l a l @ d lam þ 0 n3 A 5/k;2 Gik ð1:20Þ
2 2
_   _  _
_ h0 __ _ _
A1lam
i3 ¼ 1 þ /l a l 2d lam þ h lam
0 n3 /k Gik ð1:21Þ
2 Fig. 3. Fiber finite element mapping for initial and current configurations.
_
The expressions of A0lam
ij are left to the reader.
Remembering that the strain energy per unit of initial volume of
the shell finite element is given by Eq. (1.6) the first derivative of mapped from the non-dimensional auxiliary configuration X1,
the strain energy regarding nodal shell positions is given by: respectively, by the expressions:
_ _ _ _ _ _
@w @w @E @C 1 @w @C ~ 0i ¼ /
~xi ¼ u ~p
~ p ðgÞX ~ p ðgÞY~ p
~ 1i ¼ /
~i ¼ u
and y ð1:26Þ
_ ¼ _ : _ : _ ¼ _ : _ ð1:22Þ i i

@Y a b
@ E @ C @Y a 2 @ E @Y ab
b
where ~p
Xare coordinates of fiber nodes p in the initial configuration,
i
_ _ ~ p are coordinates of fiber nodes in the current configuration and
in which S ¼ @w
_ is the second Piola-Kirchhoff stress tensor, see Eq. Y i
@E ~ p ðgÞ are shape functions of any order evaluated in coordinates g.
/
(1.8) The first derivative of the Cauchy-Green stretch tensor regard-
In Eq. (1.26) index p = 1, . . . , N and i = 1, 2, 3 represent, respec-
ing nodal positions in Eq. (1.22) is:
tively, the fiber finite element nodes and fiber coordinate
_ _ T _ T  1 _ T _ T _1 _ 1 directions.
@C 0 @ðA1 Þ _1 _0 @A The tangent vector of the fiber and its modulus are calculated at
_ ¼ A  _ A  A þ A0  A1  _  A0
b b
@Y a @Y a @Y ba the initial configuration, respectively, as:
ð1:23Þ 0 ~ p ðgÞ p
d/
T Xi ¼ ~
X and
i
In the same way, the second derivative of the strain energy dg
regarding nodal shell positions is given by: !2 !2 !2
 0 2 ~ p ðgÞ p
d/ ~ p ðgÞ p
d/ ~ p ðgÞ p
d/
~X  ~ ~ ~
2 _ 2 _
_ _ _ _ T  ¼ X þ X þ X ð1:27Þ
@ w 1 @ w @C @C 1 @w @2 C dg 1
dg 2
dg 3

_ _ ¼ _ _ : _ : _ þ _ : _ _ ð1:24Þ
@Y ba @Y nc 4 @ E @ E @Y n @Y b 2 @ E @Y n @Y b
c a c a For the current configuration one finds:
and the second derivative of the Cauchy-Green stretch tensor ~ p ðgÞ p
d/
regarding nodal position results: T Xi ¼ Y~ i and
dg
_ _ T 2 _1 T _ _ 1 !2 !2 !2
@2 C @ ðA Þ  2 ~ p ðgÞ p
d/ ~ p ðgÞ p
d/ ~ p ðgÞ p
d/
A0  _ _  A1  A0 ~X 
_ _ ¼ T  ¼ Y~ 1 þ Y~ 2 þ Y~ 3 ð1:28Þ
b n
@Y a @Y c @Y ba @Y nc dg dg dg
_ T
_ T @ A1 _  1 From Eqs. (1.27) and (1.28), the one-dimensional Green strain is
0 @A1 _0 written as:
þ A  _  _  A
@Y ba @Y nc 0 2  2 1
~X   X0 
_ T T   ~ T 
~ ¼ 1B
E @  2
C
A ð1:29Þ
_ T @ A1 _  1 2 ~X0 
@A1 _ T 
þ A0  _  _  A0
@Y cn @Y ab
In order to proceed with the equilibrium analysis it is necessary
_ T _ T _ _ 1 to know the first derivative of strain energy regarding positions.
@ 2 A1
þ A0  A1  _ _  A0 ð1:25Þ For the one-dimensional fiber we use the Saint–Venant–Kirchhoff
@Y ba @Y cn constitutive law given by Eq. (1.7). From the energy conjugate con-
~
cept the natural internal fiber force vector, @W
~j
, related to node j and
@Y
k
4. The fiber kinematics and strain energy function direction k is calculated regarding fiber parameters as:
~ p ðgÞ p
 ~ ðgÞ
Any order curved finite elements with three degrees of freedom ~ Z d/
Yk
d/ j
@W ~E~ dg dg
~0
per node are used to discretize fibers. These elements are devel- ¼ E  2 dV ð1:30Þ
oped to consider geometrical nonlinearity using Saint–Venant– @Y j
~
k
~
V 0
~X0 
T 
Kirchhoff constitutive Law and the kinematics description shown
in Fig. 3. Developing the necessary calculations and considering the ini-
The points with coordinates ~xi in initial configuration X0 and ~ and non-dimen-
tial volume as a function of the transverse area A
~
the points with coordinates yi in the current configuration X are sional coordinate g, one writes the following simple expression:
 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814 803

Fig. 4. Fiber reinforced cantilever beam.

5. Fiber – laminated shell coupling

To be possible embed fibers at any position in the domain (plate

or shell) without increasing the degrees of freedom and without
matching nodes, the positions of fibers nodes are written as func-
tion of the positions of laminated shell nodes as:
_ _
Y~ pi ¼ /l ðnp1 ; np2 ÞY li
20 1 0 12 3
_ _
lam lam
6 _ h _ _ _ h 7_ _
þ 4@ d lam þ 0 np3 A þ /l ðnp1 ; np2 ÞY l4 @ d lam þ 0 np3 A 5/k ðnp1 ; np2 ÞY kðiþ4Þ
2 2

Fig. 5. Cross-sectional geometry: (a) present solution (b) homogenized solution. ð1:36Þ
_
  ~ ðgÞ
in which the symbols ð Þ and ðÞ
~ distinguish, respectively, shell and
~ p ðgÞ p
~ Z ~L d/
Y~ k
d/j
fiber variables and p indicates a fiber node that belongs to a shell
@W 0
dg dg
~ ~s0
F~jk ¼ j ¼ ~E
E ~  2 Ad _
element. In Eq. (1.36) the rate of thickness variation a l of node l is
~
@Y 0 ~X0 
k T  defined as the fourth degree of freedom of the shell node k and
_
~ p ðgÞ p
 ~ ðgÞ the generalized vector Gki components of node k are defined as the
Z 1
d/
Y~ k
d/j

~E~ dg dg
~ g
~J 0 ðgÞAd other three shell degrees of freedom. At this point we consider
¼ E  2 ð1:31Þ
1 ~X0  the lamina and the non-dimensional shell coordinates npk of each
T 
fiber node known.
Remembering that the strain energy stored in a laminated rein-
In the same way, the Hessian matrix components for the fiber
forced body is the sum of the strain energy stored in the matrix and
element are obtained by the second derivative of the strain energy,
fiber, see Eq. (1.2), the internal force at a shell node b following
i.e.:
Z direction a, considering both fiber and matrix contributions, is
@2W ~ @2w ~
Hjb dV~ 0 found by the conjugate energy concept as:
ka ¼ ~ j ~b
¼
~ j ~b
ð1:32Þ
@Y k@Y a V~ 0 @ Y k @ Y a Z _ Z  _ 
@W @w _ ~
@w
Developing the necessary calculations, Eq. (1.32) results in: _ ¼ _ _ dV 0 þ _ Y~ pj Y li dV~ 0 ð1:37Þ
0 1 @Y ba V0 @Y ba ~
V 0 @Y ba
Z ! !
2 ~ ~ ~ ~ ~ ~
@ W B E d/p ðgÞ ~ p d/b ðgÞ d/p ðgÞ ~ p d/j ðgÞC ~ The first integral at the last term of Eq. (1.37) is known from Eq.
¼ @  4 Ya Yk AdV 0
@ Y~ jk @ Y~ ba V~ 0 ~X0  dg dg dg dg (1.22). However, it is necessary to differentiate the fiber strain
T 
0 1 energy regarding the matrix nodal coordinates using the chain
Z ~E~ d/~ b ðgÞ d/
~ j ðgÞ C rule, that is:
B E
þ @  2 dka AdV~ 0
~
V ~X0  dg dg @w~ @w~ @ Y~ p @ Y~ p
0
T  F~ba ¼ _ ¼ p _i ¼ _i F~pi ð1:38Þ
~
@Y ba @ Y i @Y ba @Y ba
ð1:33Þ
where F~pi ¼ @@Y~w~p is the nodal internal force of the fiber given in Eq.
or ~ p i
@Y
0
! !
1 (1.31) and _
b
i
is achieved deriving Eq. (1.36). It is obvious that these
~ Z 1 ~ p ðgÞ ~ b ðgÞ d/ ~ p ðgÞ p ~ j ðgÞC
@Y a
@2W B E ~ d/ d/ d/ ~ g
¼ @  4 Y~ pa Y~ k A~J 0 ðgÞAd derivatives are not null only if the represented directions coincide.
~ j ~b
@Y k @Y a   d g d g dg dg
~
0
1
TX  For each shell node b and a = 1, 2, 3 results a diagonal 3  3
Z þ1 ~ ~ ~ b ðgÞ d/ ~ j ðgÞ matrix given by:
EE d/ ~ g
þ  2 dka~J 0 ðgÞAd _
1 ~X0  dg dg
T  @ Y~ pi _ @Y il _ _
wbai ¼ _ ¼ /l ðnp1 ; np2 Þ _ ¼ /l ðnp1 ; np2 Þdbl dai ¼ /b ðnp1 ; np2 Þdai ð1:39Þ
ð1:34Þ @Y ba @Y ab
where For each shell node b and a = 4 the following three position vec-
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tor results:
 0  2  2  2
~J 0 ðgÞ ¼ ~X  d~x1 d~x2 d~x3 2 0 12 3
T ¼ þ þ ð1:35Þ _
dg dg dg @ Y~ pi lam
6_ p p @_lam h 0 p A 7_ p p _k
wb4i ¼ _ ¼ 4/b ðn1 ; n2 Þ d þ n 5/k ðn1 ; n2 ÞY iþ4 ð1:40Þ
~ 0i , see Fig. 3.
is the differential Jacobian of u @Y b4 2 3
Eqs. (1.31) and (1.34) are solved using the Gauss–Legendre
quadrature. Finally, for each shell node b and a = 5, 6, 7 one achieves:
804 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

(a)

(b)

(c)

(d)

Fig. 6. Final configuration of the beam: (a) 3D frame solution in (cm); (b) 2D solid solution in (cm); (c) homogenized solution in (cm); (d) present solution in (cm).

5
Timoshenko-Reissner 3D frame solution
2D solid solution
0 Homogenized solution
Present solution

-5
Displacements (cm)

-10

-15

-20

-25

-30
0 30 60 90 120 150 180 210 240 270 300 330
Coordinates X (cm)

Fig. 7. Transverse displacement of the beam (cm).

 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814 805

Table 1
Maximum transverse displacement – numerical values in centimeters.

2D solid solution (reference) 3D frame solution Homogenized laminated solution Present solution
23.912 22.757 23.622 23.611
Relative difference 4.83% 1.212% 1.258%

Fig. 8. Pinched problem extract from [39] and adopted discretization.

Table 2
Physical properties used to model the problem with the homogenized formulation: (E)  106 and (G)  106.

Ply arrangement Lamina Matrix

E m G
E1 E2 E3 m12 m13 m23 G12 G13 G23
0°/0°/0° 1 40 1 1 0.25 0.25 0 0.6 0.6 0.5
2 40 1 1 0.25 0.25 0 0.6 0.6 0.5
3 40 1 1 0.25 0.25 0 0.6 0.6 0.5
90°/90°/90° 1 1 40 1 0.25 0 0.25 0.6 0.5 0.6
2 1 40 1 0.25 0 0.25 0.6 0.5 0.6
3 1 40 1 0.25 0 0.25 0.6 0.5 0.6
0°/90°/0° 1 40 1 1 0.25 0.25 0 0.6 0.6 0.5
2 1 40 1 0.25 0 0.25 0.6 0.5 0.6
3 40 1 1 0.25 0.25 0 0.6 0.6 0.5
90°/0°/90° 1 1 40 1 0.25 0 0.25 0.6 0.5 0.6
2 40 1 1 0.25 0.25 0 0.6 0.6 0.5
3 1 40 1 0.25 0 0.25 0.6 0.5 0.6

0º/0º/0º 90º/90º/90º 0º/90º/0º 90º/0º/90º

(a)

(b)
Fig. 9. Fiber–matrix ply arrangements: (a) fibers in midsurface of the considered lamina and (b) two fibers by lamina away ±0.25 from the midsurface of considered lamina.
806 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

Table 3
Physical properties used to model the problem with the present formulation. (E)  106 and (G)  106.

Ply arrangement Lamina Fibers Matrix

Nf Ef Af E m G
E1 E2 E3 m12 m13 m23 G12 G13 G23
0°/0°/0° 1 47 39 1 1 1 1 0.25 0.25 0 0.4 0.4 0.5
2 47 39 1 1 1 1 0.25 0.25 0 0.4 0.4 0.5
3 47 39 1 1 1 1 0.25 0.25 0 0.4 0.4 0.5
90°/90°/90° 1 30 39 1 1 1 1 0.25 0 0.25 0.4 0.5 0.4
2 30 39 1 1 1 1 0.25 0 0.25 0.4 0.5 0.4
3 30 39 1 1 1 1 0.25 0 0.25 0.4 0.5 0.4
0°/90°/0° 1 47 39 1 1 1 1 0.25 0.25 0 0.4 0.4 0.5
2 30 39 1 1 1 1 0.25 0 0.25 0.4 0.5 0.4
3 47 39 1 1 1 1 0.25 0.25 0 0.4 0.4 0.5
90°/0°/90° 1 30 39 1 1 1 1 0.25 0 0.25 0.4 0.5 0.4
2 47 39 1 1 1 1 0.25 0.25 0 0.4 0.4 0.5
3 30 39 1 1 1 1 0.25 0 0.25 0.4 0.5 0.4

@ Y~ pi As a rule, for fiber nodes belonging to the same shell element the
wbai ¼
@Y a
_
b sub-matrices [H]p are placed in the same macro-column, if fiber
20
_ 1 0 _ 12 3 nodes belong to different elements the sub-matrices [H]p are
lam lam
6@_lam h 0 p A _ p p _‘ @_lam h 0 p A 7_ p p placed in different macro-columns.
¼4 d þ n þ /‘ ðn1 ; n2 ÞY 4 d þ n 5/b ðn1 ; n2 Þdða4Þi
2 3 2 3 Proceeding in the same way as described for the calculation of
internal forces, the second derivative of strain energy of the rein-
ð1:41Þ forced finite element regarding the plate or shell nodal parameters
Therefore, for a fiber node p connected to a shell, the total ten- is developed as follows:
sor wbai for a single shell node is of order 7  3 and its transpose can Z _ Z  _ 
@2W @2 w _ @2w
~
be written as: _ _ ¼ _ _ _ dV 0 þ _ _ Y~ p Y l
j i
~0
dV ð1:48Þ
2 3p @Y ab @Y cn V0 @Y ba @Y nc ~
V 0 @Y ba @Y nc
n oT wb11 0 0 wb41 wb51 0 0
6 7 The first integral kernel at the right side of Eq. (1.48) is the Hes-
½K pb ¼4 0 wb22 0 wb42 0 wb62 0 5 ð1:42Þ
sian matrix of the homogeneous shell, known from Eq. (1.24). It is
0 0 wb33 wb43 0 0 wb73 necessary to observe that the kernel of the last integral is the spe-
As the adopted shell finite element has 10 nodes Eq. (1.38) can cific strain energy of a fiber derived twice regarding the plate or
be rewritten, for one fiber node, in a matrix form as: shell nodal parameters. Therefore, applying the chain rule one
_
finds:
F ¼ ½Hp  F~ ðpÞ
p
ð1:43Þ
@2w
~ @2w ~ @ Y~ qx @ Y~ gp
p
in which [H] is a 70  3 matrix and its transpose is given by: _ _ ¼ _ _ dwa dqb dpc dgn ð1:49Þ
@Y ab @Y cn @ Y~ qx @ Y~ gp @Y ba @Y nc
p T p p p p
½H ¼ ½Kp1 ½Kp2 ½K3 ½Kp4 ½ K 5 ½Kp6 ½ K 7 ½Kp8 ½ K 9 ½Kp10 2
in which @Y~@q @w~Y~ g is the Hessian matrix of the fiber element written
ð1:44Þ x p
regarding fiber nodes, see Eq. (1.32) or Eq. (1.34).
When n fiber nodes are considered, if each fiber node belongs to a Following the same strategy used to organize the internal force
different shell element, the transformation from the internal fiber procedure, Eq. (1.49) is rewritten as:
force to the internal equivalent shell force, Eq. (1.38) is written sim- _
~  ½WT
H ¼ ½W  H ð1:50Þ
ply as:
_ To finish the proposed procedure it is necessary to set the lam-
F ¼ ½W  F~ ð1:45Þ ina position and to automatically find the non-dimensional coordi-
in which [W]70n  3n is given by: nates npk of any fiber node inside any shell element.
2 3 As laminated plates or shells are the subjects of this study and
½H170  3 0  0 by the usual manufacturing process of this kind of structural ele-
6 7
6 0 ½H270  3  0 7 ments, the lamina and the third non-dimensional variable np3 for
6 7
½W ¼ 6 .. .. .. 7 ð1:46Þ any fiber are considered previously known, i.e., prescribed at the
6 7
4 . . . 5 generation process.
0 0    ½Hn70  3 70n  3n
The pair of dimensionless plate or shell variables ðnp1 ; np2 Þ associ-
ated with the physical fiber node position is solved in the following
If all fiber nodes belong to the same shell element the following nonlinear system:
matrix [W]70  3n results:
_ _
2 3 ~ p ¼ /l ðnp ; np ÞX l
X ð1:51Þ
½H170  3 i 1 2 i
6 7 _
~ p are the
6 ½H270  3 7 where /l are shape functions of the plate or shell element, X
6 7 i
½W ¼ 6 .. 7 ð1:47Þ known physical fiber nodes coordinates, generated independently
6 7 _
4 . 5 of the plate or shell mesh, and X li are the plate or shell nodes coor-
½Hn70  3 70  3n
dinates. To solve Eq. (1.51), it is expanded into a Taylor series of the
first order around a trial dimensionless coordinate ðnpt pt
1 ; n2 Þ:
 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814 807

0º/0º/0º

(a) (b)
90º/90º/90º

(a) (b)
0º/90º/0º

(a) (b)
90º/0º/90º

(a)
(b)
Fig. 10. Final configuration of 1/8 of the cylinder for small displacements: (a) homogenized solution and (b) present solution.

_  coordinates and Hij is a two-dimensional matrix. The correction of

~ p
_
pt pt
_
l @/l n;1 n2  ~p ¼ X
~ pt þ Hij Dnj
X i ffi /l ðn1 ; n2 ÞX i þ  Dnj or X i i the trial dimensionless coordinates Dni is calculated by solving
@nj  pt
ðn1 ;n2 Þ
pt the following linear system of equations:
ð1:52Þ ~p  X
Hij Dnj ¼ X ~ pt ð1:53Þ
i i

where ~ pt
X is a trial position of the fiber node calculated from the The procedure is a simple Newton–Raphson nonlinear solver
i
plate or shell element geometry and the trial dimensionless that relates all fiber nodes to the connected plate or shell element,
808 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

Table 4
Numerical results for the non-dimensional deflection of the cylinder under the applied load: linear case.

Ply arrangement Homogenized solution Present solution

1 Fiber line per lamina 2 Fiber lines per lamina
EhW/nP Relative difference (%) EhW/nP Relative difference (%)
0°/0°/0° 2875.32 3226.08 12.199 3190.56 10.964
90°/90°/90° 995.628 1120.67 12.559 1077.11 8.184
0°/90°/0° 1908.6 3131.04 64.049 2286.36 19.793
90°/0°/90° 776.124 949.764 23.373 875.58 12.814

0.30

0.25

0.20

0.15
Load

0º0º0º present solution

0º0º0º homogenized solution
90º90º90º present solution
0.10
90º90º90º homogenized solution
0º90º0º present solution
0º90º0º homogenized solution
0.05 90º0º90º present solution
90º0º90º homogenized solution

0.00
0 5 10 15 20 25 30 35 40 45 50 55 60
Displacements

Fig. 11. Equilibrium path of the point under the applied load P.

revealing the pair of dimensionless variables ðnp1 ; np2 Þ. From this frame formulation, the matrix is discretized in 10 fifth order finite
information one also knows the current fiber node positions as a elements, see Fig. 6(b). For laminated shell formulation with
function of the plate or shell nodes positions, see Eq. (1.36), and homogenized orthotropic material five laminas are employed, as
Eqs. (1.50) and (1.45) can be applied. depicted in Fig. 5(b), and 60 triangular shell elements are used to
discretize the problem, as shown in Fig. 6(c). Finally, for the devel-
oped formulation only one lamina is employed to discretize the
6. Numerical examples
matrix (with 60 triangular cubic shell elements) and 30 cubic fiber
elements are used to discretize each reinforcement, see Fig. 6(d). It
Five examples are shown to demonstrate the good behavior of
is worth noting that Fig. 6(c) and (d) are in a perspective that
the proposed formulation. The first one is used to confirm that
allows the view of the immersed fibers.
the mechanical coupling between fibers and the laminated shell
The adopted geometrical properties for the beam are:
finite elements is working properly for large displacement analy-
L = 300 cm, b = 10 cm, h = 20 cm and d = 3 cm. The matrix Young
ses. The second and the third deal with a fiber reinforced laminated
modulus and the Poisson ratio are given, respectively, by Em = 21
cylindrical shell under small and large displacements, respectively.
GPa and mm = 0. The Young modulus and the cross-sectional area
The fourth one deals with a laminated shell reinforced with ran-
for fibers are given, respectively, by Ef = 210 GPa and Af = 2 cm2.
dom fibers under small and large displacements and the last one
The total load q = 50 kN/m is applied in 10 steps.
combines long and random short fibers for large displacement
Fig. 6 shows the final configuration of the beam after deformation
analysis.
and Fig. 7 compares the transverse displacements achieved with the
four considered formulations. The maximum deflections at the
6.1. Fiber reinforced cantilever beam under large displacements loaded end of the reinforced cantilever beam as well as the relative
differences between results are shown in Table 1. As one can see
In this example, the displacement achieved at the free end of a the solutions obtained with the present formulation and the homog-
uniformly loaded beam is compared with results obtained by three enized formulation are very close to the reference solution.
different formulations: (i) a Timoshenko–Reissner 3D frame for- As expected all solutions are more flexible than the 3D frame
mulation presented by [35–38]; (ii) a 2D solid formulation pre- solution, because it does not consider the adopted cross section
sented by [39]; and (iii) a laminated equivalent formulation, strain enhancement used in the shell formulation and naturally
named homogenized formulation, presented by [34]. present in the 2D solid solution.
The beam is reinforced by four stiffeners arranged as shown in
Figs. 4 and 5. When the beam is solved by the 2D solid formulation, 6.2. Fiber reinforced laminated cylindrical shell – small displacements
named reference formulation, the matrix is discretized by 60 trian-
gular cubic solid elements and 30 cubic fiber elements for each line In this example, the orthotropic laminated cylindrical pinched
of reinforcements, see Fig. 6(a). For the Timoshenko–Reissner 3D problem presented in [34,40–42] is solved considering four
 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814 809

0º/0º/0º

(a) (b)
90º/90º/90º

(a) (b)
0º/90º/0º

(a) (b)
90º/0º/90º

(a) (b)

Fig. 12. Final configuration of 1/8 of the cylinder in large displacements: (a) homogenized solution and (b) present solution.

different ply arrangements i.e., 0°/0°/0°, 90°/90°/90°, 0°/90°/0° and shell developing small displacements provided by [41]. For more
90°/0°/90°. In these ply arrangements, 0° corresponds to the longi- details see [34].
tudinal direction and 90° corresponds to the circumferential direc- The original geometry of the problem, extracted from [40], is
tion. It is important to mention that the reference solution used shown in Fig. 8. As given by the original reference, dimensions
here to validate the introduced formulation was compared by and other physical properties do not have units. The geometry val-
[34] with an analytical solution for laminate orthotropic cylindrical ues are: L = 600, R = 300 and h = 3.0. The orthotropic laminated
810 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

(a) 0º0º0º (b) 90º90º90º

(c) 0º90º0º (d) 90º0º90º

Fig. 13. Zoom of final configuration obtained with the present formulation.

Table 5
Geometrical data used to modeling fibers: (E)  106.

Lamina Fiber set Random fibers

Ef Af n3 Case 1 Case 2 Case 3
Nf Lf Nf Lf Nf Lf
1 1 426 0.457 +0.5 608 2.54 309 5 154 10
2 426 0.457 0.5 608 2.54 309 5 154 10
2 1 426 0.457 +0.5 608 2.54 309 5 154 10
2 426 0.457 0.5 608 2.54 309 5 154 10
3 1 426 0.457 +0.5 608 2.54 309 5 154 10
2 426 0.457 0.5 608 2.54 309 5 154 10

cylindrical pinched is divided into 3 layers of thickness hl = 1.0 To allow small displacement analysis a small load is adopted,
each. P = 2.5  107. Fig. 10 shows the final configuration of the cylinder
The physical properties adopted in the analysis are specified in when fibers are spread in two lines of each lamina for all consid-
Table 2. These data are used to obtain the reference solution, ered ply arrangements. Table 4 compares the achieved results for
named homogenized solution, and follow the same values used the non-dimensional deflection (EhW/nP) under the applied load
by the original reference [41]. Taking advantage of the convergence with the homogenized solution. In the non-dimensional deflection
analysis made by [34], the reference surface of the cylinder is dis- expression, W is the final displacement under the applied load and
cretized by 800 curved triangular cubic finite elements with 3721 n = 4 is the division of the load due to the adopted discretization,
nodes and 26047 degrees of freedom as shown in Fig. 8. that is, 1/8 of the cylinder.
For the analysis performed with the present formulation, two The case in which fibers are arranged in the midsurface of lam-
ways of spread fibers inside the layers are considered: (i) the fibers inas presents larger relative difference when compared to the
are positioned in the midsurface of each lamina as shown in homogenized solution. This occurs because the contribution of
Fig. 9(a) and (ii) the fibers are located in two positions away inertia bending moment of the central lamina is crucial to the total
±0.25 from the midsurface of each considered lamina, as shown inertia of the shell. When fibers are spread in the matrix in two
in Fig. 9(b). lines, the relative difference decreases, mainly for the ply arrange-
Table 3 shows the adopted data for the matrix and the num- ments 0°/90°/0° and 90°/0°/90°. Furthermore, as the number of
ber of fibers needed in each lamina to approximate the Young fibers increase the difference between the discrete and homoge-
modulus of the equivalent or homogenized problem when fibers nized approaches reduces. Therefore, using a discrete model, as
are located in the midsurface of the considered lamina. When the one proposed in this work, may be a good alternative to solve
the fibers are spread in two lines in each lamina the cross- the difficulty in defining homogenized parameters for laminated
sectional area and the Young modulus are given, respectively, shell analysis, as for example, the shear elastic modulus.
by Af = 0.457 and Ef = 426  106. The other properties are the For all analyzed cases the results achieved with the present for-
same as shown in Table 3. Each fiber is discretized by 40 cubic mulation are more flexible than the ones achieved using the homog-
finite elements and the reference surface discretization is the enized formulation. As the homogenized problem and the reinforced
same shown in Fig. 8. problem are different, equal results are not expected.
 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814 811

Table 6 Table 7
Physical properties used to model matrix. (E)  106 and (G)  106. Numerical results for the non-dimensional deflection of the cylinder reinforced with
random fibers under the applied load.
Lamina Matrix
Homogeneous solution Present solution
E m G
Lf= 2.54 Lf= 5 Lf= 10
E1 E2 E3 m12 m13 m23 G12 G13 G23
EhW/nP EhW/nP EhW/nP EhW/nP
1 1 1 1 0.25 0.25 0 0.4 0.4 0.5
6751.2 5206.44 5283.84 5761.8
2 1 1 1 0.25 0.25 0 0.4 0.4 0.5
Relative difference (%) 22.881 21.735 14.655
3 1 1 1 0.25 0.25 0 0.4 0.4 0.5

6.4. Laminated cylindrical shell reinforced by random short fibers

6.3. Fiber reinforced laminated cylindrical shell – large displacements
This qualitative example is used to demonstrate the potential of
This example compares the behavior of the reinforced shell ana- the proposed formulation. In this case, random short fibers are
lyzed in Example 6.2 for large displacements. The homogenized spread in the laminas of the pinched problem at the rate of 1%
case is compared to the discrete reinforced case considering two by volume of the cylinder, see Table 5. Three fiber lengths are con-
lines of reinforcement. In this case, the total load P = 0.25 is applied sidered, Lf = 2.54, Lf = 5.0 and Lf = 10. For all cases, the fiber Young
in 50 steps. The physical and geometrical properties adopted in the modulus and the cross-sectional area are given respectively by
performed analyses are the same as used for the linear case pre- Ef = 426  106 and Af = 0.457. A short fiber is discretize by only
sented in Table 3. one cubic fiber element. To keep the volume rate, the total number
Fig. 11 compares the equilibrium paths obtained by the present of fibers for each considered case is given, respectively, by 3568,
formulation with that one obtained with the homogenized formu- 1854 and 924. Fibers are spread in six equal sets, two by lamina,
lation for all considered ply arrangements. As one may observe, the away ±0.25 from the midline of the considered lamina. Table 5
arrangements 0o/0o/0o and 0o/90o/0o are more flexible than the ply shows the geometrical data used to model the fibers.
arrangements 90o/90o/90o and 90o/0o/90o. The displacements The analyses in small displacements were performed consider-
obtained with the present formulations are greater than those pre- ing P = 2.5  107. For large displacements analyses, the loads
sented by the homogenized formulation. These results show that P = 0.025 and P = 0.25 were adopted. In all cases the total load is
the mechanical coupling between fibers and the laminated shell applied in 50 steps. First the homogenous case is performed and
finite element proposed in this paper is working properly for large then fibers are considered. The data used to model the matrix are
displacement analyses. shown in Table 6.
Fig. 12 shows the final configuration for all considered arrange- Fig. 14 show the final configurations for 1/8 of the cylinder in
ments and Fig. 13 shows an enlarged picture of case (b), shown in small displacements for different considered fiber lengths and
Fig. 12. Table 7 shows the non-dimensional deflection (EhW/nP) under

(a) (b)

(c) (d)
Fig. 14. Final configuration of 1/8 of the cylinder reinforced by random fibers for small displacements: (a) homogeneous case; (b) Lf = 2.54; (c) Lf = 5.0 and (d) Lf = 10.
812 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

0.030 (a) 0.30 (b)

Homogeneous shell
Reinforced with random fibres: L = 2.54
0.025 0.25 Reinforced with random fibres: L = 5.0
Reinforced with random fibres: L = 10

0.020 0.20

Load
Load

0.015 0.15

0.010 Homogeneous shell 0.10

Reinforced with random fibres: L=2.54
Reinforced with random fibres: L=5.0
Reinforced with random fibres: L=10
0.005 0.05

0.000 0.00
0 2 4 6 8 10 12 14 16 0 20 40 60 80 100 120 140 160 180 200

Displacements Displacements

Fig. 15. Equilibrium path of the point under the applied load P for the laminated shell reinforced by random fibers: (a)P = 0.025; (b)P = 0.25.

(a) (b)

(c) (d)
Fig. 16. Final configuration of 1/8 of the cylinder reinforced by random fibers for large displacements and P = 0.25: (a) homogeneous case; (b) Lf = 2.54; (c) Lf = 5.0 and (d)
Lf = 10.

the applied load and the relative differences related to homoge- Observing Fig. 15(b) one may note the presence of some snap-
neous case. through points along the equilibrium paths and the large flexibility
For all cases the addition of fibers stiffened the laminated shell, of the laminated shell reinforced by random fibers for P = 0.25. This
however, it is noted that although the volume of fibers is the same behavior strengths the need to properly know the material proper-
in all considered cases, the achieved displacements for the longer ties and dispose of formulations that allow more realistic modeling
short fiber are larger than those achieved when shorter fibers are of the problem as the one introduced in this paper.
employed. Fig. 16 show the final configurations for 1/8 of the cylinder rein-
Fig. 15(a) and (b) compare the equilibrium path of the cylinder forced with random fibers for P = 0.25 and Fig. 17 shows selected
reinforced by random fibers with the homogeneous case for views of the final configuration of 1/8 of the cylinder reinforced
P = 0.025 and P = 0.25, respectively. with random fibers for Lf = 2.54 and P = 0.25.
 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814 813

Fig. 17. Some views of the final configuration of 1/8 of the cylinder reinforced by random fibers for Lf = 2.54 and P = 0.25.

Case 1 Case 2 Case 3

Fig. 18. Final configurations for all considered cases.

0.30

0.25

0.20
Load

0.15

0.10
Case 1
Case 2
0.05
Case 3

0.00
0 2 4 6 8 10 12 14 16 18
Displacements

Fig. 19. Equilibrium path of the point under the applied load P for the laminated shell reinforced with long and random short fibers.

6.5. Laminated cylindrical shell reinforced with long and random short
fibers

In this example the ply arrangement 90°/90°/90° of Example 6.2

is reinforced by random short fibers in three different ways. In case
1 short fibers are spread in three lamina; in case 2 short fibers are
spread only in the central lamina and in case 3 long fibers are
excluded from the middle layer of case 2. The considered rates of
short fibers are 1% of the cylinder volume for case 1 and 0.33%
for the other cases.
The analysis is performed considering P = 0.025. The total load
is applied in 50 steps. The physical and geometrical properties
are the same of the previous examples for both matrix, long and
short fibers. Fig. 20. Zoom of final configuration for case 1.
814 M.S.M. Sampaio et al. / Composite Structures 119 (2015) 799–814

Fig. 18 shows the final configurations of the 1/8 of the cylinder [12] Liu Y, Soh C-K. Shear correction for Mindlin type plate and shell elements. Int J
Numer Methods Eng 2006;69:2789–806.
for all considered cases and Fig. 19 compares the equilibrium path
[13] Zhang Y, Yang C. Recent developments in finite element analysis for laminated
of the cylinder reinforced with long and random fibers for the three composite plates. Compos Struct 2009;88:147–57.
considered cases. As one can observe the proposed formulation [14] Ramesh SS, Wang CM, Reddy JN, Ang KK. A higher-order plate element for
presents coherent results, i.e., the less dense reinforcement case accurate prediction of interlaminar stresses in laminated composite plates.
Compos Struct 2009;91:337–57.
is more flexible than the other cases. Finally, Fig. 20 shows an [15] Yang H, Saigal S, Liaw D. Advances of thin shell finite elements and some
enlarged Picture of case 1 shown in Fig. 18. applications – version I. Comput Struct 1990;35:481–504.
[16] Yang HTY, Saigal S, Masud A, Kapania RK. A survey of recent shell finite
elements. Int J Numer Methods Eng 2000;47:101–27.
7. Conclusion [17] Qatu M, Asadi E, Wang W. Review of recent literature on static analyses of
composite shells: 2000–2010. Open J Compos Mater 2012;2:61–86.
[18] Qatu M, Sullivan R, Wang W. Recent research advances on the dynamic
The mathematical and numerical modeling of fiber reinforced
analysis of composite shells: 2000–2009. Compos Struct 2010;93:14–31.
laminated plates and shells is a very challenging and complex task [19] Li D, Qing G, Liu Y. A layerwise/solid-element method for the composite
because it involves a great number of variables, different fiber stiffened laminated cylindrical shell structures. Compos Struct
2013;98:215–27.
arrangements in the laminas and the need to correctly know the
[20] Zudans Z. Analyses of asymmetric stiffened shell type structures by the finite
material parameters. In this study we successfully propose and element method: 1. Flat rectangular elements. Nucl Eng Des 1968;8:367–79.
implement a new FE formulation for the inclusion of long and ran- [21] Zudans Z. Analysis of asymmetric stiffened shell type structures by the finite
dom fibers in any lamina of laminated plates and shells that element method: III. Constant transverse shear model. Nucl Eng Des
1970;11:177–94.
develop small or large deformations. The proposed formulation [22] Venkatesh A, Rao K. Analysis of laminated shells with laminated stiffeners
naturally results in the reinforced continuum behavior. This for- using rectangular shell finite elements. Comput Methods Appl Mech Eng
mulation does not increase the number of degrees of freedom of 1983;38:255–72.
[23] Venkatesh A, Rao K. Analysis of laminated shells of revolution with laminated
the original homogeneous discretization and does not require the stiffeners using a doubly curved quadrilateral finite element. Comput Struct
necessity of matching nodes in the discretization of fibers and 1985;20:669–82.
matrix. Numerical examples show the good behavior and the [24] Bhimaraddi A, Carr A, Moss P. Finite element analysis of laminated shells of
revolution with laminated stiffeners. Comput Struct 1989;33:295–305.
potential of the proposed formulation. Furthermore, this formula- [25] Liao C-L, Reddy J. Analysis of anisotropic, stiffened composite laminates using
tion may be an alternative strategy to homogenization techniques a continuum-bases shell element. Comput Struct 1990;34:805–15.
in view of the great difficulty in determining representative elastic [26] Li LY, Applegarth I, Bull JW, Bettes P, Bond TJ, Thompson PA. An auto-adaptive
finite element analysis software for stiffened plate and shell structures. Adv
properties for orthotropic or anisotropic medium, mainly when
Eng Softw 1997;28:285–91.
fibers are randomly spread in the analyzed body. In future works [27] Li LY, Bettes P, Bull JW, Bond TJ. Adaptive finite element analysis of stiffened
fiber and matrix physical nonlinearities should be investigated, shells. Adv Eng Softw 1997;28:501–7.
[28] Samanta A, Mukhopadhyay M. Finite element large deflection static analysis of
as well, the fiber–matrix debonding.
shallow and deep stiffened shells. Finite Elem Anal Des 1999;33:187–208.
[29] Barut A, Madenci E, Tessler A, Starnes Jr JH. A new stiffened shell element for
Acknowledgments geometrically nonlinear analysis of composite laminates. Comput Struct
2000;77:11–40.
[30] Prusty B. Linear static analysis of composite hat-stiffened laminated shells
The authors thank the Coordination for the Improvement of using finite elements. Finite Elem Anal Des 2003;39:1125–38.
Higher Education Personnel (CAPES); the São Paulo Research Foun- [31] Prusty B, Satsangi S. Analysis of stiffened shell for ships and ocean structures
dation (FAPESP) and the Amazon Research Foundation (FAPEAM) by finite element method. Ocean Eng 2001;28:621–38.
[32] Wagner W, Balzani C. Simulation of delamination in stringer stiffened fiber-
for the financial support of this research. reinforced composite shells. Comput Struct 2008;86:930–9.
[33] Coda HB, Paccola RR. A positional FEM formulation for geometrical non-linear
References analysis of shells. Latin Am J Solids Struct 2008;5:205–23.
[34] Coda HB, Paccola RR, Sampaio MSM. Positional description applied to the
solution of geometrically non-linear plates and shells. Finite Elem Anal Des
[1] Lo K, Christensen R, Wu E. A high-order theory of plate deformation. Part 2:
2013;67:66–75.
laminated plates. J Appl Mech 1997;44:669–76.
[35] Coda HB, Greco M. A simple FEM formulation for large deflection 2D frame
[2] Levison M. An accurate, simple theory of the statics and dynamics of elastic
analysis based on position description. Comput Methods Appl Mech Eng
plates. Mech Res Commun 1980;7:343–50.
2004;193:3541–57.
[3] Reddy J. A simple higher-order theory for laminated composite plates. J Appl
[36] Coda HB. A solid-like FEM for geometrically non-linear 3D frames. Comput
Mech 1983;51:745–52.
Methods Appl Mech Eng 2009;198(47-48):3712–22.
[4] Bert C. A critical evaluation of new plate theories applied to laminated
[37] Coda HB, Paccola RR. Improved finite element for 3D laminate frame analysis
composites. Compos Struct 1984;2:329–47.
including warping for any cross-section. Appl Math Model
[5] Reddy J. A refined nonlinear theory of plates with transverse shear
2010;34(4):1107–37.
deformation. Int J Solids Struct 1984;20:881–96.
[38] Coda HB, Paccola RR. A FEM procedure based on positions and unsconstrained
[6] Phan N, Reddy J. Analysis of laminated composite plates using a higher-order
vectors applied to non-linear dynamic of 3D frames. Finite Elem Anal Des
shear deformation theory. Int J Numer Methods Eng 1985;21:2201–19.
2011;47(4):319–33.
[7] Reddy J, Liu C. A higher-order shear deformation theory of laminated elastic
[39] Sampaio MSM, Paccola RR, Coda HB. Fully adherent fiber–matrix FEM
shells. Int J Eng Sci 1985;23:319–30.
formulation for geometrically nonlinear 2D solid analysis. Finite Elem Anal
[8] Kant T, Kommineni J. C0 finite element geometrically non-linear analysis of
Des 2013;66:12–25.
fiber reinforced composite and sandwich laminates based on a higher-order
[40] Sansour C, Kollmann F. Families of 4-nodes and 9-nodes finite elements for a
theory. Comput Struct 1992;45:511–20.
finite deformation shell theory. An assessment of hybrid stress, hybrid strain
[9] Liangxin S, Zhiyu S. The analysis of laminated composite plates based on the
and enhanced strain elements. Comput Mech 2000;24:435–47.
simple higher-order theory. Comput Struct 1992;43:831–7.
[41] Jones I. Pinched cylinder benchmarks for shear deformable laminated shells.
[10] Kant T, Kommineni J. Geometrically non-linear transient analysis of laminated
Plast Rubber Compos 2002;31(6):249–61.
composite and sandwich shells with a refined theory and C0 finite elements.
[42] Kulikov GM, Plotnikova SV. Non-linear geometrically exact assumed stress-
Comput Struct 1994;52:1243–59.
strain four-node solid-shell element with high coarse-mesh accuracy. Finite
[11] Arciniega R, Reddy J. Consistent third-order shell theory with application to
Elem Anal Des 2007;43:425–43.
composite circular cylinders. Am Inst Aeronaut Astronaut (AIAA) J
2005;43(9):2024–38.