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

^{1}Department of Civil Engineering, University of Salerno, Fisciano, Italy^{2}Department 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.

## Introduction

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 *C*^{h} a finite-element discretization of axis curve of a laminated beam, and let us assume that the elements *C*_{1}, …, *C*_{n} belong to the Lagrange family (Reddy, 1992) (see Figure 1)

**Figure 1**. Finite element discretization of the axis of a laminate curved beam (*n*_{l}: 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 *C*_{e}

In Equation (2) *N*_{I} is the shape function corresponding to the node *I* and consists of a complete polynomial of order *n* − 1; *Z*_{2I} and *Z*_{3I} are the coordinates of *I* with respect to the global frame {0, *Z*_{1}, *Z*_{2}, *Z*_{3}} (Figure 1); ${\widehat{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.

The Jacobian of the transformation from the local coordinate ξ (Figure 2) to the global coordinate *X*_{3} (Figure 1) is given by Nomizu and Kobayashi (1963)

*N*_{I, ξ} being the derivative of *N*_{I} with respect to ξ. Denoting by (·)′ the derivative with respect to *X*_{3}, we have:

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

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

where *U*_{e} is the M-dimensional vector collecting nodal, generalized, displacements of

while ** N** is the following matrix

whose blocks are diagonal submatrices

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

where

The σ × *m* submatrices *B*_{0I} in Equation (12)_{1} are given by

where

The matrices ** A** and

**, which appear in Equation (12)**

*G*_{2}, are instead given by

where

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

where Π^{h} is the discretized functional of the total potential energy Π; δΠ^{h} is the first variation of Π^{h} with increment δ** U**;

**(**

*R***, λ) (residuai vector) is the Gateaux derivative of Π**

*U*^{h}with respect to $U(R={D}_{U}{\Pi}^{h})$.

On accounting for a possibility of non-linear elastic response of the material, we assume that the elasticity matrix $\widehat{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 ($\widehat{\text{D}}$ = $\widehat{\text{D}}(U)$).

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

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

where $\widehat{S}$ is the generalized stress vector

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

and the nodal force vectors ${\widehat{Q}}_{I}^{(1)}$ and ${\widehat{Q}}_{I}^{(2)}\text{}({\widehat{u}}_{I})\text{}(I=1,\text{}2,\text{}\dots ,\text{}{n}_{n})$.

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

*N*equations

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

where *e*_{p} 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

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

The matrix *D*_{U}** R**, which is usually referred to as tangent stiffness matrix and denoted by

*K*_{T}, is given in the Appendix. Concerning the vector

*D*

_{λ}

**, from Equation (23) we deduce**

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

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.

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

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 *U*_{p} which exhibited the largest variation in the previous step. We hence set in the extended system (25)

and begin the new iteration loop by assuming the predictor $\stackrel{~}{U}={\text{}U}^{i-1}+\delta {e}_{p}$, $\stackrel{~}{\lambda}={\lambda}^{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

which must be satisfied by every variation δ** U**.

A point ** U**, λ such that Equation (36) holds for some

*U*_{1}is usually called stability point, while

*U*_{1}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

In the following we will denote *U*_{1} by ** V**. In order to exclude the trivial case

**, it is necessary to append a constraint equation**

*V*=0*l*(

**)**

*V***=**0 to system (37). Possible choices of such an equation are

*V*_{0} 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

*K*_{T}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

**becoming exceedingly large. Concerning**

*V**V*

_{0}, we set

*V*_{0} 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

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

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 *K*_{T} 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)

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

where *Q*=*Q*^{(1)} + *Q*^{(2)}(** U**), while

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

We furnish the expressions of the vectors *h*_{j} (*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 *K*_{0} (see the Appendix) and an arbitrary numeric value by λ_{0}, we consider the couple λ_{0}, ${U}_{0}={\lambda}_{0}{K}_{0}^{-1}{Q}^{(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 *K*_{T}.

If the sign of det *K*_{T} 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 $\stackrel{~}{U}={\text{U}}^{i+1}$, $\stackrel{~}{V}={V}_{0}={K}_{0}^{-1}{e}_{p}$, $\stackrel{~}{\lambda}={\lambda}^{i+1}$, $\stackrel{~}{\mu}={{e}_{p}^{T}\text{U}}^{i+1}$ (for the meaning of the index p see the beginning of section Computation of Stability Points).

Once the stability point *U*_{c}, λ_{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 *U*_{c} a vector proportional to the eigenvector *V*

ζ being a scaling factor to be determined in such a way that it results ** R**(

**,λ**

*U*_{c})

**≤**

*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 *H*_{1} and *H*_{2} along the directions *X*_{1} and *X*_{2}, respectively. We denoted the cross-section area by *A* = *H*_{1}*H*_{2}, the moments of inertia by *I*_{1} and *I*_{2}, the polar moment of inertia by *I*_{G} and the De Saint Venant torsional rigidity by *J*_{t}

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

No Warping (NW): *w* = 0

Warping Function (W1): *w* = *w*_{11}*X*_{1}*X*_{2}

Warping Function (W3): $w={w}_{20}{X}_{1}^{2}+{w}_{11}{X}_{1}{X}_{2}+{w}_{02}{X}_{2}^{2}+$

${w}_{30}{X}_{1}^{3}+{w}_{21}{X}_{1}^{2}{X}_{2}+{w}_{12}{{X}_{1}X}_{2}^{3}+{w}_{03}{X}_{2}^{3}$.

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)

where *R* is the initial radius of the centerline. In particular, the ratios *H*_{1}/*H*_{2} = 2, *R*/*H*_{2} = 20 were considered.

According to Donnel' s theory (see e.g., Brush and Almroth, 1975), the first bifurcation point occurs at a load level ${q}_{bif}=\text{}4\text{}E\text{}{I}_{1}/{R}^{3}$ 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 ${q}_{bif}=\text{}3\text{}E\text{}{I}_{1}/{R}^{3}$ 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 (*H*_{1}/*H*_{2} = 2, *R*/*H*_{2} = 20, 20 element mesh).

Denoting by *D*_{p} the least diagonal term of the factorized tangent stiffness matrix and by *D*_{op} the corresponding term in the factorized initial stiffness matrix *K*_{0}, in our numerical experiments we set η = 0, η = (*D*_{op} − *D*_{p}) × 10 and η = (*D*_{op} − *D*_{p}) × 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\text{}E{I}_{1}/{R}^{3}$.

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 *X*_{1} and *X*_{2} of the axial strain *E*_{33}, which are instead present in our model. In all the subsequent numerical applications we η = (*D*_{op} − *D*_{p}) × 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 (*X*_{3} = *L*). Table 4 compares the bifurcation loads computed by using the present theory with those corresponding to Euler's theory (${P}_{Eul}={\pi}^{2}\text{}E{I}_{1}/{L}^{2}$) and Timoshenko's theory, for several the values of the ratio *L*/*H*_{2}, assuming *H*_{1}/*H*_{2} = 2. Upon putting ρ = χ* P*_{Eul}/*GA* (χ = 1.2 shear correction factor), we considered both the exact values deriving from Timoshenko's Theory

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

Concerning the comparisons between our and classical theories, we stressed the influences of warping and quadratic terms in *X*_{1}, *X*_{2} of the axial strain *E*_{33} (Γ 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*/*H*_{2} (*H*_{1}/*H*_{2} = 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 (*H*_{1}/*H*_{2} = 0.1, *L*/*H*_{2} = 10). The kinematical boundary conditions of the first case are

Three positions of the load were considered: load at the centroid, load at the extrados (*X*_{1} = 0, *X*_{2} = − *H*_{2}/2) and load at the intrados (*X*_{1} = 0, *X*_{2} = *H*_{2}/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 (*H*_{1}/*H*_{2} = 0.1, *L*/*H*_{2} = 10, 20 element mesh).

In the context of the present theory, we adopted a bilinear warping function (*W*1) and assumed both a linear and a (geometrically) nonlinear pre-buckling behavior (the first was treated by discarding the part *K*_{L} 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 *GJ*_{t}.

It has to be remarked that the present theory assumes a torsional rigidity equal to *GI*_{G} (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 (*v*_{1} = *v*_{2} = *v*_{3} = ϕ_{3} = 0 for *X*_{3} = 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 (*H*_{1}/*H*_{2} = 0.1, *R*/*H*_{2} = 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.

## Acknowledgments

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: https://www.frontiersin.org/articles/10.3389/fbuil.2018.00057/full#supplementary-material

## References

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

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

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

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

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.

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

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

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

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

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

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

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

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

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.

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

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

Ö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

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

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

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

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

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

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

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

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

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

## Notation

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

*D*

_{U}

*f*the vector which represents the Gateaux derivative of

*f*with respect to

*U*${V}^{T}{D}_{U}F={\text{lim}}_{\gamma \to 0}\frac{1}{\gamma}\left[f\left(U+\gamma \text{}V,\lambda \right)-f\left(U,\lambda \right)\right],\forall V\in {\Re}^{N}$

and by *D*_{λ}*f* the partial derivative of *f* with respect to λ.

Similarly, given a vector function ** R**(

**,λ)**

*U***:**ℜ

^{N}

**×**ℜ → ℜ

^{N}, we denote by

*D*

_{U}

**the**

*R**N*×

*N*matrix which represents the Gateaux derivative of

**with respect to**

*R*

*U*${D}_{U}R\text{}V={\mathrm{lim}}_{\gamma \to 0}\frac{1}{\gamma}\left[R\left(U+\gamma \text{}V,\lambda \right)-R\left(U,\lambda \right)\right],\forall V\in {\Re}^{N},$

and by *D*_{λ}** R** the partial derivative of

**with respect to λ.**

*R*Finally, in the case of a scalar function *f*(** U**, λ), we denote by ${D}_{U}^{2}f$ the bilinear operator ${D}_{U}^{2}f\text{}{V}_{1}{V}_{2}={V}_{2}^{T}{D}_{U}\left({D}_{U}f\right)\text{}{V}_{1},\forall {V}_{1},{V}_{2}\in {\Re}^{N}.$

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, NorwayReviewed by:

Georgios A. Drosopoulos, University of KwaZulu-Natal, South AfricaFrancesco 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, f.fraternali@unisa.it