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

Abstract

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₁, …, C_n belong to the Lagrange family (Reddy, 1992) (see Figure 1)

Figure 1

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₁, Z₂, Z₃} (Figure 1); 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

The Jacobian of the transformation from the local coordinate ξ (Figure 2) to the global coordinate X₃ (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₃, 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)₃, 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)₁ are given by

where

The matrices A and G, which appear in Equation (12)₂, 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(U, λ) (residuai vector) is the Gateaux derivative of Π^h with respect to .

On accounting for a possibility of non-linear elastic response of the material, we assume that the elasticity matrix , 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 ( = ).

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 is the generalized stress vector

while Q⁽¹⁾ and Q⁽²⁾(U) are the global force vectors, which derive from the assembly of the element vectors

and the nodal force vectors and .

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_UR, which is usually referred to as tangent stiffness matrix and denoted by K_T, is given in the Appendix. Concerning the vector D_λR, from Equation (23) we deduce

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 , , 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₁ is usually called stability point, while U₁ 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₁ 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

V₀ 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 V becoming exceedingly large. Concerning V₀, we set

V₀ 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⁽¹⁾ + Q⁽²⁾(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₀ (see the Appendix) and an arbitrary numeric value by λ₀, we consider the couple λ₀, 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 , , , (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₁ and H₂ along the directions X₁ and X₂, respectively. We denoted the cross-section area by A = H₁H₂, the moments of inertia by I₁ and I₂, 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₁₁X₁X₂

Warping Function (W3):

.

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₁/H₂ = 2, R/H₂ = 20 were considered.

According to Donnel' s theory (see e.g., Brush and Almroth, 1975), the first bifurcation point occurs at a load level 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 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 η.

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₀, 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 .

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₁ and X₂ of the axial strain E₃₃, 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₃ = L). Table 4 compares the bifurcation loads computed by using the present theory with those corresponding to Euler's theory () and Timoshenko's theory, for several the values of the ratio L/H₂, assuming H₁/H₂ = 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₁, X₂ of the axial strain E₃₃ (Γ 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 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₁/H₂ = 0.1, L/H₂ = 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₁ = 0, X₂ = − H₂/2) and load at the intrados (X₁ = 0, X₂ = H₂/2). The last two cases obviously give rise to a deformation-dependent loading (Q⁽²⁾(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).

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 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₁ = v₂ = v₃ = ϕ₃ = 0 for X₃ = 0, πR in the first case). The load consists of a uniform radial dead load along the centerline.

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:

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

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

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_Uf the vector which represents the Gateaux derivative of f with respect to U

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_UR the N × N matrix which represents the Gateaux derivative of R with respect to U

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 the bilinear operator

References

• 1

  AscioneL.BerardiV. P.GiordanoA.SpadeaS. (2013). Local buckling behavior of FRP thin-walled beams: a mechanical model. Comp Struct. 98, 111–120. 10.1016/j.compstruct.2012.10.049

  • CrossRef
  • Google Scholar

• 2

  AscioneL.GiordanoA.SpadeaS. (2011). Lateral buckling of pultruded FRP beams. Compos B Eng. 42, 819–824. 10.1016/j.compositesb.2011.01.015

  • CrossRef
  • Google Scholar

• 3

  BathozI. L.DhattG. (1979). Incremental displacement algorithms for non-linear problems. Int. J. Num. Meth. Engng. 14, 1262–1267. 10.1002/nme.1620140811

  • CrossRef
  • Google Scholar

• 4

  BencardinoF.RizzutiL.SpadeaG.SwamyR. N. (2012). Implications of test methodology on post-cracking and fracture behavior of steel fibre reinforced concrete. Compos. B Eng. 46, 31–38. 10.1016/j.compositesb.2012.10.016

  • CrossRef
  • Google Scholar

• 5

  BrushD. O.AlmrothB. O. (1975). Buckling of Bars, Plates and Shells. New York, NY: McGraw-Hill.

  • Google Scholar

• 6

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

  • Google Scholar

• 7

  FeoL.FraternaliF. (2000). On a moderate rotation theory of thin-walled composite beams. Compos. B. Eng. 31, 141–158. 10.1016/S1359-8368(99)00079-7

  • CrossRef
  • Google Scholar

• 8

  FeoL.MancusiG. (2010). Modeling shear deformability of thin-walled composite beams with open cross-section. Mech. Res. Commun.37, 320–325. 10.1016/j.mechrescom.2010.02.005

  • CrossRef
  • Google Scholar

• 9

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

  • CrossRef
  • Google Scholar

• 10

  FraternaliF.AngelilloM.FortunatoA. (2002). A lumped stress method for plane elastic problems and the discrete-continuum approximation. Int. J. Solids Struct. 39, 6211–6240. 10.1016/S0020-7683(02)00472-9

  • CrossRef
  • Google Scholar

• 11

  FraternaliF.CianciaV.ChechileR.RizzanoG.FeoL. (2011). Incarnato, L. Experimental study of the thermo-mechanical properties of recycled PET fiber reinforced concrete. Comp. Struct. 93, 2368–2374. 10.1016/j.compstruct.2011.03.025

  • CrossRef
  • Google Scholar

• 12

  FraternaliF.FarinaI.PolzoneC.PagliucaE.FeoL. (2012). On the use of R-PET strips for the reinforcement of cement mortars. Compos. B. Eng. 46, 207–210. 10.1016/j.compositesb.2012.09.070

  • CrossRef
  • Google Scholar

• 13

  FraternaliF.NegriM.OrtizM. (2010). On the convergence of 3D free discontinuity models in variational fracture. Int. J. Fract.166, 3–11. 10.1007/s10704-010-9462-0

  • CrossRef
  • Google Scholar

• 14

  FraternaliF.SpadeaS.AscioneL. (2013). Buckling behavior of curved composite beams with different elastic response in tension and compression. Compos. Struct. 100, 280–289, 10.1016/j.compstruct.2012.12.021

  • CrossRef
  • Google Scholar

• 15

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

  • Google Scholar

• 16

  MarkkulaS.MaleckiH. C.ZupanM. (2013). Uniaxial tension and compression characterization of hybrid CNS-glass fiber-epoxy composites. Comp. Struct. 95, 337–345. 10.1016/j.compstruct.2012.07.013

  • CrossRef
  • Google Scholar

• 17

  MascoloI.FulgioneM.PasquinoM. (in press). Lateral torsional buckling of compressed open thin walled beams: experimental confermations. Int. J. Mansory Res. Innovat.

  • Google Scholar

• 18

  MascoloI.PasquinoM. (2016). Lateral-torsional buckling of compressed highly variable cross section beams. Curved. Al. Layer. Struct. 3443, 146–153. 10.1515/cls-2016-0012

  • CrossRef
  • Google Scholar

• 19

  NomizuK.KobayashiS. (1963). Foundations of Differential Geometry. New York, NY: Wiley.

  • Google Scholar

• 20

  ÖzütokA.MadenciE. (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. 10.1016/j.ijmecsci.2017.06.013

  • CrossRef
  • Google Scholar

• 21

  PrathapG.BhashyamG. (1982). Reduced integration and the shear-flexible beam element. lnt.J. Num. Meth. Engng. 18, 195–210. 10.1002/nme.1620180205

  • CrossRef
  • Google Scholar

• 22

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

  • Google Scholar

• 23

  ReddyJ. N.AverillR. C. (1990). On the behavior of plate elements based on the first-order shear deformation theory. Engng. Comput. 7, 57–74. 10.1108/eb023794

  • CrossRef
  • Google Scholar

• 24

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

  • CrossRef
  • Google Scholar

• 25

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

  • CrossRef
  • Google Scholar

• 26

  SchmidtB.FraternaliF.OrtizM. (2009). Eigenfracture: an eigendeformation approach to variational fracture. Siam. J. Multiscale Model. Simul. 7, 1237–1266. 10.1137/080712568

  • CrossRef
  • Google Scholar

• 27

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

  • CrossRef
  • Google Scholar

• 28

  SpenceA.JepsonA. D. (1985). Folds in the solution of two parameter systems and their calculation, Pan. I, SIAMJ. Num. Anal. 22, 347–368. 10.1137/0722021

  • CrossRef
  • Google Scholar

• 29

  TimoshenkoS. P.GereJ. M. (1961). Theory of Elastic Stability. New York, NY: McGraw-Hill.

  • Google Scholar

• 30

  TrefftzE. (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.

  • Google Scholar

• 31

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

  • CrossRef
  • Google Scholar