Nonlinear bifurcation analysis of statically loaded free-form curved beams using isogeometric framework and pseudo-arclength continuation

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Nonlinear bifurcation analysis of statically loaded free-form curved beams using isogeometric framework and pseudo-arclength continuation Ali Hashemian a ,∗, Seyed Farhad Hosseini b a b

Department of Mechanical Engineering, Hakim Sabzevari University, Sabzevar, Iran Sun-Air Research Institute, Ferdowsi University of Mashhad, Mashhad, Iran

ARTICLE

INFO

Keywords: Isogeometric analysis Geometrically nonlinear curved beams Static bifurcation Pseudo-arclength continuation Timoshenko beam theory

ABSTRACT This paper presents an isogeometric analysis (IGA) integrated with pseudo-arclength continuation technique in order to predict the elastic bifurcation behavior of geometrically nonlinear free-form curved beams under static loading. The presented formulation is based on the Timoshenko beam theory and accounts for shear deformations in the large-displacement analysis of curved beams. The full equilibrium path of the buckled structures including common types of static instabilities (e.g., the saddle–node and pitchfork bifurcations) can be tracked by the presented approach. Thanks to the exact geometry representation of the IGA framework, the objective of this paper is to extend the current literature of the nonlinear buckling analysis of curved beams to a broader range of free-form geometries with engineering applications. Some practical case studies with different loading scenarios are investigated and the results are validated against the commercial finite element software and well-established benchmark examples.

1. Introduction Curved beams are among the most commonly used structural members that have a wide variety of applications in aerospace, civil and mechanical designs. In contrast to straight beams that generally have uncoupled longitudinal and transverse displacement, the geometric curvature in curved beams makes the displacements coupled. In addition, when the applied transverse load on the beam becomes large, the beam develops internal reaction forces that affect its deformation, and the linear load–deflection relationship will not be valid anymore [1]. As a result, the beam stiffness will itself be a function of displacement field, which is commonly referred to as geometric nonlinearity [2–5]. In the cases of softening-type nonlinearity where increasing the applied load decreases the stiffness of the beam, the curved beam is prone to lose its stability and suffer a static elastic bifurcation (i.e. buckling). Two common types of static bifurcation which are recurrent in a rich variety of curved beams are saddle–node and pitchfork bifurcations [5,6]. In the saddle–node (or fold) bifurcation, there is no stable equilibrium state past the bifurcation point (which is generally referred to as limit or turning point). This type of instability can be observed in shallow curved beams in the form of snap-through or snap-back behavior when the system is statically loaded in a force or displacement-controlled manner. On the other hand, in the pitchfork bifurcation, when the fundamental equilibrium state of the structure suffers a loss of stability at a critical point, the system bifurcates smoothly into either one of the

two stable equilibrium states. This type of instability can be observed in non-shallow curved structures with imperfections. Studying the structural response of buckled curved beams is of interest to many researchers in recent years. Some research works like, e.g., those performed by Tadjbakhsh [7], Attard et al. [8] and Luu et al. [9] focused on the linear buckling analysis of parabolic beams. Including the geometric nonlinearities due to large deformations in the analysis, Lo [10], Addessi et al. [11], Pi and Pi and Bradford [12, 13], Kim et al. [14,15] and Zhu et al. [16] investigated the nonlinear buckling phenomenon in the curved beams. However, their researches are limited to constant curvature circular arches. For non-circular curved beams, Nieh et al. [17] and Luu and Luu and Lee [18] studied the static bifurcation in elliptic beams subject to a concentrated load at the middle, where the latter work employs the Timoshenko beam theory to account for shear deformations in the analysis. For parabolic beams, we can refer to researches by Cai et al. [19] and Zhu et al. [20] that explored different loading scenarios in the analysis. In some other researches, like e.g., [21,22] the nonlinear vibration analysis of beams are also studied. A summary of the above-mentioned literature is represented in Table 1. It can be inferred that most of research works in this area are focusing on geometries of conic sections for which analytical solutions are derived in some occasions. However, based on the geometrical aspects of modern engineering and architectural design, extending the domain of static bifurcation analysis to a broader range of free-form curved beams with variable curvature is increasingly

∗ Corresponding author. E-mail address: [email protected] (A. Hashemian).

https://doi.org/10.1016/j.ijnonlinmec.2019.03.002 Received 16 October 2018; Received in revised form 20 February 2019; Accepted 3 March 2019 Available online 26 March 2019 0020-7462/© 2019 Elsevier Ltd. All rights reserved.

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Table 1 Some selected research works concerning the buckling analysis of statically loaded curved beams. Research by

Geometric nonlinearity

Variable curvature

Geometry type

Shear deformable theory

Different loading scenarios

Solution method

Remarks

Tadjbakhsh [7]

—-

✔

Parabolic

—-

✔

Perturbation

Attard et al. [8]

—-

✔

Parabolic

✔

✔

Weighted residual method (WRM)

Luu et al. [9]

—-

✔

Parabolic

✔

✔

IGA

Lo [10]

✔

—-

Circular

✔

✔

FEA integrated with continuation

Addessi et al. [11]

✔

—-

Circular

✔

—-

Analytical

Pi and Bradford [12,13]

✔

—-

Circular

—-

—-

Analytical

Kim et al. [14]

✔

—-

Circular

✔

✔

Analytical

Kim et al. [15]

✔

—-

Circular

✔

✔

FEA

Zhu et al. [16]

✔

—-

Circular

✔

✔

Nieh et al. [17]

—-

✔

Elliptic

—-

—-

Trapezoid with Richardson extrapolation Analytical

Luu and Lee [18]

✔

✔

Elliptic

✔

—-

IGA integrated with continuation

Cai et al. [19]

✔

✔

Parabolic

—-

—-

Analytical

Zhu et al. [20]

✔

✔

Parabolic

✔

—-

Trapezoid with Richardson extrapolation

Limited to inextensible beam theory Simple basis functions of WRM are not suitable for free-form beams Geometric nonlinearity is not included Requiring a fine mesh and high number of iterations/increments to converge Only hinged–hinged circular arches under a concentrated load are studied Only hinged–hinged circular arches under a concentrated load are studied Limited to circular arches with constant curvature Limited to circular arches with constant curvature Limited to circular arches with constant curvature Geometric nonlinearity is not included Limited to elliptic arches under concentrated loading Limited to parabolic arches under distributed loading Limited to parabolic arches under distributed loading

essential. Since analytical approaches can rarely find a closed-form solution for such geometries, some methods are based on developing the finite element analysis (FEA) that, on the other hand, needs a fine mesh and high number of iterations/increments to converge. This is generally because FEA approach discretizes the geometry into lower order and zero continuity elements. In order to resolve this difficulty, the isogeometric analysis (IGA), introduced by Hughes et al. [23], with the advantages like, e.g., higher order elements with arbitrary continuity, higher accuracies and lower computational effort can be employed. Isogeometric analysis has recently proved itself as an applicable model for bridging the gap between computer-aided design (CAD) and finite element analysis. The main characteristic of the IGA framework is that the shape functions of the solution space in FEA are considered the same as the basis functions defining the CAD geometry. Thanks to exact geometry representation by non-uniform rational Bsplines (NURBS) along with the mesh smoothing capabilities, the IGA framework has a superior computational performance compared to the conventional FEA. This method has been successfully implemented to a wide range of engineering problems with applications in, e.g., structural mechanics [24–28], fluid mechanics [29–31], heat transfer [32– 34] and eigenvalue problems [35,36]. Isogeometric analysis of freeform curved beams has also attracted many researchers. Some research works in this area can be revisited here as follows. Bouclier et al. [37] and Zhang et al. [38] studied the locking-free isogeometric analysis of Timoshenko curved beams. Isogeometric shape optimization of the

curved beams are studied by Hosseini et al. [39] and Nagy et al. [40]. Focusing on large deformation analysis, Bauer et al. [41] and Hosseini et al. [42], by presenting a nonlinear iterative solution, considered the geometric nonlinearity in the analysis of Euler–Bernoulli and Timoshenko curved beams, respectively. Some other research works related to the isogeometric analysis of curved beams can be found in [43–47]. It is also interesting to note that, in recent years, isogeometric analysis has been successfully applied to buckling analysis of different types of structures (see, e.g., [9,18,48–51]). In order to follow the equilibrium path completely in nonlinear problems, the path-following or so-called continuation techniques are frequently employed by analysts (see, e.g., [52–55]). These methods have been originally introduced by Riks [56,57] and Wempner [58], and then modified by other researchers (e.g., Crisfield [59] and Feng et al. [60,61]) to cover all practical aspects of the solution. Among different continuation schemes in the literature, in this paper, the pseudo-arclength continuation technique integrated with the isogeometric framework is presented in order to predict the static elastic bifurcation behavior in geometrically nonlinear free-form curved beams. The presented IGA formulation is based on the Timoshenko beam theory including shear deformations in the large-displacement analysis of curved beams. Considering the importance of free-form geometries in the modern engineering applications, the main contribution of this research lies in implementation of IGA in extending the current literature of the nonlinear buckling analysis of curved beams to a larger range 2

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 2. DOFs of an arbitrary point on the midline of a free-form Timoshenko curved beam.

Fig. 1. Kinematics of a curved beam in 𝑠-𝑧 coordinate system based on the Timoshenko beam theory.

of practical free-form geometries with variable curvature. Additionally, various beam geometries under different loading conditions are studied in this paper and validated against commercial FEA software and well-established benchmark examples. The remainder of this article is organized as follows. Section 2 presents the isogeometric formulation of geometrically nonlinear curved beams based on the Timoshenko beam theory. In Section 3, the pseudo-arclength continuation technique is introduced based on the IGA formulation in order to follow the full equilibrium path of the buckled beams. Section 4 deals with the case studies, followed by Section 5 that draws some conclusions.

Fig. 3. A cubic B-spline curve with eight control points.

Table 2 B-spline representation of the geometry and displacement fields of curved beams. Geometry 𝑥(𝜉) =

𝑛 ∑

Displacement fields over the 𝑘th element

𝑦 (𝜉) =

∑

∑

𝑘+𝑝−1

𝑁𝑖,𝑝 (𝜉) 𝑋𝑖

𝑢𝑒 (𝜉) =

𝑖=0 𝑛

𝑁𝑖,𝑝 (𝜉) 𝑈𝑖

𝑖=𝑘−1 𝑘+𝑝−1

𝑁𝑖,𝑝 (𝜉) 𝑌𝑖

𝑤𝑒 (𝜉) =

𝑖=0

∑

𝑁𝑖,𝑝 (𝜉) 𝑊𝑖

𝑖=𝑘−1 𝑘+𝑝−1

𝜑𝑒 (𝜉) =

2. Isogeometric analysis of geometrically nonlinear curved beams

∑

𝑁𝑖,𝑝 (𝜉) Φ𝑖

𝑖=𝑘−1

2.1. Kinematics of the problem as follows [42,63,64] where 𝜀, 𝜒 and 𝛾 denote the membrane, flexural and shear strains, respectively, and 𝜅 = 1∕𝜌 is the curvature of the beam. Note that to write the equations more compactly, here and in the following, the superscript ′ refers to differentiation with respect to 𝑠. 1 1 𝜀 = 𝑢′ − 𝜅𝑤 + 𝑢′2 + 𝑤′2 2 2 (2) 𝜒 = 𝜑′

The present formulation of geometrically nonlinear curved beams is in accordance with Timoshenko beam theory, allowing for shear deformation. It implies that plane sections normal to the midline remain plane after deformation, while they are not required to remain normal to the midline. In addition, the formulation is based on the assumptions of large displacements, small strains and small to moderate rotations [1]. The displacement field of a curved beam with variable curvature in the curvilinear coordinate system (𝑠-𝑧) is demonstrated in Fig. 1 where the 𝑠-curvilinear coordinate coincides with the midline of the beam, 𝜌 is the radius of curvature, 𝑢 and 𝑤 are tangential and normal displacements of the midline, respectively, 𝜓 is the rotation of the transverse straight line on the cross section, and 𝑤′ = 𝑑𝑤∕𝑑𝑠. Before proceeding with the displacement field, it is important to mention that the rotation due to flexure (𝜓) cannot be considered as a nodal degree of freedom (DOF) for curved beam elements as it is not equal to the physical rotation of the section with respect to the global coordinate system. The appropriate choice of nodal DOF, instead, would be 𝜑 = 𝜓 + 𝑢∕𝜌 for Timoshenko curved beams [62]. Hence, the displacement fields of an arbitrary point on a free-form curved beam in the curvilinear coordinate system (Fig. 2) are defined as:

𝛾 = 𝜅𝑢 + 𝑤′ − 𝜑 2.2. Governing equations The total potential energy (Π) of a curved beam, as an elastic member, is defined as follows where U and W are the elastic strain energy and work of external loads, respectively. Π=U−W 𝑙

1 2 ∫0 𝑙( ) W= 𝑢𝑓𝑠 + 𝑤𝑓𝑧 + 𝜑𝑚 𝑑𝑠 ∫0 U=

𝑢 (𝑠, 𝑧) = 𝑢 (𝑠) − 𝑧𝜑(𝑠) 𝑤 (𝑠, 𝑧) = 𝑤 (𝑠)

( ) 𝐸𝐴𝜀2 + 𝐸𝐼𝜒 2 + 𝐾𝑠 𝐺𝐴𝛾 2 𝑑𝑠

(3) (4) (5)

In the above equations, 𝑓𝑠 , 𝑓𝑧 and 𝑚 are the external distributed longitudinal, transverse and moment loads, respectively. The material and cross-sectional properties of the beam are characterized by 𝐸, 𝐺, 𝐴 and 𝐼, while 𝑙 is the length of the beam and 𝐾𝑠 is the shear correction

(1)

The nonlinear strain–displacement relations including the second order von Karman type nonlinearities are, therefore, defined by Eq. (2) 3

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 4. Cubic basis functions, 𝑁𝑖,3 (𝜉), of the B-spline curve of Fig. 3.

Fig. 5. Characterization of a curved beam element in IGA.

coefficient. For a system in equilibrium, it is necessary and sufficient for the total potential energy to be stationary. In other words, the first variation of the total potential energy should be equal to zero [65–67].

𝛿Π = 𝛿U − 𝛿W = 0

and a numerical approach is essential. For this purpose, the isogeometric analysis is presented in order to discretize and solve the governing equations. 2.3. A brief review of B-spline curves

(6)

As a result, substituting Eqs. (4) and (5) into Eq. (6) and utilizing the calculus of variations, the weak forms of the displacement fields of the curved beam (𝑢, 𝑤, 𝜑) can be described as follows, where 𝑄𝑠 , 𝑄𝑧 and 𝑄𝑚 are boundary loads. ) [ ( 1 1 𝐸𝐴𝛿𝑢′ 𝑢′ − 𝜅𝑤 + 𝑢′2 + 𝑤′2 2 2 ) ( 1 1 + 𝐸𝐴𝑢′ 𝛿𝑢′ 𝑢′ − 𝜅𝑤 + 𝑢′2 + 𝑤′2 2) 2 ] ( + 𝐾𝑠 𝐺𝐴𝜅𝛿𝑢 𝜅𝑢 + 𝑤′ − 𝜑 − 𝑓𝑠 𝛿𝑢 𝑑𝑠 − 𝛿𝑢𝑄𝑠 |𝑙0 = 0 ( ) 𝑙[ 1 1 −𝐸𝐴𝜅𝛿𝑤 𝑢′ − 𝜅𝑤 + 𝑢′2 + 𝑤′2 ∫0 2 2 ( ) 1 1 + 𝐸𝐴𝑤′ 𝛿𝑤′ 𝑢′ − 𝜅𝑤 + 𝑢′2 + 𝑤′2 2 2 ( ) ] + 𝐾𝑠 𝐺𝐴𝛿𝑤′ 𝜅𝑢 + 𝑤′ − 𝜑 − 𝑓𝑧 𝛿𝑤 𝑑𝑠 − 𝛿𝑤𝑄𝑧 |𝑙0 = 0

B-spline curves play a central role in the isogeometric analysis of free-form curved beams since they serve both for representing the geometry and for expressing the displacement fields (i.e. the shape functions of the field variables of interest, 𝑢, 𝑤 and 𝜑, are the same as the basis functions defining the geometry). B-splines are generally interpreted as unweighted or non-rational form of NURBS curves that can model a wide variety of free-form geometries. They are also well consistent with commercial CAD software and employed in different engineering problems (see, e.g., [68–72]). Starting with the curve definition, a B-spline curve of degree 𝑝, as a piecewise continuous function with 𝑛 + 1 control points 𝑷 𝑖 = [𝑋𝑖 , 𝑌𝑖 ] is expressed as follows [73]:

𝑙

∫0

𝑙

∫0

[ ( ) ] 𝐸𝐼𝜑′ 𝛿𝜑′ − 𝐾𝑠 𝐺𝐴𝛿𝜑 𝜅𝑢 + 𝑤′ − 𝜑 − 𝑚𝛿𝜑 𝑑𝑠 − 𝛿𝜑𝑄𝑚 |𝑙0 = 0

(7)

(8)

𝑪 (𝜉) =

𝑛 ∑

𝑁𝑖,𝑝 (𝜉)𝑷 𝑖

(10)

𝑖=0

where 𝑪 (𝜉) = [𝑥 (𝜉) , 𝑦 (𝜉)] is a vector-valued function whose components are represented separately as mappings of the curve parameter 𝜉 ∈ [0, 1] onto the geometric space, where the univariate parameter space is characterized by the knot vector Ξ with 𝑛 − 𝑝 + 1 non-zero knot spans. In addition, 𝑁𝑖,𝑝 (𝜉) as the 𝑖th basis function of degree 𝑝 is defined by the Cox–de Boor recursion formula. The knot vector and

(9)

Up to this stage, we have formulated the governing equations of equilibrium for a geometrically nonlinear Timoshenko curved beam. However, the aforementioned equations cannot be solved analytically 4

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 6. Pseudo-arclength continuation scheme in the vicinity of a limit point of an equilibrium path.

Fig. 7. A clamped–clamped shallow arch beam under central point load.

basis functions are expressed as [73]: Ξ = [0, 0, … , 0, 𝜉𝑝+1 , 𝜉𝑝+2 , … , 𝜉𝑛 , 1, 1, … , 1] ⏟⏞⏞⏟⏞⏞⏟ ⏟⏞⏞⏟⏞⏞⏟ 𝑝+1

{

𝑁𝑖,0 (𝜉) = 𝑁𝑖,𝑝 (𝜉) =

(11)

𝑝+1

1

𝜉𝑖 ≤ 𝜉 < 𝜉𝑖+1

0

otherwise

(12)

𝜉𝑖+𝑝+1 − 𝜉 𝜉 − 𝜉𝑖 𝑁 𝑁 (𝜉) + (𝜉) 𝜉𝑖+𝑝 − 𝜉𝑖 𝑖,𝑝−1 𝜉𝑖+𝑝+1 − 𝜉𝑖+1 𝑖+1,𝑝−1

As an instance, Fig. 3 depicts a cubic B-spline curve with eight control points and five non-zero knot spans, i.e. Ξ = [0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1]. For this curve, which can be thought of as the geometry of a curved beam in IGA, the first and last control points are coincident with the starting and ending points on the curve. The cubic basis functions, 𝑁𝑖,3 (𝜉), which are of order two continuous (𝐶 2 ) at knots, are also demonstrated in Fig. 4. Finally, it is interesting to note that if only a set of data points representing the geometry of a curve is available, the B-spline expression of the curved beam can be found by a curve fitting technique [73–75]. 2.4. Isogeometric discretization Referring to Fig. 5, within the 𝑘th element Ω𝑒 ∶ 𝑠 ∈ [𝑠𝑘−1 , 𝑠𝑘 ] on the curved beam that is equivalent to 𝜉 ∈ [𝜉𝑝+𝑘−1 , 𝜉𝑝+𝑘 ] on the parameter space, the local support property of B-spline curves necessitates having 𝑝 + 1 non-zero (i.e. active) basis functions, namely 𝑁𝑘−1,𝑝 (𝜉) to 𝑁𝑘+𝑝−1,𝑝 (𝜉) where the corresponding control points are 𝑷 𝑘−1 to 𝑷 𝑘+𝑝−1 [73]. Table 2 presents the geometry and displacement fields of curved beams in the isogeometric framework where 𝑥 and 𝑦 denote the position of an arbitrary point on the beam’s midline, [𝑋𝑖 , 𝑌𝑖 ] is the position of the 𝑖th control point of the geometry in the 𝑥–𝑦 plane, and 𝑈𝑖 , 𝑊𝑖 and Φ𝑖 are the control points of the field variables. Substituting displacement fields of Table 2 in the weak forms of Eqs. (7) to (9) and integrating over the element domain, the force– displacement relation at the element level will be obtained:

Fig. 8. Bifurcation analysis of the shallow arch example under a central point load and with different numbers of control points: (a) equilibrium path validated against Lo [10] and ABAQUS, (b) convergence study at limit points A and D.

𝑲 𝑒 𝒒𝑒 = 𝑭 𝑒 5

(13)

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 9. IGA results of clamped–clamped shallow arch under a central point load: (a) effect of different beam thicknesses, (b) undeformed vs. deformed shapes at various stages of the solution for 𝐻 = 3∕16 in.

Fig. 10. IGA results of clamped–hinged shallow arch under a central point load: (a) effect of different beam thicknesses, (b) undeformed vs. deformed configurations at various stages of the solution for 𝐻 = 3∕16 in.

6

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

In the above integrals, the integration is with respect to the curvilinear coordinate 𝑠, while the integrands are in terms of the curve parameter 𝜉. As a result, the transformation 𝑑𝑠 = 𝐽 𝑑𝜉 from the curvilinear coordinate on the geometric space into the curve parameter should be performed by the Jacobian (𝐽 ), defined as: √ ( )2 ( )2 𝜕𝑦 𝜕𝑥 𝑑𝑠 + (18) = 𝐽= 𝑑𝜉 𝜕𝜉 𝜕𝜉 Hence, each component of the stiffness matrix and load vector can be determined using the standard Gauss–Legendre quadrature rule as a numerical integration technique. For this purpose, as depicted in Fig. 5, another transformation from the curve parameter space 𝜉 into the parent element 𝜉̃ is needed as follows: ) )] ( 1 [( (19) 𝜉= 𝜉𝑖+1 − 𝜉𝑖 𝜉̃ + 𝜉𝑖+1 + 𝜉𝑖 2 Given the definition of the Jacobian, one can obtain: ( ) 1 𝜕𝑥 𝜕 2 𝑦 𝜕𝑦 𝜕 2 𝑥 𝜅= − (20) 3 2 2 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝐽 𝑑𝑁𝑖 𝑑𝑁𝑖 𝑑𝜉 1 𝑑𝑁𝑖 𝑁𝑖′ = = = (21) 𝑑𝑠 𝑑𝜉 𝑑𝑠 𝐽 𝑑𝜉

Fig. 11. Effect of different beam thicknesses on the IGA results of the clamped–clamped shallow arch under a distributed load.

where 𝒒 𝑒 is the vector of control points of the field variables at the element level: [ ]𝑇 ⎧ 𝑼 ⎫ ⎧ 𝑈 ,…,𝑈 ⎫ 𝑘−1 𝑘+𝑝−1 ⎪ ⎪ ⎪ [ ]𝑇 ⎪ 𝑒 (14) 𝒒 = ⎨ 𝑾 ⎬ = ⎨ 𝑊𝑘−1 , … , 𝑊𝑘+𝑝−1 ⎬ ]𝑇 ⎪ ⎪ Φ ⎪ ⎪ [Φ , … , Φ ⎩ ⎭ ⎩ ⎭ 𝑘−1 𝑘+𝑝−1

where ∑ 𝑑 𝑟 𝑁𝑖 𝜕𝑟 𝑥 = 𝑋 𝑟 𝜕𝜉 𝑑𝜉 𝑟 𝑖 𝑖=0 𝑛

∑ 𝑑 𝑟 𝑁𝑖 𝜕𝑟 𝑦 = 𝑌 𝜕𝜉 𝑟 𝑑𝜉 𝑟 𝑖 𝑖=0 𝑛

and ⎡𝑲 11 𝑲 𝑒 = ⎢𝑲 21 ⎢ 31 ⎣𝑲

𝑲 12 𝑲 22 𝑲 32

𝑲 13 ⎤ 𝑲 23 ⎥ ⎥ 𝑲 33 ⎦

and

[ ] 𝑁𝑖,𝑝−1 (𝜉) 𝑁𝑖+1,𝑝−1 (𝜉) 𝑑𝑁𝑖 =𝑝 − 𝑑𝜉 𝜉𝑖+𝑝 − 𝜉𝑖 𝜉𝑖+𝑝+1 − 𝜉𝑖+1

(15)

⎧𝑭 1 ⎫ ⎪ ⎪ 𝑭 𝑒 = ⎨𝑭 2 ⎬ ⎪𝑭 3 ⎪ ⎩ ⎭

𝐾𝑖𝑗23

=

𝐹𝑖1 = 𝐹𝑖2 = 𝐹𝑖3 =

𝐾𝑗𝑖32

∫Ωe ∫Ωe ∫Ωe

=

∫Ωe

(23)

Having included geometric nonlinearities in the analysis of curved beams, the components of the stiffness matrix of Eq. (16) are functions of the unknown displacement field variables. In other words, the assembled force–displacement relationship becomes 𝑲 (𝒒) 𝒒 = 𝑭 . Therefore, an iterative/incremental method for solving nonlinear system of equations should be employed to find the unknown displacements. However, the conventional Newton–Raphson approach can only be helpful in the pre-buckling region of the equilibrium path of the structure as it breaks down on the bifurcation points. For overcoming this difficulty, the arclength continuation techniques should be employed. In this paper, we suggest using the pseudo-arclength formulation as described in the following section.

The explicit forms of the stiffness matrix and load vector at the element level are as follows, noting that for brevity, we write 𝑁𝑖 for 𝑁𝑖,𝑝 (𝜉) and 𝑑𝑁 𝑁𝑖′ for 𝑑𝑠𝑖 where 𝑖, 𝑗 = 𝑘 − 1, … , 𝑘 + 𝑝 − 1. ( 3 𝐾𝑖𝑗11 = 𝐸𝐴𝑁𝑖′ 𝑁𝑗′ + 𝐸𝐴𝑢′ 𝑁𝑖′ 𝑁𝑗′ ∫Ωe 2 ) 1 + 𝐸𝐴𝑢′2 𝑁𝑖′ 𝑁𝑗′ + 𝐾𝑠 𝐺𝐴𝜅 2 𝑁𝑖 𝑁𝑗 𝑑𝑠 2 ( 1 (16) 𝐾𝑖𝑗22 = 𝐸𝐴𝜅 2 𝑁𝑖 𝑁𝑗 − 𝐸𝐴𝜅𝑤′ 𝑁𝑖 𝑁𝑗′ − 𝐸𝐴𝜅𝑤′ 𝑁𝑖′ 𝑁𝑗 ∫Ωe 2 ) 1 + 𝐸𝐴𝑤′2 𝑁𝑖′ 𝑁𝑗′ + 𝐾𝑠 𝐺𝐴𝑁𝑖′ 𝑁𝑗′ 𝑑𝑠 2 ( 1 12 𝐾𝑖𝑗 = −𝐸𝐴𝜅𝑁𝑖′ 𝑁𝑗 + 𝐸𝐴𝑤′ 𝑁𝑖′ 𝑁𝑗′ − 𝐸𝐴𝜅𝑢′ 𝑁𝑖′ 𝑁𝑗 ∫Ωe 2 ) 1 + 𝐸𝐴𝑢′ 𝑤′ 𝑁𝑖′ 𝑁𝑗′ + 𝐾𝑠 𝐺𝐴𝜅𝑁𝑖 𝑁𝑗′ 𝑑𝑠 2 ( 1 𝐾𝑖𝑗21 = −𝐸𝐴𝜅𝑁𝑖 𝑁𝑗′ − 𝐸𝐴𝜅𝑢′ 𝑁𝑖 𝑁𝑗′ + 𝐸𝐴𝑤′ 𝑁𝑖′ 𝑁𝑗′ ∫Ωe 2 ) 1 + 𝐸𝐴𝑤′ 𝑢′ 𝑁𝑖′ 𝑁𝑗′ + 𝐾𝑠 𝐺𝐴𝜅𝑁𝑖′ 𝑁𝑗 𝑑𝑠 2 ( ) 𝐾𝑖𝑗33 = 𝐸𝐼𝑁𝑖′ 𝑁𝑗′ + 𝐾𝑠 𝐺𝐴𝑁𝑖 𝑁𝑗 𝑑𝑠 ∫Ωe 𝐾𝑖𝑗13 = 𝐾𝑗𝑖31 =

(22)

3. IGA Integrated with pseudo-arclength continuation The assembled form of the force–displacement relationship is a system of 𝓃 nonlinear equations in terms of 𝓃 unknowns through the [ ]𝑇 vector 𝒒 = 𝑞1 , 𝑞2 , … , 𝑞𝓃 where 𝓃 is the total number of degrees of freedom. 𝒇 = 𝑲𝒒 − 𝑭 = 0

(24)

By defining the force vector as 𝑭 = 𝜆𝑭 where 𝑭 is the fixed external load vector and 𝜆 is the load multiplier, the above-mentioned equilibrium equation can be written as follows noting that the set (𝒒, 𝜆) describes one of the equilibrium states on the equilibrium path.

−𝐾𝑠 𝐺𝐴𝜅𝑁𝑖 𝑁𝑗 𝑑𝑠 𝒇 (𝒒, 𝜆) = 0 (25) ( ∗ ∗) If 𝒒 , 𝜆 is a known equilibrium state, according to the implicit func( ) tion theorem, there exists a neighborhood of 𝒒 ∗ , 𝜆∗ on the equilibrium path such as (𝒒 (𝜎) , 𝜆(𝜎)) that can be tracked by the continuation of the solution where 𝜎 is the arclength parameter along the path (Fig. 6). Utilizing the pseudo-arclength continuation scheme [5,6], the next equilibrium state (𝒒 (𝜎) , 𝜆(𝜎)) of the buckled curved beam is sought as the intersection of the line normal to the tangent vector passing through

−𝐾𝑠 𝐺𝐴𝑁𝑖′ 𝑁𝑗 𝑑𝑠 ∫Ωe

( ) 𝑁𝑖 𝑓𝑠 𝑑𝑠 + 𝑄𝑠 𝑙𝑒 − 𝑄𝑠 (0) ( ) 𝑁𝑖 𝑓𝑧 𝑑𝑠 + 𝑄𝑧 𝑙𝑒 − 𝑄𝑧 (0)

(17)

( ) 𝑁𝑖 𝑚𝑑𝑠 + 𝑄𝑚 𝑙𝑒 − 𝑄𝑚 (0) 7

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 12. Bifurcation behavior of the semi-circular arch: (a) load vs. vertical displacement validated against Addessi et al. [11] and ABAQUS, (b) load vs. horizontal displacement, (c) undeformed vs. deformed configurations of the buckled beam.

8

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 13. Effect of various imperfection values on the bifurcation points of the equilibrium path of semi-circular arch.

Fig. 14. A free-form Tschirnhausen beam under a distributed normal load.

( ∗ ) 𝒒 + 𝒒 𝜎 Δ𝜎, 𝜆∗ + 𝜆𝜎 Δ𝜎 and the equilibrium path (see Fig. 6) where Δ𝜎 is the step length of the incremental solution and 𝒒 𝜎 and 𝜆𝜎 are the derivatives of 𝒒 and 𝜆 with respect to the arclength parameter. Consequently, the new equilibrium state can be found by solving the following system of 𝓃 + 1 nonlinear equations in terms of 𝓃 + 1 unknowns (𝒒, 𝜆) in a predicting–correcting manner. { 𝒇 (𝒒, 𝜆) = 0 (26) ( )𝑇 ( ) 𝒈 (𝒒, 𝜆) = 𝒒 − 𝒒 ∗ 𝒒 𝜎 + 𝜆 − 𝜆∗ 𝜆𝜎 − Δ𝜎 = 0 ( ) In the above equations, 𝒒 𝜎 and 𝜆𝜎 should be calculated from 𝒒 ∗ , 𝜆∗ as follows: ( )−1∕2 𝜆𝜎 = ± 1 + 𝒛𝑇 𝒛 (27)

Hence, by applying the above equation to the stiffness components of Eq. (16), one can write: ( 3 𝑇𝑖𝑗11 = 𝐾𝑖𝑗11 + 𝐸𝐴𝑢′2 𝑁𝑖′ 𝑁𝑗′ + 𝐸𝐴𝑢′ 𝑁𝑖′ 𝑁𝑗′ − 𝐸𝐴𝜅𝑤𝑁𝑖′ 𝑁𝑗 ∫Ωe 2 ) 1 + 𝐸𝐴𝑤′2 𝑁𝑖′ 𝑁𝑗′ 𝑑𝑠 2 ( 1 1 𝑇𝑖𝑗22 = 𝐾𝑖𝑗22 + 𝐸𝐴𝑢′ 𝑁𝑖′ 𝑁𝑗′ + 𝐸𝐴𝑢′2 𝑁𝑖′ 𝑁𝑗′ − 𝐸𝐴𝜅𝑤′ 𝑁𝑖 𝑁𝑗′ ∫Ωe 2 2 ) ′ ′ ′2 ′ ′ − 𝐸𝐴𝜅𝑤𝑁𝑖 𝑁𝑗 + 𝐸𝐴𝑤 𝑁𝑖 𝑁𝑗 𝑑𝑠 ) ( 1 1 𝐸𝐴𝑤′ 𝑁𝑖′ 𝑁𝑗′ + 𝐸𝐴𝑢′ 𝑤′ 𝑁𝑖′ 𝑁𝑗′ 𝑑𝑠 𝑇𝑖𝑗12 = 𝐾𝑖𝑗12 + (31) ∫Ωe 2 2 ) ( 1 1 𝑇𝑖𝑗21 = 𝐾𝑖𝑗21 + − 𝐸𝐴𝜅𝑢′ 𝑁𝑖 𝑁𝑗′ + 𝐸𝐴𝑢′ 𝑤′ 𝑁𝑖′ 𝑁𝑗′ 𝑑𝑠 ∫Ωe 2 2

(28)

𝒒 𝜎 = 𝜆𝜎 𝒛

𝑇𝑖𝑗33 = 𝐾𝑖𝑗33

where the plus and minus signs determine the direction of the continuation and 𝒛 is the solution of the following linear algebraic equations noting that 𝒇 𝒒 is the so-called tangent stiffness matrix and 𝒇 𝜆 is the derivative vector of the equilibrium equation with respect to the control parameter 𝜆. ( ) 𝒇 𝒒 𝒒 ∗ , 𝜆∗ 𝒛 = −𝒇 𝜆 (𝒒 ∗ , 𝜆∗ ) (29)

𝑇𝑖𝑗13 = 𝑇𝑗𝑖31 = 𝐾𝑖𝑗13 = 𝐾𝑗𝑖31 𝑇𝑖𝑗23 = 𝑇𝑗𝑖32 = 𝐾𝑖𝑗23 = 𝐾𝑗𝑖32 Returning to the nonlinear system of equations (26), a Newton–Raphson scheme is needed to be applied in order to find the equilibrium state (𝒒, 𝜆) in an iterative manner as follows where the superscript 𝑡 is the iteration number and 0 < 𝑟 ≤ 1 is the relaxation parameter chosen such )‖ ‖ ( )‖ ‖ ( that ‖𝒇 𝒒 𝑡+1 , 𝜆𝑡+1 ‖ < ‖𝒇 𝒒 𝑡 , 𝜆𝑡 ‖. ‖ ‖ ‖ ‖

Considering the isogeometric formulation of Eqs. (14)–(16), the tangent stiffness matrix 𝑻 𝑒 = 𝒇 𝒒𝑒 at the 𝑘th element is to be calculated from the components of 𝑲 𝑒 as follows [1] where 𝛼, 𝛽 = 1, 2, 3 and 𝑞𝑗1 = 𝑈𝑗 , 𝑞𝑗2 = 𝑊𝑗 , 𝑞𝑗3 = Φ𝑗 . 𝑇𝑖𝑗𝛼𝛽 = 𝐾𝑖𝑗𝛼𝛽 +

3 𝑘+𝑝−1 𝛼𝑚 ∑ ∑ 𝜕𝐾𝑖ℎ 𝑚=1 ℎ=𝑘−1

𝜕𝑞𝑗𝛽

𝑞ℎ𝑚

(30)

9

𝒒 𝑡+1 = 𝒒 𝑡 + 𝑟Δ𝒒 𝑡+1

(32)

𝜆𝑡+1 = 𝜆𝑡 + 𝑟Δ𝜆𝑡+1

(33)

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 15. Mechanical behavior of the Tschirnhausen beam under uniform loading: (a) and (b) force–displacement graphs validated against ABAQUS, (c) depiction of undeformed shapes at various stages of the solution.

Fig. 16. Effect of different slenderness ratios on (a) vertical and (b) horizontal displacement of midspan of Tschirnhausen beam.

Once Δ𝜆𝑡+1 and Δ𝒒 𝑡+1 are known, the iterations are continued until the required convergence is achieved in the Newton–Raphson procedure. After the set (𝒒, 𝜆) is computed, we can move on to determine another equilibrium state on the equilibrium path. More mathematical details on the pseudo-arclength continuation can be found in [5,6].

In the above equations, Δ𝜆𝑡+1 and Δ𝒒 𝑡+1 can be determined from: ( ) 𝒈 𝒒 𝑘 , 𝜆𝑘 + 𝒛𝑇1 𝒒 𝜎 Δ𝜆𝑡+1 = − (34) 𝜆𝜎 + 𝒛𝑇2 𝒒 𝜎 Δ𝒒 𝑡+1 = 𝒛1 + 𝒛2 Δ𝜆𝑡+1

(35)

where 𝒛1 and 𝒛2 are the solutions of the following linear systems of equations: ( ) ( ) 𝒇 𝒒 𝒒 𝑡 , 𝜆𝑡 𝒛1 = −𝒇 𝒒 𝑡 , 𝜆𝑡 (36) ( 𝑡 𝑡) ( 𝑡 𝑡) 𝒇 𝒒 𝒒 , 𝜆 𝒛2 = −𝒇 𝜆 𝒒 , 𝜆 (37)

4. Case studies The presented methodology of this article will be investigated by different case studies in this section. In the first two benchmarked examples, namely circular shallow and high arches, the IGA results 10

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 19. An L-shaped beam under uniform distributed loading.

Fig. 17. Hinged–hinged logarithmic spiral beam under a clockwise bending moment at the top end.

are validated against FEA results published by Lo [10] and Addessi et al. [11], respectively. In order to illustrate the applicability of the presented framework, the next two examples are selected from the wellknown free-form geometries, namely Tschirnhausen and logarithmic spiral beams. The bifurcation behavior of an L-shaped beam as a geometry with a sharp corner is also investigated as the last example. The results in all case studies are compared with FEA computations of the commercial ABAQUS software, obtained by generating appropriate fine meshes of 500 quadratic elements. It is also important to mention that all examples are modeled by cubic curves, considering that geometries of the third or higher degrees generally alleviate the locking issues of the IGA approach [76–78]. The numbers of control points, which govern the number of DOFs and may be variable for different examples, are selected in such a manner as to reach a suitable convergence in results.

Fig. 20. Effect of decreasing 𝜖 on modeling an L-shaped beam with a sharp corner.

4.1. Circular shallow arch

The vertical displacement at the apex (point of load), as the output of the analysis, is computed using the continuation technique integrated by the IGA approach with different number of control points and plotted versus the applied load in Fig. 8a. As it is clear in the figure, the shallow arch represents a snap-through behavior due to saddle–node bifurcation. In other words, after limit point A is reached, there is no

The large deflection of a circular shallow arch is analyzed as the first example. In order to validate the presented formulation, a clamped– clamped arch under a concentrated load at the middle is modeled with dimensional and material data provided in [10] (see Fig. 7).

Fig. 18. Mechanical behavior of the logarithmic spiral beam under a clockwise bending moment at the top end: (a) nonlinear moment–rotation diagram validated against ABAQUS, (b) depiction of undeformed shapes of the beam.

11

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 21. Mechanical behavior of the L-shaped beam for 𝑎 = 0.15 L : normalized vertical (a) and horizontal (b) displacements of point of load vs. normalized load, (c) deformed shapes of the structure at some arbitrary equilibrium states.

Fig. 10a. As demonstrated in the figure, having asymmetric boundary conditions, a snap-back behavior can also be seen in the equilibrium path. More precisely, the beam represents a snap-through instability at limit point A under force control, while if the loading is performed under a controlled displacement, it reveals a snap-back instability at turning point B. As also described for the symmetric case, by increasing the beam thickness, the structure becomes stiffer and the saddle–node bifurcation behavior is gradually alleviated. The deformed configurations of the beam at various stages of the solution, for the case 𝐻 = 3∕16 in, are also demonstrated in Fig. 10b. The instability of shallow arches can also be studied under distributed loading. Fig. 11 depicts the equilibrium paths of the current example with clamped–clamped boundary conditions and different beam thicknesses subject to a uniform distributed loading 𝜔0 . The results again show that increasing the thickness would result in a stable response in the beam. In this figure, the vertical axis is also scaled by the ratio of moments of inertia (𝐼1 ∕𝐼).

adjacent equilibrium state and the system snaps through an additional stable post-buckling configuration. There is a similar interpretation for limit point D during the unloading procedure noting that the interval A to D with negative slope is an unstable equilibrium state and cannot be reached under force control. The FEA results of Lo [10] and ABAQUS software are also superimposed onto the equilibrium path. In addition, a convergence study is performed in order to find the appropriate number of control points by comparing the computation errors of the proposed methodology at limit points A and D with respect to reference FEA solution of ABAQUS. This comparison shows that by modeling the arch with 50 control points, a suitable convergence in IGA results is obtained (see Fig. 8b). The results of this example prove the validity of the developed methodology in static bifurcation analysis of curved beams with geometric nonlinearity. In order to assess another involving parameter changing the bifurcation behavior of the problem, we investigate the effect of different beam thicknesses. For this purpose, Fig. 9a depicts the load–displacement results of the current example with different thicknesses where the vertical axes is scaled by the basic moment of inertia 𝐼1 of the case 𝐻 = 3∕16 in. It is inferred from the figure that a shallow arch with lower thickness (i.e. lower stiffness) is more prone to the snap-through instability, while if the beam is thick enough, it can be totally stable. The deformed shapes of the beam at various stages of the solution are also demonstrated in Fig. 9b for the case 𝐻 = 3∕16 in. For showing the effect of different boundary conditions, the clamped–hinged arch is subject to a central point load and the force– displacement diagrams for different thickness values are compared in

4.2. Semi-circular high arch The mechanical behavior of non-shallow or so-called high arches is distinctly different from that of shallow arches. Despite shallow arches that are often prone to snap-through instabilities, high arches or non-shallow curved beams in general, may exhibit lateral buckling because they are very sensitive to imperfections. In the presence of a slight imperfection, they usually show a pitchfork bifurcation behavior at a branch point before reaching the expected limit points on their 12

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

Fig. 22. Effect of changing the position of applied load on the (a) vertical and (b) horizontal displacements of point of load in the L-shaped beam example.

fundamental path. Generally, the rank of the 𝓃 × (𝓃 + 1) augmented [ ] matrix 𝒇 𝒒 𝒇 𝜆 is used to determine if the bifurcation point is a limit or a branch point, where 𝓃 is the number of DOFs. This matrix has a rank of 𝓃 at saddle–node and 𝓃 − 1at pitchfork bifurcation points [6]. In the second case study of this paper, a semi-circular high arch of radius 𝑅 with hinged–hinged boundary conditions subject to a point load 𝑃 at the apex is considered. As demonstrated in Fig. 12a, a pitchfork bifurcation is detected at point B and the equilibrium state deviates from its fundamental path as a result of a small imperfection. The imperfection is added to the model by means of a slight horizontal load 𝑄 = 0.001𝑃 . The equilibrium path generated by the developed methodology shows a great consistency with FEA results of Addessi et al. [11] and ABAQUS software. Passing the branch point B, two stable buckled states exist (see Fig. 12b). Although the boundary conditions are symmetric, the beam has an asymmetric deformed shape and buckles to either right or left. By increasing the applied load, the buckled stable state becomes unstable at the limit point D that

corresponds to a saddle–node bifurcation. Showing such response is due to the fact that the beam is sufficiently shallow at this stage and behaves like a shallow arch. The deformed configurations of the buckled beam are also depicted in Fig. 12c. It is worth noticing that in the real-world structures, an imperfection can be interpreted as a slight deviation in the position or direction of the central load, or in boundary conditions. As illustrated in Fig. 13, imperfections up to 0.01𝑃 has a slight effect on the branch point, but with no effect on the limit point. However, increasing the magnitude of 𝑄 to greater values will significantly change the system response. 4.3. Tschirnhausen free-form beam The Tschirnhausen beam as a free-form geometry with variable curvature is the third example of this article. In order to predict the bifurcation behavior, the beam is considered to be under a constant distributed load per unit length 𝜔0 that is always normal to the curve 13

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16

take the form of Eq. (40) where 𝜖 is a small quantity [44]. As depicted in Fig. 20, by decreasing 𝜖 toward zero the geometry approaches its exact predefined shape. Hence, a convergence check is required to obtain the appropriate value of 𝜖, which is selected as 0.001 in this example.

(see Fig. 14). The Tschirnhausen geometry is generated using the following parametric relations that guarantee a fixed rise to span ratio for the beam. It is interesting to note that in order to apply loads in the Cartesian 𝑥–𝑦 configuration, we need to employ a rotation matrix where the rotation angle is calculated by obtaining the unit tangent vector (i.e., the norm of the first derivative of geometry). Conversely, to transform the displacement results from the curvilinear coordinate system to the Cartesian frame, the inverse (i.e., transpose) of the rotation matrix should be employed. ( ) √ 𝑥 = 3𝑎 3 − 𝑡2 ( ) (38) 0≤𝑡≤ 3 𝑦 = 𝑎𝑡 3 − 𝑡2

Ξ = [0, … , 0.5 − 𝜖, 0.5, 0.5 + 𝜖, … , 1]

In order to analyze the load–displacement behavior of this example, a downward point load is exerted at a distance 𝑎 from the corner. When the point of load is near the corner, because of the shape of the structure, the buckling phenomenon is predictable. In Fig. 21, the normalized vertical and horizontal displacements (𝑤∕𝐿 and 𝑢∕𝐿, respectively) of point of load are plotted versus the normalized load exerted at distance 𝑎 = 0.15𝐿 on the beam. The results are also validated against those computed by the FEA software. The deformed shapes of the structure at some arbitrary equilibrium states are also demonstrated in the figure. It is also interesting to note that by changing the position of the applied load, the beam shows a different geometric nonlinearity trend. As demonstrated in Fig. 22, for larger values of 𝑎, the deformed beam becomes gradually stiffer by increasing the applied load and exhibits a hardening-type nonlinearity. As a result, no saddle–node bifurcation will be observed for larger 𝑎∕𝐿 ratios. The results of this study show the superiority of the representation of IGA in solving nonlinear problems of beams with sharp corners.

The beam in this example is assumed to have a circular cross section and hinged boundary conditions at both ends. The mechanical behavior of the Tschirnhausen beam is predicted by the presented IGA approach and demonstrated in Fig. 15. In parts (a) and (b), the normalized displacements of the midspan are considered as the output and plotted versus the normalized distributed force. The results, which are also validated against FEA computations in ABAQUS, show that a snapthrough instability is inevitable for the structure. The deformed shapes of the beam at various stages of the solution are also depicted in Fig. 15c. It is inferred from the equilibrium path that the beam deflects like a shallow curved beam and as a result, we expect a dependency in response to beam slenderness. Fig. 16 investigates the effect of different slenderness ratios on the bifurcation behavior of the Tschirnhausen beam. The results show that the instability will be lightened for thicker and shorter beams. In this comparison, we employed the span √ 𝐿 instead of the beam length in calculation of the slenderness ratio (𝛼 = 𝐴𝐿2 ∕𝐼) .

5. Conclusions This paper aims at presenting an efficient isogeometric framework for static bifurcation analysis of geometrically nonlinear curved beams. Considering the importance of free-form geometries in the modern engineering applications, the main contribution of this research lies in extending the current literature of the nonlinear buckling analysis of curved beams to a larger range of practical free-form geometries with variable curvature. In this regard, the isogeometric analysis (IGA) is integrated with pseudo-arclength continuation technique in order to predict the full equilibrium path of the structures including two common types of static instabilities, namely the saddle–node and pitchfork bifurcations. The developed formulation is based on the Timoshenko beam theory and accounts for shear deformations in the large-displacement analysis of curved beams. The advantage of the presented method over conventional finite element solution is that IGA helps reduce the number of elements and therefore unknown DOFs in the analysis, enhancing the convergence rate and reducing the computational effort of the continuation technique. In addition, the isoparametric formulation of the IGA framework (i.e., using the same basis functions for both geometry and solution), increases the accuracy of the nonlinear bifurcation analysis, particularly, in the case of free-form curved beams. The validity of the presented formulations is successfully tested in the well-established shallow and non-shallow circular arch examples. Additionally, the Tschirnhausen and logarithmic spiral beams as freeform curved geometries and an L-shaped beam as a geometry with a sharp corner are also investigated under various loading conditions noting that in all examples, good agreements with FEA results of the commercial ABAQUS software are achieved. Finally, as future works, it is suggested to extend the current bifurcation analysis of this article to spatial beams with large deformations or to bivariate problems like free-form shells.

4.4. Logarithmic spiral free-form beam Spiral curved shapes in various forms are frequently used in engineering applications. Therefore, the logarithmic spiral beam as a free-form geometry with variable curvature is selected as the fourth case study of this article. The beam geometry is parametrically defined by Eq. (39) in the 𝑥–𝑦 coordinate system. A curve fitting algorithm should then be employed to generate the appropriate geometry for IGA. As depicted in Fig. 17, the hinged–hinged logarithmic spiral beam is made of a steel tube and supposed to be under a clockwise bending moment at the top end. 𝑥 = 2.41 cos 𝑡 𝑒0.405𝑡 0 ≤ 𝑡 ≤ 2.5𝜋 𝑦 = 2.41 sin 𝑡 𝑒0.405𝑡

(40)

(39)

The capability of the current work in dealing with different loading scenarios is investigated in this example. An important point to consider is that, in the presented IGA framework, the Timoshenko beam element has a rotational DOF as illustrated earlier in Fig. 2, so that we can input the bending moment directly in the load vector. The results of the static bifurcation analysis for the logarithmic spiral beam are shown in Fig. 18. The static equilibrium path, which is in a good agreement with FEA results (see Fig. 18a), shows a few numbers of turning points as a result of saddle–node instability. It is worth mentioning that the continuation beyond point G is computationally possible but is not physically meaningful since it would result in self-intersection (see Fig. 18b).

References

4.5. L-shaped beam with a sharp corner

[1] J.N. Reddy, An Introduction to Nonlinear Finite Element Analysis, Oxford University Press, Oxford, UK, 2004. [2] M.A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, John Wiley & Sons, New York, NY, 1991. [3] R. Levy, W.R. Spillers, Analysis of Geometrically Nonlinear Structures, second ed., Springer Dordrecht, Netherlands, 2003. [4] Z.P. Bažant, L. Cedolin, Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories, World Scientific, Singapore, 2010.

In this example, an L-shaped beam undergoing a point load (see Fig. 19) is analyzed as a geometry with a sharp corner. The geometry of the problem is constructed by combining the B-spline definitions of two equal straight lines in horizontal and vertical directions. In order for the sharp corner to be modeled as close to reality as possible and also to have a smooth field transition between two lines, the knot vector should 14

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16 [36] D. Wang, Q. Liang, H. Zhang, A superconvergent isogeometric formulation for eigenvalue computation of three dimensional wave equation, Comput. Mech. 57 (2016) 1037–1060. [37] R. Bouclier, T. Elguedj, A. Combescure, Locking free isogeometric formulations of curved thick beams, Comput. Methods Appl. Mech. Engrg. (2012) 245–246, 144–162. [38] G. Zhang, R. Alberdi, K. Khandelwal, On the locking free isogeometric formulations for 3-D curved Timoshenko beams, Finite Elem. Anal. Des. 143 (2018) 46–65. [39] S.F. Hosseini, B. Moetakef-Imani, S. Hadidi-moud, B. Hassani, Pre-bent shape design of full free-form curved beams using isogeometric method and semi-analytical sensitivity analysis, Struct. Multidiscip. Optim. 58 (2018) 2621–2633. [40] A.P. Nagy, M.M. Abdalla, Z. Gürdal, Isogeometric sizing and shape optimisation of beam structures, Comput. Methods Appl. Mech. Engrg. 199 (2010) 1216–1230. [41] A.M. Bauer, M. Breitenberger, B. Philipp, R. Wüchner, K.U. Bletzinger, Nonlinear isogeometric spatial Bernoulli beam, Comput. Methods Appl. Mech. Engrg. 303 (2016) 101–127. [42] S.F. Hosseini, A. Hashemian, B. Moetakef-Imani, S. Hadidimoud, Isogeometric analysis of free-form Timoshenko curved beams including the nonlinear effects of large deformations, Acta Mech. Sinica 34 (2018) 728–743. [43] A.-T. Luu, N.-I. Kim, J. Lee, Isogeometric vibration analysis of free-form Timoshenko curved beams, Meccanica 50 (2015) 169–187. [44] S.F. Hosseini, B. Moetakef-Imani, S. Hadidi-Moud, B. Hassani, The effect of parameterization on isogeometric analysis of free-form curved beams, Acta Mech. 227 (2016) 1983–1998. [45] X. Zhang, Y. Xia, Q. Hu, P. Hu, Efficient isogeometric formulation for vibration analysis of complex spatial beam structures, Eur. J. Mech. A Solids 66 (2017) 212–231. [46] G. Balduzzi, S. Morganti, F. Auricchio, A. Reali, Non-prismatic Timoshenko-like beam model: Numerical solution via isogeometric collocation, Comput. Math. Appl. 74 (2017) 1531–1541. [47] S.F. Hosseini, A. Hashemian, A. Reali, On the application of curve reparameterization in isogeometric vibration analysis of free-from curved beams, Comput. Struct. 209 (2018) 117–129. [48] S. Shojaee, N. Valizadeh, E. Izadpanah, T. Bui, T.-V. Vu, Free vibration and buckling analysis of laminated composite plates using the NURBS-based isogeometric finite element method, Compos. Struct. 94 (2012) 1677–1693. [49] T. Le-Manh, J. Lee, Postbuckling of laminated composite plates using NURBS-based isogeometric analysis, Compos. Struct. 109 (2014) 286–293. [50] T.T. Yu, S. Yin, T.Q. Bui, S. Hirose, A simple FSDT-based isogeometric analysis for geometrically nonlinear analysis of functionally graded plates, Finite Elem. Anal. Des. 96 (2015) 1–10. [51] T. Yu, S. Yin, T.Q. Bui, C. Liu, N. Wattanasakulpong, Buckling isogeometric analysis of functionally graded plates under combined thermal and mechanical loads, Compos. Struct. 162 (2017) 54–69. [52] P. Ribeiro, Hierarchical finite element analyses of geometrically non-linear vibration of beams and plane frames, J. Sound Vib. 246 (2001) 225–244. [53] B. Krauskopf, H.M. Osinga, J. Galan-Vioque, Numerical Continuation Methods for Dynamical Systems, Springer and Canopus Publishing Limited, New York, NY, 2007. [54] R.M.J. Groh, D. Avitabile, A. Pirrera, Generalised path-following for well-behaved nonlinear structures, Comput. Methods Appl. Mech. Engrg. 331 (2018) 394–426. [55] B.S. Cox, R.M.J. Groh, D. Avitabile, A. Pirrera, Exploring the design space of nonlinear shallow arches with generalised path-following, Finite Elem. Anal. Des. 143 (2018) 1–10. [56] E. Riks, The Application of Newton’s method to the problem of elastic stability, J. Appl. Mech. 39 (1972) 1060–1065. [57] E. Riks, An incremental approach to the solution of snapping and buckling problems, Int. J. Solids Struct. 15 (1979) 529–551. [58] G.A. Wempner, Discrete approximations related to nonlinear theories of solids, Int. J. Solids Struct. 7 (1971) 1581–1599. [59] M.A. Crisfield, A fast incremental/iterative solution procedure that handles "snap-through", Comput. Struct. 13 (1981) 55–62. [60] Y.T. Feng, D. Perić, D.R.J. Owen, Determination of travel directions in path-following methods, Math. Comput. Modelling 21 (1995) 43–59. [61] Y.T. Feng, D. Perić, D.R.J. Owen, A new criterion for determination of initial loading parameter in arc-length methods, Comput. Struct. 58 (1996) 479–485. [62] P. Raveendranath, G. Singh, G. Venkateswara Rao, A three-noded shear-flexible curved beam element based on coupled displacement field interpolations, Internat. J. Numer. Methods Engrg. 51 (2001) 85–101. [63] K. Pan, J. Liu, Geometric nonlinear formulation for curved beams with varying curvature, J. Theor. Appl. Mech. Lett. 2 (2012) 063006. [64] K. Pan, J. Liu, Geometric nonlinear dynamic analysis of curved beams using curved beam element, Acta Mech. Sinica 27 (2011) 1023–1033. [65] J.M.T. Thompson, A general theory for the equilibrium and stability of discrete conservative systems, Z. Angew. Math. Phys. 20 (1969) 797–846. [66] C.H. Yoo, S. Lee, Stability of Structures: Principles and Applications, Butterworth-Heinemann, Boston, MA, 2011. [67] J.N. Reddy, Energy Principles and Variational Methods in Applied Mechanics, third ed., John Wiley & Sons, Inc., Hoboken, NJ, 2017.

[5] W. Lacarbonara, Nonlinear Structural Mechanics, Springer, New York, NY, 2013. [6] A.H. Nayfeh, B. Balachandran, Applied Nonlinear Dynamics, Wiley-Interscience, New York, NY, 1995. [7] I.G. Tadjbakhsh, Stability and optimum design of arch-type structures, Int. J. Solids Struct. 17 (1981) 565–574. [8] M.M. Attard, J. Zhu, D.C. Kellermann, In-plane buckling of prismatic funicular arches with shear deformations, Arch. Appl. Mech. 84 (2014) 693–713. [9] A.-T. Luu, N.-I. Kim, J. Lee, Bending and buckling of general laminated curved beams using NURBS-based isogeometric analysis, Eur. J. Mech. A Solids 54 (2015) 218–231. [10] S.H. Lo, Geometrically nonlinear formulation of 3D finite strain beam element with large rotations, Comput. Struct. 44 (1992) 147–157. [11] D. Addessi, W. Lacarbonara, A. Paolone, Linear vibrations of planar pre-stressed arches undergoing static bifurcations, in: Proc. Proceedings of the 6th European Conference on Structural Dynamics, EURODYN 2005, Paris, France, 2005. [12] Y.L. Pi, M.A. Bradford, Non-linear buckling and postbuckling analysis of arches with unequal rotational end restraints under a central concentrated load, Int. J. Solids Struct. 49 (2012) 3762–3773. [13] Y.L. Pi, M.A. Bradford, Nonlinear analysis and buckling of shallow arches with unequal rotational end restraints, Eng. Struct. 46 (2013) 615–630. [14] M.-Y. Kim, S.-B. Kim, N.-I. Kim, Spatial stability of shear deformable curved beams with non-symmetric thin-walled sections. I: Stability formulation and closed-form solutions, Comput. Struct. 83 (2005) 2525–2541. [15] M.-Y. Kim, N.-I. Kim, S.-B. Kim, Spatial stability of shear deformable curved beams with non-symmetric thin-walled sections. II: F. E. solutions and parametric study, Comput. Struct. 83 (2005) 2542–2558. [16] J. Zhu, M.M. Attard, D.C. Kellermann, In-plane nonlinear buckling of circular arches including shear deformations, Arch. Appl. Mech. 84 (2014) 1841–1860. [17] K.Y. Nieh, C.S. Huang, Y.P. Tseng, An analytical solution for in-plane free vibration and stability of loaded elliptic arches, Comput. Struct. 81 (2003) 1311–1327. [18] A.-T. Luu, J. Lee, Non-linear buckling of elliptical curved beams, Int. J. Non-Linear Mech. 82 (2016) 132–143. [19] J. Cai, Y. Xu, J. Feng, J. Zhang, In-plane elastic buckling of shallow parabolic arches under an external load and temperature changes, J. Struct. Eng. 138 (2012) 1300–1309. [20] J. Zhu, M.M. Attard, D.C. Kellermann, In-plane nonlinear buckling of funicular arches, Int. J. Struct. Stab. Dyn. 15 (2015) 1450073. [21] W. Lacarbonara, H.N. Arafat, A.H. Nayfeh, Non-linear interactions in imperfect beams at veering, Int. J. Non-Linear Mech. 40 (2005) 987–1003. [22] N. Mohamed, M.A. Eltaher, S.A. Mohamed, L.F. Seddek, Numerical analysis of nonlinear free and forced vibrations of buckled curved beams resting on nonlinear elastic foundations, Int. J. Non-Linear Mech. 101 (2018) 157–173. [23] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD finite elements NURBS exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [24] B. Hassani, A.H. Taheri, N.Z. Moghaddam, An improved isogeometrical analysis approach to functionally graded plane elasticity problems, Appl. Math. Model. 37 (2013) 9242–9268. [25] R. Bouclier, T. Elguedj, A. Combescure, On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D, Comput. Mech. 52 (2013) 1085–1112. [26] G. Zhang, K. Khandelwal, Modeling of nonlocal damage-plasticity in beams using isogeometric analysis, Comput. Struct. 165 (2016) 76–95. [27] S. Liu, T. Yu, T.Q. Bui, Size effects of functionally graded moderately thick microplates: A novel non-classical simple-FSDT isogeometric analysis, Eur. J. Mech. A Solids 66 (2017) 446–458. [28] A. Norouzzadeh, R. Ansari, Nonlinear dynamic behavior of small-scale shelltype structures considering surface stress effects: An isogeometric analysis, Int. J. Non-Linear Mech. 101 (2018) 174–186. [29] Y. Bazilevs, V.M. Calo, T.J.R. Hughes, Y. Zhang, Isogeometric fluid–structure interaction: Theory, algorithms, and computations, Comput. Mech. 43 (2008) 3–37. [30] Y. Ghaffari Motlagh, H.T. Ahn, T.J.R. Hughes, V.M. Calo, Simulation of laminar and turbulent concentric pipe flows with the isogeometric variational multiscale method, Comput. & Fluids 71 (2013) 146–155. [31] B. Bastl, M. Brandner, J. Egermaier, K. Michálková, E. Turnerová, IgA-Based solver for turbulence modelling on multipatch geometries, Adv. Eng. Softw. 113 (2017) 7–18. [32] M. Yoon, S.-H. Ha, S. Cho, Isogeometric shape design optimization of heat conduction problems, Int. J. Heat Mass Transfer 62 (2013) 272–285. [33] Z.-P. Wang, S. Turteltaub, M. Abdalla, Shape optimization and optimal control for transient heat conduction problems using an isogeometric approach, Comput. Struct. 185 (2017) 59–74. [34] Z. An, T. Yu, T.Q. Bui, C. Wang, N.A. Trinh, Implementation of isogeometric boundary element method for 2-D steady heat transfer analysis, Adv. Eng. Softw. 116 (2018) 36–49. [35] T.J.R. Hughes, J.A. Evans, A. Reali, Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems, Comput. Methods Appl. Mech. Engrg. 272 (2014) 290–320. 15

A. Hashemian and S.F. Hosseini

International Journal of Non-Linear Mechanics 113 (2019) 1–16 [73] W. Tiller Piegl, The NURBS Book, second ed., Springer-Verlag, New York, NY, 1997. [74] B.M. Imani, S.A. Hashemian, NURBS-based profile reconstruction using constrained fitting techniques, J. Mech. 28 (2012) 407–412. [75] A. Hashemian, S.F. Hosseini, An integrated fitting and fairing approach for object reconstruction using smooth NURBS curves and surfaces, Comput. Math. Appl. 76 (2018) 1555–1575. [76] R. Echter, M. Bischoff, Numerical efficiency locking and unlocking of NURBS finite elements, Comput. Methods Appl. Mech. Engrg. 199 (2010) 374–382. [77] D.J. Benson, Y. Bazilevs, M.C. Hsu, T.J.R. Hughes, Isogeometric shell analysis: The Reissner–Mindlin shell, Comput. Methods Appl. Mech. Engrg. 199 (2010) 276–289. [78] S.J. Lee, K.S. Park, Vibrations of Timoshenko beams with isogeometric approach, Appl. Math. Model. 37 (2013) 9174–9190.

[68] S.F. Hosseini, B. Moetakef-Imani, Innovative approach to computer-aided design of horizontal axis wind turbine blades, J. Comput. Design Eng. 4 (2017) 98–105. [69] A. Hashemian, S.F. Hosseini, S.N. Nabavi, Kinematically smoothing trajectories by NURBS reparameterization – an innovative approach, Adv. Robot. 31 (2017) 1296–1312. [70] A. Hashemian, B.M. Imani, Surface fairness: A quality metric for aesthetic assessment of compliant automotive bodies, J. Eng. Des. 29 (2018) 41–64. [71] A. Hashemian, B.M. Imani, A new quality appearance evaluation technique for automotive bodies including effect of flexible parts tolerances, Mech. Based Des. Struct. Mach. 46 (2018) 157–167. [72] H. Saghi, A. Hashemian, Multi-dimensional NURBS model for predicting maximum free surface oscillation in swaying rectangular storage tanks, Comput. Math. Appl. 76 (2018) 2496–2513.

16