2018 edition, Scopus 2019

Frontiers journals are at the top of citation and impact metrics

Original Research ARTICLE

Front. Built Environ., 23 October 2018 |

A Finite Element Analysis of the Stability of Composite Beams With Arbitrary Curvature

  • 1Department of Civil Engineering, University of Salerno, Fisciano, Italy
  • 2Department of Structures for Engineering and Architecture, University of Naples Federico II, Naples, Italy

A finite element approximation of a theory recently proposed for the geometrically nonlinear analysis of laminated curved beams is developed. The application of the given finite element model to the computation of stability points and post-buckling behavior of beams with arbitrary curvature is also carried out, on taking into account the influences of shear deformation and warping effects on the in-plane and out-plane responses of the beam. The stability analysis is performed through a path-following procedure and a bordering algorithm. Several numerical results are given and comparisons with classical beam theories and other theories available in the relevant literature are established. The given results highlight that the proposed finite element model is well suited to study the stability of structures that incorporate laminated composite beams, such as, e.g., light-weight roof structures and arch bridges.


Laminated composite (fiber-reinforced) structures are increasingly used in a wide range of engineering applications (naval, aeronautical, automation, mechanical, civil, medical engineering, etc.), due to the outstanding proprieties that such structures may exhibit when a proper design of the material and the lamination scheme are employed: light weight, high stiffness-to-weight and tensile strength-to-weight ratios, high damping, excellent corrosion, thermal and high impact resistance (Fraternali et al., 2011, 2012; Bencardino et al., 2012). Nowadays, laminated composite structures play a crucial role in the production of various innovative structures or products, which include: light-weight roof structures, arch bridges, impact energy mitigation and vibration isolation devices, just to name a few examples. In order to capture the puzzling mechanical response of such structures, various theoretical and numerical approaches have been proposed in the literature, including zig-zag displacement-based theories and stress-based methods, with special attention on the modeling of anisotropy, warping, fracture and damage (Feo and Fraternali, 2000; Roberts and Al-Ubaidi, 2001; Fraternali et al., 2002, 2010, 2011, 2012; Fraternali, 2007; Schmidt et al., 2009; Feo and Mancusi, 2010; Bencardino et al., 2012; Markkula et al., 2013; Viera et al., 2013; Özütok and Madenci, 2017). Shear deformation and warping effects may be rather important in composite beams. At this regard (Özütok and Madenci, 2017) by Özütok and Madenci analyses the effects of a non-linear distribution of the shear stress through the beam thickness within a higher-order shear deformation theory; Roberts and Al-Ubaidi develop in Roberts and Al-Ubaidi (2001) an approximate theory for assessing the influence of shear deformation on restrained torsional warping of pultruded bars; while Viera et al. develop a thin-walled beam model able to simulate warping and higher order effects in Viera et al. (2013). For what concerns stability phenomena, which are of peculiar interest in the case of composite beams due to the characteristic slenderness of such structures, it is worth mentioning the non-linear elastic approaches proposed in Ascione et al. (2011, 2013); Fraternali et al. (2013); Mascolo and Pasquino (2016); Özütok and Madenci (2017), and Mascolo et al. (in press). The recent study illustrated in Fraternali et al. (2013) presents a geometrically nonlinear theory of laminated curved beams, which assumes that cross-section rotations and shear strains are moderately large, while axial strains are infinitesimal. Due to its minor complexity with respect to the finite elasticity theory, the model presented in Fraternali et al. (2013) is particularly convenient for computing the first stability point of a composite laminated beam and studying its behavior near such a point.

In the present paper, we develop a comprehensive finite-element approximation of the mechanical model given in Fraternali et al. (2013), which is founded upon the use of Lagrangian isoparametric elements (section Finite Element Model). The adopted model proves to be a robust and versatile tool that allows to model the geometrically non-linear response and the buckling behavior of laminated composite beams with arbitrary curvature. It takes into account both shear deformations and warping effects, which are essential to accurately predict in-plane and out-plane buckling loads. The proposed stability analysis employs the path-following procedure proposed by Bathoz and Dhatt (1979) and the algorithm for computing stability points proposed by Simo and Wriggers (1990) (section Finite Element Analysis of the Stability of Composite Curved Beams). The accuracy of the proposed finite-element model is assessed by presenting different numerical results relative to the stability of isotropic and composite beams and establishing comparisons with the corresponding results of classical beam theories and other theories available in the relevant literature (section Numerical Results). We end with concluding remarks and directions for future work in section Concluding Remarks.

Finite Element Model

Let us denote by Ch a finite-element discretization of axis curve of a laminated beam, and let us assume that the elements C1, …, Cn belong to the Lagrange family (Reddy, 1992) (see Figure 1)

Ch=e=1neCe.    (1)

Figure 1. Finite element discretization of the axis of a laminate curved beam (nl: numbers of layers).

On adopting an isoparametric finite-element approximation (Reddy, 1992) we use the same shape functions to approximate both the geometry and the (generalized) displacement field over the generic element Ce

Z2eh=I=1nNIZ2I, Z3eh=I=1nNIZ3I, u^eh=I=1nNIu^I.    (2)

In Equation (2) NI is the shape function corresponding to the node I and consists of a complete polynomial of order n − 1; Z2I and Z3I are the coordinates of I with respect to the global frame {0, Z1, Z2, Z3} (Figure 1); u^I is the generalized displacement vector relative to the same node.

In particular, for a four-node Lagrangian element, we represent in Figure 2 the transformation (2)1,2 which maps the master element onto a curved (cubic) element.


Figure 2. Standard (A) and distorted (B) four-node Lagrangian element.

The Jacobian of the transformation from the local coordinate ξ (Figure 2) to the global coordinate X3 (Figure 1) is given by Nomizu and Kobayashi (1963)

J=dX3dξ=(dZ2ehdξ)2+(dZ3ehdξ)2   =(I=1nNI,ξZ2I)2+(I=1nNI,ξZ3I)2    (3)

NI, ξ being the derivative of NI with respect to ξ. Denoting by (·)′ the derivative with respect to X3, we have:

NI=dNIdX3=J1N1, ξ.    (4)

By making use of Equation (2)1,2, we obtain the following approximation of the curvature radius R (Nomizu and Kobayashi, 1963)

Reh=[(Z2eh),ξ2+(Z3eh),ξ2]3/2[(Z2eh),ξ(Z3eh),ξξ(Z2eh),ξξ(Z3eh),ξ]=[(I=1nNI,ξZ2I)2+(I=1nNI,ξZ3I)2]3/2[(I=1nNI,ξZ2I)(J=1nNJ,ξξZ3J)(I=1nNI,ξξZ2I)(J=1nNJ,ξZ3J)].    (5)

Coming back to Equation (2)3, we now observe that it can be written in the following compact form

u^eh= NUe    (6)

where Ue is the M-dimensional vector collecting nodal, generalized, displacements of

Ce(Ue=[u^1T, u^2T, ,u^nT]T),    (7)

while N is the following matrix

N[m×M]=[N1, N2, , Nn],    (8)

whose blocks are diagonal submatrices

NI[m×m]=diag(NI NI  NI).    (9)

Equation (6) leads us to obtain the following approximation of the generalized strains

E^eh(Ue)=E^e(1)h(Ue)+12E^e(2)h(Ue, Ue)    (10)


E^e(1)h(Ue)=BoUeE^e(2)h(Ue, δUe)=BL(Ue)δUe    (11)
B0 and BL(Ue) being the following σ×M matricesB0=[B01, B02, , B0n ],BL(Ue)=A(Ue)G.    (12)

The σ × m submatrices B0I in Equation (12)1 are given by

BoI=[BoIEνBoIEφ00BoI Θϕ000000BoIww00BoIww]   ,    (13)


BoIEν=BoI Θϕ=[N1000N1N1/R0N1/RN1],BoIEφ=[0-N10N100000],BoIww[mw×mw]=diag(NI NI NI),BoIww[mw×mw]=diag(N1 N1 N1).    (14)

The matrices A and G, which appear in Equation (12)2, are instead given by

A(Ue)[σ×9]=[A(Ue)0 ],        G[9×M]=[G1, G2, , Gn],      (15)


GI[9×m]=[G¯I, 0],G¯I[9×6]=[N1000000N1000000N1000000N1000000N1000000N1000N1000000N1000000N1],    (16)
A(Ue)[9×9]=[0ϕ3ϕ3Rϕ320(v2+v3R+ϕ12)000ϕ3000ϕ32(v1+ϕ22)000v1(v2+v3R2)(v2R+v3R2)000000(ϕ3+ϕ2R)000(v1R+ϕ2R+ϕ32)(ϕ22+ϕ3R)0ϕ32(v1+ϕ22)0(ϕ3+ϕ2R)(ϕ3R+ϕ2R2)(ϕ32+ϕ22R)(v22+v3R2+ϕ12R)ϕ12ϕ320(v2v3Rϕ12)000(ϕ22+ϕ32R)ϕ12ϕ12Rϕ22ϕ1200000(ϕ3R+ϕ2R2)(ϕ2R+ϕ3R2)0(ϕ2+ϕ3R)(ϕ3+ϕ2R)0000(ϕ3R+ϕ2R2)0ϕ10(ϕ3+ϕ2R)00000ϕ1R(ϕ2+ϕ3R)ϕ10]     (17)

Finite Element Analysis of the Stability of Composite Curved Beams

Path-Following Procedure

Let us consider the variational formulation of the equilibrium equations. The use of Equations (6, 10, 11) allows us to set such equations into the following discrete form

δΠh=δUTR(U,λ)             =e=1ne{δUeT11[B0T+BLT(Ue)]D^[B0+12BL(Ue)]UeJdξ}             λe=1ne{δUeT11NT[q^e(1)+q^e(2)(Ue)]Jdξ}              λI=1nn{δu^IT[Q^I(1)+Q^I(2)(u^I)]}=0    (18)

where Πh is the discretized functional of the total potential energy Π; δΠh is the first variation of Πh with increment δU; R(U, λ) (residuai vector) is the Gateaux derivative of Πh with respect to U(R=DUΠh).

On accounting for a possibility of non-linear elastic response of the material, we assume that the elasticity matrix D^, whose elements are the resultants and the resultant moments of the local elastic moduli (Fraternali et al., 2013), depends on the deformation of the beam (D^ = D^(U)).

Equation (18) can also be written in the following compact form

δUTR(U,λ)=δUT{K(U)Uλ[Q(1)+Q(2)(U)]}=0,    (19)

where K(U) is the N × N global (secant) stiffness matrix, which derives from the assembly of the element stiffness matrices Ke (e = 1, 2, …, ne). Such matrices are defined by the equations

Ke(Ue)(Ue)=11[B0T+BLT(Ue)]S^(Ue)Jdξ,    (20)

where S^ is the generalized stress vector

S^(Ue)=D^[B0+12BL(Ue)]Ue,    (21)

while Q(1) and Q(2)(U) are the global force vectors, which derive from the assembly of the element vectors

Qe(1)=11NTq^e(1)Jdξ,          Qe(2)(Ue)=11NTq^e(2)(Ue)Jdξ,    (22)

and the nodal force vectors Q^I(1) and Q^I(2) (u^I)  (I=1, 2, , nn).

Due to the arbitrariness of δU, Equation (19) is equivalent to the following nonlinear system of N equations

R(U,λ)=K(U)U λ[Q(1)+ Q(2)(U)]=0,    (23)

which can be solved by employing one of the algorithms known in literature as path-following methods (for an overview of such methods (see e.g., Riks, 1972). The basic idea of path-following methods is to append a constraint Equation f(U, λ) = 0 to (23). Here, following Bathoz and Dhatt (1979), we adopt a displacement control and assume

f(U)=epTUu=0,    (24)

where ep is the vector of ℜN which has only the pth component different from zero and equal to 1, while μ is a prescribed value of the pth component of U. Therefore, we are led to solve the extended system

R( U,λ)={R( U,λ) epTU μ}=0,    (25)

The linearization of Equation (25) by the Newton-Raphson method gives the system of incremental equilibrium equations

R(U+Δ U,λ+Δλ)= R(U,λ)+[DURDλR epT0]{ΔUΔλ}=0,    (26)

The matrix DUR, which is usually referred to as tangent stiffness matrix and denoted by KT, is given in the Appendix. Concerning the vector DλR, from Equation (23) we deduce

DλR=[ Q(1)+ Q(2)(U)].    (27)

The non-symmetric system (26), that we rewrite in the form

[KT(Q(1)+ Q(2)) epT0]{Δ UΔλ}={ R( U,λ) epTU μ},    (28)

can be solved by a procedure known in literature as bordering algorithm (see e.g., Keller, 1977). Such algorithm and the overall procedure for the solution of the extended system (25) are described in Table 1.


Table 1. Bordering algorithm.

In particular, if the first predictor U satisfies the constraint (24), Equation (31) reduces to

Δλ=epT Δ URepT Δ UQ    (34)

We point out that, since we proceed by displacement control, we apply the above iterative procedure in incremental steps. Within the generic step, say the ith one, we increment by δ the displacement component Up which exhibited the largest variation in the previous step. We hence set in the extended system (25)

μ=epT Ui1+δ    (35)

and begin the new iteration loop by assuming the predictor U~= Ui-1+δep, λ~=λi-1, which satisfies Equation (24).

Computation of Stability Points

In the current section we get a finite-element approximation of the problem of computing stability points based on the Trefftz criterion (Trefftz, 1930).

Within the previous settings, we obtain the following discrete equation

DU2eh(U,λ)U1δU=δUTKT(U,λ)U1=0    (36)

which must be satisfied by every variation δU.

A point U, λ such that Equation (36) holds for some U1 is usually called stability point, while U1 is called buckling mode (or eigenvector) associated with U, λ.

Due to the arbitrariness of δU, Equation (36) is equivalent to the system of N Equations

KT( U,λ)U1=0    (37)

In the following we will denote U1 by V. In order to exclude the trivial case V=0, it is necessary to append a constraint equation l(V)=0 to system (37). Possible choices of such an equation are

V1 = 0    (38)
epTVV0 = 0    (39)

V0 being a fixed (non-zero) value of the pth component of V. In this work we make use of Equation (39) and, after having reduced KT to an upper triangular matrix (by Gauss elimination technique), we identify the index p with the equation number where the lowest diagonal term of the reduced stiffness matrix appears. In this way we prevent the pth component of V becoming exceedingly large. Concerning V0, we set

V0=epTV0 V0    (40)

V0 being the initial approximation to V.

Stability points can be classified in limit (or turning) and bifurcation points (Budiansky, 1974). Following Spence and Jepson (Spence and Jepson, 1985) we can distinguish between the two cases by using the following criteria

Bifurcation point:    VT(Q(1)+ Q(2)(U))=0    (41)
Limit point:               VT(Q(1)+ Q(2)(U))0    (42)

An efficient procedure for the computation of stability points has been proposed by Simo and Wriggers (1990). It consists of solving the extended system

R(U,V,λ,μ)={R(U,λ)KT(U,λ)VepTVV0epTUμ}=0    (43)

which derives from the addition of Equations (37, 39) to the system of non-linear equilibrium Equations (25).

Since the tangent stiffness matrix becomes progressively ill-conditioned as the solution approaches the stability point, where KT is singular, from a numerical point of view it is convenient to transform the extended system (43) into the following equivalent form Simo and Wriggers (1990)

Rη(U,V,λ,μ)={R(U,λ)+η(epTUμ)epKT(U,λ)V+η(epTVV0)epepTVV0epTUμ}=0,    (44)

η being an arbitrary positive number. The linearization of Equation (44) by the Newton-Raphson method leads us to obtain the set of incremental equations

[KTη0QηepDU(KTV)KTηDλ(KTV)00TepT00epT0T01]{ΔUΔVΔλΔμ}={R(U,λ)+η(epTUμ)epKTη(U,λ)VηV0epepTVV0epTUμ}    (45)

where Q=Q(1) + Q(2)(U), while

KTη(U,λ)=KT(U,λ)+ηepepT (p not summed)    (46)

is a rank-one updated stiffness matrix. The solution of the non-symmetric system (45) can be obtained by a bordering algorithm similar to that described in the previous section Path-Following Procedure. Table 2 shows this algorithm and the global procedure for the solution of the extended system (44).


Table 2. Bordering algorithm.

We furnish the expressions of the vectors hj (j = 1, …, 4), which appear in Equations (50), in the Appendix.

Overall Algorithm for Stability Analysis

We compute the pre-buckling and post-buckling equilibrium paths of a laminated beam by combining the procedures described in Sections Path-Following Procedure and Computation of Stability Points.

More precisely, denoted the initial stiffness matrix by K0 (see the Appendix) and an arbitrary numeric value by λ0, we consider the couple λ0, U0=λ0K0-1Q(1) as the initial predictor of the first equilibrium state.

We hence correct this predictor as described in section Path-Following Procedure and keep computing equilibrium states along the primary path. We check for the sign of the tangent stiffness matrix determinant in correspondence of each state, which is a simple operation since the solution of the extended system (25) requires the factorization (i.e., the triangular decomposition) of KT.

If the sign of det KT changes between two successive states, say i and i + 1, a stability point has passed. We hence switch from the path-following procedure to the procedure for computing stability points.

In particular we assume U~=Ui+1, V~=V0=K0-1ep, λ~=λi+1, μ~=epTUi+1 (for the meaning of the index p see the beginning of section Computation of Stability Points).

Once the stability point Uc, λc has been computed, we check if it is a limit or a bifurcation point. In the case of a limit point we come back to the path-following procedure to complete the primary path. In the case of a bifurcation point, we switch to the secondary (or bifurcated) path by adding to Uc a vector proportional to the eigenvector V

U=Uc+ζ V V    (55)

ζ being a scaling factor to be determined in such a way that it results R(Uc) tol.

We then follow the secondary path using the path-following procedure and arrest the calculations when the cross-section rotations or the shear strains are more than moderately large or the axial strains are more than infinitesimal (18).

Numerical Results

We present in this section several numerical results relative to the evaluation of the stability points and to the post-buckling behavior of straight and curved beams.

In all the examples we supposed the beam cross-section to be rectangular with lengths H1 and H2 along the directions X1 and X2, respectively. We denoted the cross-section area by A = H1H2, the moments of inertia by I1 and I2, the polar moment of inertia by IG and the De Saint Venant torsional rigidity by Jt

I1=H1H2312,I2=H13H212,IG=I1+I2,Jt=H13H23.    (56)

Concerning the expression of the warping function w, we examined the following three cases:

No Warping (NW): w = 0

Warping Function (W1): w = w11X1X2

Warping Function (W3): w=w20X12+w11X1X2+w02X22+


We used cubic Lagrangian finite elements and a four-point Gauss quadrature formula to compute the tangent stiffness matrix and its derivatives. By this choice we avoided the numerical inconvenient known in literature as shear and membrane (or inplane) locking (see e.g., Prathap and Bhashyam, 1982; Reddy and Averill, 1990).

We always assumed that the external loads retain their directions during the deformation of the beam (dead loading).

Bifurcation Points of Isotropic, Straight and Curved Beams

In order to assess the accuracy of our numerical model, we firstly present some results concerned with bifurcation points of isotropic straight and curved beams. They can be compared with those available in the relevant literature and corresponding to classical beam theories (see e.g., Timoshenko and Gere, 1961; Brush and Almroth, 1975). We assumed a ratio E/G = 0.385 between Young's and shear moduli.

The first example deals with a circular ring submitted to a radial dead load of intensity q, which is uniformly distributed along the centerline. We discretized one half of the ring by 20 finite elements imposing the following boundary conditions (no warping was considered)

v1=v2=v3=ϕ1=ϕ2=ϕ3=0 for X3=0,    v1=v3=ϕ1=ϕ2=ϕ3=0, for  X3=πR,    (57)

where R is the initial radius of the centerline. In particular, the ratios H1/H2 = 2, R/H2 = 20 were considered.

According to Donnel' s theory (see e.g., Brush and Almroth, 1975), the first bifurcation point occurs at a load level qbif= 4 E I1/R3 and the buckling mode corresponds to an ovalization of the ring.

It has to be remarked that the first bifurcation point occurs at a sensibly different load level qbif= 3 E I1/R3 if the external load remains orthogonal to the axis during the deformation of the ring (see e.g., Timoshenko and Gere, 1961; Brush and Almroth, 1975).

Table 3 shows the numerical convergence of the solutions of the extended system (44) for increasing values of the parameter η.


Table 3. Convergence behavior or the first bifurcation point or a circular ring submitted to a radial dead load q uniformly distributed along the centerline (H1/H2 = 2, R/H2 = 20, 20 element mesh).

Denoting by Dp the least diagonal term of the factorized tangent stiffness matrix and by Dop the corresponding term in the factorized initial stiffness matrix K0, in our numerical experiments we set η = 0, η = (DopDp) × 10 and η = (DopDp) × 1000.

As already observed by Simo and Wriggers (1990), for η = 0 the results exhibit oscillations near the bifurcation point, while for η > 0 they converge in a stable way to the value 3.9888 EI1/R3.

Table 3 also shows that in the latter case the solution is rather insensitive to the value of η.

The small difference existing between our and Donnel's bifurcation load can be justified observing that Donnel's theory does not account for shear deformation and for the quadratic terms in X1 and X2 of the axial strain E33, which are instead present in our model. In all the subsequent numerical applications we η = (DopDp) × 10.

The second example we considered deals with the classical case of a simply supported straight beam loaded by a compressive force at one end (X3 = L). Table 4 compares the bifurcation loads computed by using the present theory with those corresponding to Euler's theory (PEul=π2 EI1/L2) and Timoshenko's theory, for several the values of the ratio L/H2, assuming H1/H2 = 2. Upon putting ρ = χ PEul/GA (χ = 1.2 shear correction factor), we considered both the exact values deriving from Timoshenko's Theory

PbifPEul=(4+3ρ)+(4+3ρ)2+16ρ2ρ    (58)

and the approximate ones (see e.g., Timoshenko and Gere, 1961)

PbifPEul=11+ρ.    (59)

Concerning the comparisons between our and classical theories, we stressed the influences of warping and quadratic terms in X1, X2 of the axial strain E33 (Γ terms). It has to be remarked that our theory turns into a first-order shear deformation theory with no shear correction factor (χ = 1) in absence of warping and Γ terms.


Table 4. First bifurcation load of a simply supported, axially loaded straight beam for several ratios L/H2 (H1/H2 = 2, 20 element mesh).

Table 4 shows that our results corresponding to a cubic warping function and absence of Γ terms closely approximate those by Timoshenko. Furthermore, the influence of Γ terms is found to be appreciable (up to 2.7%) for thick beams. We took into account Γ deformation terms in all the successive numerical applications.

We complete this first group of numerical results by showing some further examples concerned with lateral buckling of straight bars and semicircular arches.

We firstly considered a narrow simply supported beam transversally loaded at the middle point and a narrow cantilever transversally loaded at the free end (H1/H2 = 0.1, L/H2 = 10). The kinematical boundary conditions of the first case are

v1=v2=v3=ϕ3=0 for X3=0,          v1=v2=ϕ3=0 for X3=L,    (60)

Three positions of the load were considered: load at the centroid, load at the extrados (X1 = 0, X2 = − H2/2) and load at the intrados (X1 = 0, X2 = H2/2). The last two cases obviously give rise to a deformation-dependent loading (Q(2)(U) ≠ 0).

Table 5 shows a comparison between the first bifurcation points computed within the present theory and those corresponding to the classical Prandtl's theory (See e.g., Timoshenko and Gere, 1961).


Table 5. First lateral bifurcation point of a simply supported beam loaded by a transverse force Q at the middle point and of a cantilever loaded by a transverse force Q at the free end (H1/H2 = 0.1, L/H2 = 10, 20 element mesh).

In the context of the present theory, we adopted a bilinear warping function (W1) and assumed both a linear and a (geometrically) nonlinear pre-buckling behavior (the first was treated by discarding the part KL of the tangent stiffness matrix, see the Appendix).

We also generalized Prandtl's theory in order to include a nonlinear pre-buckling behavior. This was obtained by using our model, neglecting shear deformation, warping and Γ terms and making the torsional stiffness equal to GJt.

It has to be remarked that the present theory assumes a torsional rigidity equal to GIG (as in the De Saint Venant theory of torsion without warping) and takes into account warping effects by introducing in the displacement field a polynomial warping function.

It is evident that our numerical results agree better with those by Prandtl when the warping is unconstrained, as in the case of the simply supported beam. On the contrary, in the second example (cantilever), our and Prandtl's results somewhat differ due to the presence of the restraint which doesn't allow the warping of the built-in end.

Further on we remark that the assumption of a linear pre-buckling behavior, as in the original Prandtl's theory, is less accurate for the cantilever than for the simply supported beam. Indeed, in the first case the in-plane displacements are sensibly higher than in the second case and produce moderate rotations.

Finally, Table 6 shows the bifurcation loads of hinged and clamped narrow semicircular arches (v1 = v2 = v3 = ϕ3 = 0 for X3 = 0, πR in the first case). The load consists of a uniform radial dead load along the centerline.


Table 6. First lateral bifurcation load of a hinged and a clamped semicircular arch loaded by a radial dead load q uniformly distributed along the centerline (H1/H2 = 0.1, R/H2 = 10, 20 element mesh).

The results of the present theory are compared with those given in Timoshenko and Gere (1961), which are relative to the classical beam theory. In particular, the first ones correspond to the choice of a bilinear warping function (W1).

Concluding Remarks

This work has developed a finite element model of the moderate rotation theory (MRT) of laminated composite beams proposed in Fraternali et al. (2013), and its application to the computation of nonlinear equilibrium paths and stability points of a variety of numerical examples. The proposed model describes laminated composite beams with arbitrary curvature of the beam axis, and takes into account shear deformation, warping effects, in-plane and out-of-plane instability. For the straight and curved beams examples analyzed in the present study, we conclude the following:

(i) the moderate rotation theory (MRT) and the classical beams theories correlate very well in almost all isotropic cases;

(ii) the influence of warping effects on the bifurcation load is generally pretty high for beams made up of composite materials.

Future work will apply the model presented in this work to a wide collection of technically relevant case-studies, with special emphasis on the study of the stability of light-weight roof structures and arch bridges, which make use of laminated composite beams (Fraternali et al., 2013).

Author Contributions

IM led the numerical part of the study. MM supervised the numerical part of the study. AA led the mechanical modeling part of the study. FF supervised all tasks.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


IM, AA, and FF gratefully acknowledge financial support from the Italian Ministry of Education, University and Research (MIUR) under the Departments of Excellence grant L.232/2016.

Supplementary Material

The Supplementary Material for this article can be found online at:


Ascione, L., Berardi, V. P., Giordano, A., and Spadea, S. (2013). Local buckling behavior of FRP thin-walled beams: a mechanical model. Comp Struct. 98, 111–120. doi: 10.1016/j.compstruct.2012.10.049

CrossRef Full Text | Google Scholar

Ascione, L., Giordano, A., and Spadea, S. (2011). Lateral buckling of pultruded FRP beams. Compos B Eng. 42, 819–824. doi: 10.1016/j.compositesb.2011.01.015

CrossRef Full Text | Google Scholar

Bathoz, I. L., and Dhatt, G. (1979). Incremental displacement algorithms for non-linear problems. Int. J. Num. Meth. Engng. 14, 1262–1267. doi: 10.1002/nme.1620140811

CrossRef Full Text | Google Scholar

Bencardino, F., Rizzuti, L., Spadea, G., and Swamy, R. N. (2012). Implications of test methodology on post-cracking and fracture behavior of steel fibre reinforced concrete. Compos. B Eng. 46, 31–38. doi: 10.1016/j.compositesb.2012.10.016

CrossRef Full Text

Brush, D. O., and Almroth, B. O. (1975). Buckling of Bars, Plates and Shells. New York, NY: McGraw-Hill.

Budiansky, B. (1974). “Theory of buckling and post-buckling behavior of elastic structures,” in Advances in Applied Mechanics, eds Chia-Shun Yih (New York, NY: Academic Press), 1–65.

Google Scholar

Feo, L., and Fraternali, F. (2000). On a moderate rotation theory of thin-walled composite beams. Compos. B. Eng. 31, 141–158. doi: 10.1016/S1359-8368(99)00079-7

CrossRef Full Text | Google Scholar

Feo, L., and Mancusi, G. (2010). Modeling shear deformability of thin-walled composite beams with open cross-section. Mech. Res. Commun. 37, 320–325. doi: 10.1016/j.mechrescom.2010.02.005

CrossRef Full Text | Google Scholar

Fraternali, F. (2007). Free discontinuity finite element models in two-dimensions for in-plane crack problems. Theor. Appl. Fract. Mech. 47, 274–282. doi: 10.1016/j.tafmec.2007.01.006

CrossRef Full Text | Google Scholar

Fraternali, F., Angelillo, M., and Fortunato, A. (2002). A lumped stress method for plane elastic problems and the discrete-continuum approximation. Int. J. Solids Struct. 39, 6211–6240. doi: 10.1016/S0020-7683(02)00472-9

CrossRef Full Text | Google Scholar

Fraternali, F., Ciancia, V., Chechile, R., Rizzano, G., and Feo, L. (2011). Incarnato, L. Experimental study of the thermo-mechanical properties of recycled PET fiber reinforced concrete. Comp. Struct. 93, 2368–2374. doi: 10.1016/j.compstruct.2011.03.025

CrossRef Full Text | Google Scholar

Fraternali, F., Farina, I., Polzone, C., Pagliuca, E., and Feo, L. (2012). On the use of R-PET strips for the reinforcement of cement mortars. Compos. B. Eng. 46, 207–210. doi: 10.1016/j.compositesb.2012.09.070

CrossRef Full Text | Google Scholar

Fraternali, F., Negri, M., and Ortiz, M. (2010). On the convergence of 3D free discontinuity models in variational fracture. Int. J. Fract. 166, 3–11. doi: 10.1007/s10704-010-9462-0

CrossRef Full Text | Google Scholar

Fraternali, F., Spadea, S., and Ascione, L. (2013). Buckling behavior of curved composite beams with different elastic response in tension and compression. Compos. Struct. 100, 280–289, doi: 10.1016/j.compstruct.2012.12.021

CrossRef Full Text | Google Scholar

Keller, H. B. (1977). “Numerical solution of bifurcation and nonlinear eigenvalue problems,” in Applications of Bifurcation Theory, eds P.H. Rabinowitz (New York, NY: Academic Press), 359–384.

Google Scholar

Markkula, S., Malecki, H. C., and Zupan, M. (2013). Uniaxial tension and compression characterization of hybrid CNS-glass fiber-epoxy composites. Comp. Struct. 95, 337–345. doi: 10.1016/j.compstruct.2012.07.013

CrossRef Full Text | Google Scholar

Mascolo, I., Fulgione, M., and Pasquino, M (in press). Lateral torsional buckling of compressed open thin walled beams: experimental confermations. Int. J. Mansory Res. Innovat.

Mascolo, I., and Pasquino, M. (2016). Lateral-torsional buckling of compressed highly variable cross section beams. Curved. Al. Layer. Struct. 3443, 146–153. doi: 10.1515/cls-2016-0012

CrossRef Full Text | Google Scholar

Nomizu, K., and Kobayashi, S. (1963). Foundations of Differential Geometry. New York, NY: Wiley.

Google Scholar

Özütok, A., and Madenci, E. (2017). Static analysis of laminated composite beams based on higher-order shear deformation theory by using mixed-type finite element method. Int. J. Mech. Sci. 130, 234–243. doi: 10.1016/j.ijmecsci.2017.06.013

CrossRef Full Text | Google Scholar

Prathap, G., and Bhashyam, G. (1982). Reduced integration and the shear-flexible beam element. lnt.J. Num. Meth. Engng. 18, 195–210. doi: 10.1002/nme.1620180205

CrossRef Full Text | Google Scholar

Reddy, J. N. (1992). An Introduction to the Finite Element Method, 2nd Edn. New York, NY: McGraw-Hill.

Google Scholar

Reddy, J. N., and Averill, R. C. (1990). On the behavior of plate elements based on the first-order shear deformation theory. Engng. Comput. 7, 57–74. doi: 10.1108/eb023794

CrossRef Full Text | Google Scholar

Riks, E. (1972). The application of newton's method to the problem of elastic stability. J. Appl. Mech. 39, 1060–1066. doi: 10.1115/1.3422829

CrossRef Full Text | Google Scholar

Roberts, T. M., and Al-Ubaidi, H. (2001). Influence of shear deformation on restrained torsional warping of pultruded FRP bars of open cross-section. Thin-Walled Struct. 39, 395–414. doi: 10.1016/S0263-8231(01)00009-X

CrossRef Full Text | Google Scholar

Schmidt, B., Fraternali, F., and Ortiz, M. (2009). Eigenfracture: an eigendeformation approach to variational fracture. Siam. J. Multiscale Model. Simul. 7, 1237–1266. doi: 10.1137/080712568

CrossRef Full Text | Google Scholar

Simo, J. C., and Wriggers, P. (1990). A general procedure for the direct computation of turning and bifurcation points. int. J. Num. Meth. Engng. 30, 155–176. doi: 10.1002/nme.1620300110

CrossRef Full Text | Google Scholar

Spence, A., and Jepson, A. D. (1985). Folds in the solution of two parameter systems and their calculation, Pan. I, SIAM J. Num. Anal. 22, 347–368. doi: 10.1137/0722021

CrossRef Full Text | Google Scholar

Timoshenko, S. P., and Gere, J. M. (1961). Theory of Elastic Stability. New York, NY: McGraw-Hill.

Google Scholar

Trefftz, E. (1930). “Über die Ableitung der Stäbilitats-kriterien des elastischen Gleichgewichtes aus der elasticitäts theorie endlicher Deformationen,” in Proceedings of the Third international Congress for Applied Mechanics (Stockholm), 44–50.

Viera, R. F., Virtuoso, F. B. E., and Pereira, E. B. R. (2013). A higher order thin-walled beam model including warping and shear modes. Int. J. Mech. Sci. 66, 67–82. doi: 10.1016/j.ijmecsci.2012.10.009

CrossRef Full Text | Google Scholar


We report below the list of the main notations used in the text. Throughout this paper we use boldface character to denote numerical vectors and matrices; we also use the superscript T to denote the transpose of a vector or a matrix.

Given a scalar function f(U,λ):N×ℜ → ℜ, we denote by DUf the vector which represents the Gateaux derivative of f with respect to U

VTDUF=limγ01γ[f(U+γ V,λ)f(U,λ)],VN

and by Dλf the partial derivative of f with respect to λ.

Similarly, given a vector function R(U,λ):N×ℜ → ℜN, we denote by DUR the N × N matrix which represents the Gateaux derivative of R with respect to U

DUR V=limγ01γ[R(U+γ V,λ)R(U,λ)],VN,

and by DλR the partial derivative of R with respect to λ.

Finally, in the case of a scalar function f(U, λ), we denote by DU2f the bilinear operator DU2f V1V2=V2TDU(DUf) V1,V1,V2N.

Keywords: composite laminates, curved beams, fiber-reinforced composites, buckling-analysis, finite element method

Citation: Mascolo I, Modano M, Amendola A and Fraternali F (2018) A Finite Element Analysis of the Stability of Composite Beams With Arbitrary Curvature. Front. Built Environ. 4:57. doi: 10.3389/fbuil.2018.00057

Received: 28 July 2018; Accepted: 28 September 2018;
Published: 23 October 2018.

Edited by:

Vagelis Plevris, OsloMet–Oslo Metropolitan University, Norway

Reviewed by:

Georgios A. Drosopoulos, University of KwaZulu-Natal, South Africa
Francesco Tornabene, Università degli Studi di Bologna, Italy

Copyright © 2018 Mascolo, Modano, Amendola and Fraternali. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Fernando Fraternali,