Isogeometric Analysis of the Nonlinear Deformation of Planar Flexible Beams with Snap-back

Acta Mechanica Solida Sinica, Vol. 29, No. 4, August, 2016 Published by AMSS Press, Wuhan, China

ISSN 0894-9166

Isogeometric Analysis of the Nonlinear Deformation of Planar Flexible Beams with Snap-back⋆⋆ Zheng Huang1,2

Zeng He1,2

Wen Jiang1,2⋆

Hou Qiao1,2

Hui Wang1,2

(1 Department of Mechanics, Huazhong University of Science and Technology, Wuhan 430074, China) ( Hubei Key Laboratory for Engineering Structural Analysis and Safety Assessment, Wuhan 430074, China) 2

Received 10 October 2015, revision received 4 March 2016

ABSTRACT Based on the continuity of the derivatives of the Non-Uniform Rational B-Splines (NURBS) curve and the Jaumann strain measure, the present paper adopted the position coordinates of the control points as the degrees of freedom and developed a planar rotation-free Euler-Bernoulli beam element for isogeometric analysis, where the derivatives of the field variables with respect to the arc-length were expressed as the sum of the weighted sum of the position coordinates of the control points, and the NURBS basis functions were used as the weight functions. Furthermore, the concept of bending strip was used to involve the rigid connection between multiple patches. Several typical examples with geometric nonlinearities were used to demonstrate the accuracy and effectiveness of the proposed algorithm. The presented formulation fully accounts for the geometric nonlinearities and can be used to study the snap-through and snap-back phenomena of flexible beams.

KEY WORDS isogeometric analysis, NURBS, rotation-free beam, bending strip method, snapback

I. Introduction Flexible beams, which have a wide range of engineering applications, such as helicopter rotor blades, aircraft wings, wind-turbine blades, robot manipulators, and slender space structures[1] , are a class of slender structures that can undergo large elastic deformations. With the application of flexible structures, the nonlinear phenomena occurred in the system has attracted more and more attention. Pai and Nayfeh[2] proposed new concepts of local displacements and local engineering stresses and strains, a new interpretation and manipulation of the virtual local rotations. Their theory, which was verified by using a multiple shooting method[3] , fully accounted for large rotations, large displacements, initial curvatures, as well as extensionality, and contained most of the existing theories as special cases. This theory was further developed by using the mechanics-based variables to eliminate the discontinuous and sequence-dependent rotational variables[4] . Coda and Greco[5] presented a geometrically nonlinear formulation based on the position description for large deflection two-dimensional frame analysis. Then, Coda[6] proposed a solid-like finite element formulation without using the finite rotation schemes to solve the geometric nonlinear three-dimensional inhomogeneous frames. The description of the exact geometry of undeformed and deformed configurations is usually required in those theories, especially for naturally curved beams. Isogeometric analysis (IGA), which uses B-Splines or Non-Uniform Rational B-Splines (NURBS) instead of polynomials to describe both the exact geometry and the filed variables ⋆ ⋆⋆

Corresponding author. E-mail: [email protected] Project supported by the National Natural Science Foundation of China (Nos. 11572132 and 11572137).

· 380 ·

ACTA MECHANICA SOLIDA SINICA

2016

in analysis process, was first introduced by Hughes et al. in 2005[7], and was rapidly extended to different fields of computational mechanics, engineering and sciences[8–11] . It has already been noted that increasing the accuracy of the geometry representation entails a significant increase in the accuracy of the result[12, 13] . One of the major advantages of NURBS is that it can exactly represent all conic sections[7, 14] and there exist many efficient and numerically stable algorithms to generate the NURBS objects[15] . Another major advantage is that the concept of k-refinement allows the use of basis functions of higher continuity[14] . Several studies have been performed to develop the one-dimensional structural elements by using the concept of IGA in the last few years[16–20], where the locking phenomena were mainly concerned. Weeger et al.[21] analyzed the vibration of nonlinear Euler-Bernoulli beams with von K´arm´an strain-displacement relationship by means of the isogeometric finite elements. Based on the Green-Lagrange strain measure, Raknes et al.[22] developed an isogeometric rotation-free cable formulation. Cazzani et al.[23] presented a plane-curved beam based on the Timoshenko model and NURBS. Lan et al.[24] integrated the NURBS geometry and the rational absolute nodal coordinate formulation for the dynamic analysis of flexible bodies. This paper attempts to establish a planar rotation-free Euler-Bernoulli beam model and focuses on the snap-back phenomena of flexible structures. In §II, a flexible beam model based on the position description is presented with the Euler-Bernoulli assumption. In §III, the isogeometric finite element discretization is completed. Then, several typical examples with geometric nonlinearities are used to demonstrate the accuracy and effectiveness of the proposed formula in §IV.

II. Description of Planar Beam Structure We consider a naturally curved planar-beam, as is shown in Fig.1. Three coordinate systems are adopted to describe the initial and current configurations. The ab system is a right-hand inertial Cartesian coordinate system used for the reference configuration; the xy system is an orthogonal curvilinear coordinate system with the x axis connecting the reference point of each cross-section of the initial beam; and the ξη system is a local orthogonal curvilinear coordinate system attached to the current configuration. We let ia and ib denote the unit direction vectors of the ab coordinate system; 0 i1 and 0 i2 denote the unit direction vectors of the xy coordinate system; and t i1 and t i2 denote the unit direction vectors of the ξη coordinate system. Moreover, s ∈ [0, L] ⊂ R is the initial arc-length along the x axis from the end of the beam to the observed reference point, and L ∈ R is the initial length of the beam.

Fig. 1. Initial and current configurations of planar-beam with an arbitrary initial curvature.

2.1. Description of initial configuration The undeformed position vector 0 R of the reference line of the initial configuration is given by 0 0

R(s) = 0 X(s)ia + 0 Y (s)ib

(1)

Then, the unit tangent vector i1 of the initial reference line is defined as d0 R(s) 0 ′ 0 i1 = = X ia + 0 Y ′ ib (2) ds ′ Hereafter, without special statement, the prime (·) ≡ d(·)/ds denotes the derivative with respect to s. Using the orthogonality property of the tangent vector 0 i1 and the normal vector 0 i2 , i.e. 0 i1 · 0 i2 = 0, we obtain 0 i2 = 0 Y ′ ia − 0 X ′ ib (3)

Vol. 29, No. 4 Zheng Huang et al.: Nonlinear Deformation of Planar Flexible Beams with Snap-back

Following the Frenet-Serret formulas, we have    d 0 i1 0 = −0 k ds 0 i2

0

k 0

0

i1 0 i2



· 381 ·

(4)

0

k =0 Y ′0 X ′′ −0 X ′0 Y ′′ Here, k is the initial curvature with respect to axis x and is a function of s only.

(5)

0

2.2. Description of current configuration The position vector t R of the reference line of the current configuration is given by t

R(s) = t X(s)ia + t Y (s)ib

(6) t

The deformed arc length is denoted by sd , and then the unit tangent vector i1 of the current reference line can be obtained as t ′ dt R(s, t) ds R (s, t) dt R(s, t) t (7) = = t ′ i1 = dsd ds dsd || R (s, t)|| considering dsd /ds ≡ 1 + e(s) and p e = ||t R′ (s, t)|| − 1 = t X ′2 +t Y ′2 − 1 (8) where e(s) is a scalar that represents the elongation (axial strain) of the reference line. Similarly, using the orthogonality property of the tangent vector t i1 and the normal vector t i2 , i.e. t i1 · t i2 = 0, we can obtain t ′ t ′ Y X t i2 = ia − ib (9) 1+e 1+e Following the Frenet-Serret formulas, we have    t  t d t i1 0 k i1 = (10) t t t i − k 0 i2 ds 2 Y ′t X ′′ −t X ′t Y ′′ (1 + e)2 t where k is the current curvature with respect to axis ξ and is a function of s only. t

t

k=

(11)

2.3. Measurement of strains The Jaumann (Jaumann-Boit-Chauchy or nomial) strains are defined using the right stretch tensor of the deformation gradient. The Jaumann strains are objective engineering strains and have directions along the deformed system configuration[2] and are adopted here to describe the nonlinear deformation. The transverse shear deformation of a very thin beam is negligibly small, therefore, the cross sections of the beam can be assumed to be flat and perpendicular to the reference line before and after the deformation, which results in the so-called Euler-Bernoulli beam theory. The cross sections can be any shape and the position vector of a certain point on the observed cross-section of the initial and current configurations can be written as 0 r(s) = 0 R(s) + y 0 i2 (12) t

r(s) = t R(s) + y t i2 (13) respectively. Here, 0 R(s) and t R are the position vectors of the observed points on the undeformed and deformed reference line as defined in §2.1 and §2.2, respectively. In the Euler-Bernoulli beam theory, only the axial strain ε11 is considered, and hence the Jaumann strains can be obtained as follows: dt r(s) t d0 r(s) 0 · i1 − · i1 = e − yk, ε22 = ε12 = 0 ds ds where k = (t k − 0 k). The corresponding Jaumann stresses are obtained to be ε11 =

σ11 = Eε11 = E(e − yk) where E is Young’s modulus.

(14)

(15)

· 382 ·

ACTA MECHANICA SOLIDA SINICA

2016

2.4. Principal of virtual displacement A geometrically exact beam theory can be derived using the principal of virtual displacement, i.e. δΠ = δW

(16)

where δΠ is the variation of elastic energy, and δW is the non-conservative work due to external loads. The variation of elastic energy δΠ is given by Z LZ Z L Z T T δΠ = {δε} {σ}dAds = {δψ} [D]{ψ}ds, [D] = [S]T E[S]dA (17) 0

A

0

A



where A denotes the cross-sectional area, L is the initial beam length, [S] = 1 The variation of non-conservative work δW is written as Z L δW = (q1 δX + q2 δY + q3 δθ)ds

 −y , and {ψ} = {e, k}T . (18)

0

where q1 and q2 are distributed forces along axes a and b, respectively, q3 is the distributed bending moment, and δθ is the virtual rotation[2] defined by  t   t  δ i1 0 δθ i1 = (19) t δ t i2 −δθ 0 i2

III. Isogeometric Finite Element Method In this section, we review some basic concepts of IGA. For a detailed introduction and computational implementations, the authors refer the readers to Refs.[7,14] and [15,25]. 3.1. B-splines and NURBS In a one-dimensional parametric space, a knot vector is a set of non-decreasing coordinates, denoted as Ξ = {ξ1 , ξ2 , . . . , ξn+p+1 }, where ξi ∈ R is the i-th knot, i is the knot index, p is the polynomial order, and n is the number of basis functions used to construct the B-Spline curve. The B-Spline basis functions are defined recursively starting with piecewise constants  1 (ξi ≤ ξ < ξi+1 ) Ni,0 (ξ) = (20) 0 (otherwise) For p = 1, 2, 3, . . ., they are defined by Ni,p (ξ) =

ξi+p+1 − ξ ξ − ξi Ni,p−1 (ξ) + Ni+1,p−1 (ξ) ξi+p − ξi ξi+p+1 − ξi+1

(21)

A B-Spline curve of order p for the open knot vector Ξ with control points B i ∈ Rd , (i = 1, . . . , n) has the form n X C(ξ) = Ni,p (ξ)B i (22) i=1

Similarly, the NURBS basis functions Ri,p of order p with weights wi ∈ R, (i = 1, 2, . . . , n) for the open knot vector Ξ and the NURBS curve are given by Ni,p (ξ)wi Ri,p (ξ) = Pn ˆi=1 Nˆi,p (ξ)wˆi C(ξ) =

n X

Ri,p (ξ)B i

(23)

(24)

i=1

In this paper, we will use the NURBS basis functions and NURBS curves for the open uniform knot vectors on [0, 1] with NURBS of order p ≥ 2 and inner knots of multiplicity k = 1.

Vol. 29, No. 4 Zheng Huang et al.: Nonlinear Deformation of Planar Flexible Beams with Snap-back

· 383 ·

3.2. Discretization formula The generalized position vector P of a point on the reference line can be expressed as P (ξ) ≡ {X(ξ), Y (ξ)}T =

n X

Ri,p (ξ)B i

(25)

i=1

Considering that a control point B i only influences the curve locally in the knot interval [ξi , ξi+p+1 ), according to the element connectivity mapping[14] , one can relate the control points and basis functions to one element. Therefore, Eq.(25) associated with the jth element can be written as P (j) (ξ) =

p X

(j)

(j)

(26)

= [Φ(j) ]{δB (j) }

(27)

Ri,p (ξ)B i

i=0

Then, the variation of {ψ} is given by {δψ} = {δe, δk} =

p X

(j)

[Φ(j,i) ] · δB i

i=0

(j)

(j,i)

T T where [Φ(j) ] = {[Φ(j,0) ], . . . , [Φ(j,p) ]}, {δB (j) } = {(δB 0 )T , . . . , (δB (j) of p ) } , and the entries Φrs (j,i) [Φ ] (i = 0, . . . , p) are listed in the Appendix. Then, the variation of elastic energy δΠ can be obtained Z L δΠ = {δψ}T [D]{ψ}ds = {δB}T [K]{B} (28) 0

[K (j) ]{B (j) } =

Z

[Φ(j) ]T [D]{ψ (j) }ds

(29)

L(j)

where Ne is the total number of elements, L(j) is the length of jth element, [K (j) ] is the element stiffness matrix, {B (j) } is the element generalized position vector of control points, [K] and {B} are the global stiffness matrix and the global generalized position vector of control points, respectively. The variation of non-conservative work due to external distributed loads is given by Z L δW = (q1 δX + q2 δY + q3 δθ) ds = {δB}T {R} (30) 0

where {R} is the global loading vector, and the elemental loading vector {R(j) } is given by Z Z p X {R(j) } = [F (j,i) ] · {Rc(j,i) }T ds = [F (j) ]{Rc(j) }ds L(j) i=0

(j,i)

(j,i)

(j,i)

(j,i)

(31)

L(j)

(j,i)

in which {Rc } = {q1 , q2 , q3 }T , and the entries Frs of [F (j,i) ] (i = 0, . . . , p) are listed in the Appendix. Substituting Eqs.(28) and (30) into Eq.(16), one can obtain the following equation [K]{B} = {R}

(32)

Here, both [K] and {R} are the nonlinear functions of {B}. In the next subsection, Eq.(32) will be solved by using an incremental/iterative method based on the modified Riks method[26] . 3.3. Incremental/Iteration method To derive the incremental form of Eq.(32) for static analysis using the incremental/iterative method, we define ¯ + {∆B (j) } {B (j) } = {B} (33) (j) ¯ denotes an equilibrium state and {∆B } denotes the generalized position vector when the where {B} loads continue to increase. Substituting Eq.(33) into {ψ} and [Φ] yields (j) ¯ + [Φ]{∆B ¯ {ψ (j) } = {ψ} }

(34)

· 384 ·

ACTA MECHANICA SOLIDA SINICA

2016

¯ + [Ψ ]{∆B (j) } [Φ(j) ] = [Φ] (j,i)

(0) (p) [Ψ ] = {[Ψrs ], . . . , [Ψrs ]},

(i) Ψrs =

(35) (j,i)

∂Φmr

Dmn ψ¯n = (j)

∂B i,s

∂ 2 ψ¯m (j)

(j)

∂B i,r ∂B i,s

Dmn ψ¯n

Then, substituting Eqs.(33), (34), and (35) into Eq.(29) yields  (j)  ⌣ (j) (j) (j) (j) [K ]{B } = [K ]{B } (j) + K {∆B (j) } ¯ 

⌣ (j)

where K



(37)

}={B}

{B

(36)

is the elemental tangent stiffness matrix and is given by ⌣ (j)

[K

]=

Z

¯ T [D][Φ] ¯ + [Ψ ])ds ([Φ]

(38)

L(j)

[K (j)

(j)

]{B

(j)

(j)

}

¯ {B (j) }={B}

=

Z

¯ ¯ T [D]{ψ}ds [Φ]

(39)

L(j)

¯ c } + {∆Rc } into Eq.(31), we can obtain Substituting {Rc } = {R Z (j) (j) ¯ + [f (B)]{∆R ¯ ¯ {R(j) } = [F (j) ]{Rc(j) }ds = {R} } c } + [g(B)]{∆B

(40)

L(j)

¯ = [f (B)]{ ¯ R ¯ c }, {R}

¯ = [f (B)]

Z

L(j)

[F (j) ]ds,

¯ = [g(B)]

h  i (j) (j) ∂frs /∂B i,s · Rc,s

¯ {B (j) }={B}

(41)

Furthermore, substituting Eqs.(37) and (40) into Eq.(32), we can obtain the following incremental governing equation:   ˆ − [g(B)] ¯ ¯ [K] {∆B} = [f (B)]{∆R (42) c} ¯ is very tedious and unnecessary, Actually, the derivation of the specific expressions of [Ψ ] and [g(B)] so it will be calculated numerically using the Richardson extrapolation method.

IV. Numerical Examples R The isogeometric finite element formulation presented above has been implemented by using MATLAB and taking advantage of the program package on NURBS provided by Spink[15] . In order to verify the proposed formulation, five numerical examples will be given in this section. The first two are used to verify the accuracy of the presented algorithm, and then three representative examples involving buckling and snap-back will be investigated to demonstrate the effectiveness of the presented formulation. Both the initial straight and curved beam structures have been studied and the large deformation can also be simulated.

4.1. Cantilever subjected to an end moment In this classical Euler beam problem, an initial horizontal cantilever is bent into a circular arc by a bending moment at the free end. Here, L, I, and E are the initial length of the beam, the moment of inertia, and Young’s modulus, respectively. The dimensionless bending moment of the free end is defined by n = M L/(2πEI), where M represents the bending moment actually applied. Ten elements with the polynomial order of the NURBS curve p = 3 and C 2 -continuous basis functions are used to discrete the beam. The analytical solution formula can be written as L L sin θ, Y |s=L = (1 − cos θ) (43) θ θ where θ is the rotation angle of the free end of the beam. Figure 2(a) shows the horizontal and vertical positions of the free end varying with the rotation angle of the free end and compares the present numerical results with the analytical solution. The deformed configurations of the cantilever corresponding to different load levels are shown in Fig.2(b). X|s=L =

Vol. 29, No. 4 Zheng Huang et al.: Nonlinear Deformation of Planar Flexible Beams with Snap-back

· 385 ·

Fig. 2. (a) Comparison of the horizontal and vertical positions of the free end of cantilever obtained from the IGA formulation with the corresponding analytical results; (b) Equilibrium configurations under different free end moment obtained from the IGA formulation.

Fig. 3. Convergence of the presented method with p = 3 and p = 4. The normalized errors of the Euclidean norm versus number of elements.

Next we set n equal to 0.25, and then we select p = 3 and p = 4 as the NURBS orders in the numerical test. Convergence plots of the configuration of the cantilever, normalized with respect to the exact solution, are reported in Fig.3 for different NURBS orders. One can see that the high order continuity of NURBS basis functions has shown distinct advantages in the aspect of convergence. As shown in Fig.3, the proposed isogeometric finite element method provides a good approximation of the solution also in the case of relatively coarse meshes.

4.2. Cantilever subjected to a transverse end load We consider the same cantilever as in example 4.1, which, however, is subjected to a transverse 1 shear end load F = π 2 EI/L2 . Table 1 reports the X and Y position at the free end of the cantilever 8 for various NURBS orders p and element numbers ne . In the case of low NURBS order, few elements can get good precision. Good results can be obtained with the NURBS order increased, even though only one element is used in the calculation. The normalized absolute error curves of horizontal and vertical positions varying with arc length shown in Fig.4 demonstrate that the accuracy of the presented algorithm is high enough. In fact, ten 2-node finite elements can also achieve the same values of X and Y positions at the free end. However, the accuracy of the non-nodal value is reduced. Even in this case, the total number of degrees of freedom of the presented method (p = 10, ne = 1) is less than that of ANSYS with ten 2-node finite elements, while the decrease of the amount of accuracy is relatively negligible.

· 386 ·

ACTA MECHANICA SOLIDA SINICA

2016

Table 1. Cantilever subjected to a transverse end load: X and Y positions at the free end versus the NURBS order p and element number ne

End Position X/L Y /L 3 1 0.9460 0.2816 3 2 0.9178 0.3576 3 4 0.9201 0.3563 IGA 4 1 0.9283 0.3359 6 1 0.9201 0.3565 10 1 0.9201 0.3567 Reference∗ 0.9021 0.3565 ∗ The reference solution is obtained from the ANSYS with one hundred 2-node finite elements (Beam188). p

ne

Fig. 4. Comparison of the deformed configuration obtained from IGA and the reference results obtained from the finite element commercial software ANSYS.

4.3. Buckling of Euler beam We consider the same cantilever as in example 4.1, which, however, is subjected to an axial load. The proposed algorithm can thus be employed to calculate large deformation. The beam is divided into twenty elements; and other parameters adopted here are the same as those in the previous example, but an initial imperfection is required to start the calculation. We assume that there is an initial linear distribution transverse deflection along the axial direction, where the deflection of the free end is L/1000. The load level is defined as n = −F/Fcr , where Fcr = π 2 EI/(4L2 ) is the Euler critical load for the first buckling mode, and F is the axial load actually applied. As demonstrated in Fig.5, because of

Fig. 5. The first buckling mode of the cantilever.

Vol. 29, No. 4 Zheng Huang et al.: Nonlinear Deformation of Planar Flexible Beams with Snap-back

· 387 ·

the initial imperfect, when n varies from 0 to 1, although the lateral deflection is relatively large, the axial deformation is still extremely small. However, both the horizontal and vertical displacements are rapidly growing larger with only a small increase of the load if n > 1. According to the load-position curves shown in Fig.5(a), one can obtain that the first order buckling mode is very close to the critical load when n = 1.

Fig. 6. Lee’s frame subjected to a concentrated load.

4.4. Lee’s frame with snap-back The Lee’s Frame[27] is depicted in Fig.6. The frame subjected to a concentrated load is widely used as a benchmark problem because it shows a distinguished snap-back behavior, i.e. the load drop phenomenon, which means that when the displacement reaches the maximum, the load will fall quickly. Points A and B are two hinges, C is a rigid connection, and D is the loading point. For the rotation-free beam formula, in order to maintain the angle of the two tangent vectors and thereby fix the rotations at junction C, a curved bending strip (a curved beam element) should be applied at the junction, as suggested by Kiendl et al.[28] . This bending strip extends from one control point of the AC patch to another on the BC patch via the common control point, and it is a quadrant with the order p = 3 and the inner knots of multiplicity k = 1, as is shown in Fig.6. Both the horizontal displacement U and the vertical displacement V of the loading point D are plotted in Fig.7(a), and the corresponding configurations of three stages, marked with points a, b, and c in Fig.7(a), are shown in Fig.7(b). The elastic critical load obtained from the presented method is Fa = 18.52EI/L2 , which has a difference of 0.13% when compared with the result of Fa = 18.55EI/L2 in Lee et al.[27] . The snap-back will occur in the vertical displacement V beyond point a; while the horizontal displacement U will decrease beyond point b when the load increases. The corresponding critical load Fb = −9.476EI/L2 is close to the

Fig. 7. (a) The horizontal and vertical displacements of the loaded point D and three points a, b, and c on the equilibrium path; (b) the configurations corresponding to points a, b, and c on the equilibrium path.

· 388 ·

ACTA MECHANICA SOLIDA SINICA

2016

result Fb = −9.621EI/L2 in Fujii[29] . 4.5. DaDeppo’s arch A pinned-clamped circular arch subjected to a concentrated load F on crown E is shown in Fig.8. Since only an arc with the central angle less than 180◦ can be constructed from a single quadratic NURBS element with the order p = 2 and the inner knots of multiplicity k = 1, the present arc should be constructed from multiple smaller arcs; however, the basis function can be no more than C 0 -continuous in the junction of those smaller arcs[14] . In the current case, the connection control point of two adjacent arcs and its two neighbor control points should be kept on a straight line in order to maintain the continuity of the bending moment on the junction point. The central angle of the arch is 215◦ and will be divided into three arcs with the central angles of 90◦ , 35◦ , and 90◦ , respectively. The moment of the connection point of two adjacent arcs is transferred by using the bending strip method. Different from the previous example, the bending strip (a straight beam element) is straight in this example, as is shown in Fig.8. Both the horizontal displacement U and the vertical displacement V of the loaded point E are plotted in Fig.9(a), and the corresponding configurations of three stages on the equilibrium path, marked with a, b, and c, are shown in Fig.9(b). The vertical displacement V shows a relatively large rigidity before point a, and the horizontal displacement U demonstrates even stiffer as the load increases. However, the DaDeppo’s arch exhibits instability without softening beyond point a. The difference between the critical load of Fa = 9.197EI/R2 obtained from the presented method and the result of Fa = 9.122EI/R2 obtained in Fujii[29] is only 0.82%. The critical load of point b is Fb = 5.274EI/R2 . As can be seen from Fig.9(b), the structure will be self-contact after point b, which cannot happen in the real world. This example further describes the snap-back occurring in the initial curved beam structures.

Fig. 8. DaDeppo’s arch subjected to a concentrated load.

Fig. 9. (a) The horizontal and vertical displacements of the loading point E; and (b) the configurations corresponding to stages a, b, and c on the equilibrium path in (a).

Vol. 29, No. 4 Zheng Huang et al.: Nonlinear Deformation of Planar Flexible Beams with Snap-back

· 389 ·

V. Conclusions In the present paper, a planar rotation-free Euler-Bernoulli beam model, which fully accounts for the geometric nonlinearities, is developed following the principal of virtual displacement based on the isogeometric finite element discretization. The presented formula uses the position coordinates of the control points as the field variables and adopted the Jaumann strain measure. In order to consider the rigid connection of the structures, the bending strip, both a curved connecting element with the order p = 3 and the inner knots of multiplicity k = 1 and a straight connecting element with the order p = 2 and the inner knots of multiplicity k = 1 are introduced. The numerical examples demonstrate that the presented formula and nonlinear solution technique are very accurate and efficient. The beam formula proposed in this paper is simple and easy to be extended to derive the cable or other beam theories. In addition, the high order continuous property of NURBS can also be applied to the isogeometric collocation method which is a new approach developed recently.

References [1] Pai,P.F., Highly Flexible Structures: Modeling, Computation and Experimentation. Reston, Virginia: AIAA, 2007. [2] Pai,P.F. and Nayfeh,A.H., A fully nonlinear theory of curved and twisted composite rotor blades accounting for warping and three-dimensional stress effects. International Journal of Solids and Structures, 1994, 31(9): 1309-1340. [3] Pai,P.F., Large-deformation analysis of flexible beams. International Journal of Solids and Structures, 1996, 33(9): 1335-1353. [4] Pai,P.F., Geometrically exact beam theory without Euler angles. International Journal of Solids and Structures, 2011, 48: 3075-3090. [5] Coda,H.B. and Greco,M., A simple FEM formulation for large deflection 2D frame analysis based on position description. Computer Methods in Applied Mechanics and Engineering, 2004, 193: 3541-3557. [6] Coda,H.B., A solid-like FEM for geometrically non-linear 3D frames. Computer Methods in Applied Mechanics and Engineering, 2009, 198: 3712-3722. [7] Hughes,T.J.R., Cottrel,J.A. and Bazilevs,Y., Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194: 41354195. [8] Cottrel,J.A., Reali,A., Bazilevs,Y. and Hughes,T.J.R., Isogeometric analysis of structure vibrations. Computer Methods in Applied Mechanics and Engineering, 2006, 195: 5257-5259. [9] Bazilevs,Y., Calo,V.M., Tezduyar,T.E. and Hughes,T.J.R., Isogeometric fluid-structure interaction: theory, algorithms, and computations. Computational Mechanics, 2008, 43: 3-37. [10] Gomez,H., Calo,V.M., Bazilevs,Y. and Hughes,T.J.R., Isogeometric analysis of Hilliard phase-filed model. Computer Methods in Applied Mechanics and Engineering, 2009, 198: 1726-1741. [11] Wall,W.A., Frenzel,M.A. and Cyron,C., Isogeometric structural shape optimization. Computer Methods in Applied Mechanics and Engineering, 2010, 199: 290-300. [12] Ibrahimbegovic,A., On finite element implementation of geometrically nonlinear Reissner’s beam theory: three-dimensional curved beam elements. Computer Methods in Applied Mechanics and Engineering, 1995, 122: 11-26. [13] Benson,D.J., Bazilevs,Y., Hsu,M.-C. and Hughes,T.J.R., Isogeometric shell analysis: the Reissner-Mindlin shell. Computer Methods in Applied Mechanics and Engineering, 2010, 199: 276-289. [14] Cottrel,J.A., Hughes,T.J.R. and Bazilevs,Y., Isogeometric Analysis: toward integration of CAD and FEA. John Wiley & Sons, 2009. [15] Spink,M. NURBS toolbox, http://www.aria.uklinux.net/nurbs.php3, http://octave.sourceforge. net/nurbs/index.html. [16] Echter,R. and Bischoff,M., Numerical efficiency, locking and unlocking of NURBS finite elements. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5-8): 374-382. [17] Bouclier,R., Elguedj,T. and Combescure,A., Locking free isogeometric formulations of curved thick beams. Computer Methods in Applied Mechanics and Engineering, 2012, 245-246: 144-162. [18] Beir˜ ao da Veiga,L., Lovadina,C. and Reali,A., Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods. Computer Methods in Applied Mechanics and Engineering, 2012, 241-244: 38-51. [19] Cuomo,M. and Greco,L., Isogeometric analysis of space rods: Considerations on stress locking. In: ECCOMAS 2012, 2012: 5094-5112.

· 390 ·

ACTA MECHANICA SOLIDA SINICA

2016

[20] Greco,L. and Cuomo,M., B-spline interpolation of Kirchhoff-Love space rods. Computer Methods in Applied Mechanics and Engineering, 2013, 256: 251-269. [21] Weeger,O., Wever,U. and Simeon,B., Isogeometric analysis of nonlinear Euler-Bernoulli beam vibrations. Nonlinear Dynamics, 2013, 72: 813-835. [22] Raknes,S.B., Deng,X., Bazilevs,Y., Benson,D.J., Mathisen,K.M. and Kvamsdal,T., Isogeometric rotationfree bending-stabilized cables: Statics, dynamics, bending strips and coupling with shells. Computer Methods in Applied Mechanics and Engineering, 2013, 263: 127-143. [23] Cazzani,A. and Malagu,M., Isogeometric analysis of plane-curved beams. Mathematics and Mechanics of Solids, 2014, 1081286514531265. [24] Lan,P., Yu,Z., Du,L. and Lu,N. Integration of non-uniform rational B-splines geometry and rational absolute nodal coordinates formulation finite element analysis. Acta Mechanica Solida Sinica, 2014, 27(5): 486-495. [25] De Falco,C., Reali,A. and V´ azquez,R., GeoPDEs: a research tool for Isogeometric Analysis of PDEs. Advances in Engineering Software, 2011, 42(12): 1020-1034. [26] Riks,E., An incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures, 1979, 15: 524-551. [27] Lee,S.L., Manuel,F.S. and Rossow,E.C., Large deflections and stability of elastic frames. ASCE Journal of Engineering Mechanics, 1968, 504(94): 521-547. [28] Kiendl,J., Bazilevs,Y., Hsu,M.C, Wuchner,R. and Bletzinger,K.U., The bending strip method for isogeometric analysis of Kirchhoff-Love shell structures comprised of multiple patches. Computer Methods in Applied Mechanics and Engineering, 2010, 199: 2403-2416. [29] Fujii,F., Scheme for elastics with snap-back and looping. ASCE Journal of Engineering Mechanics, 1989, 115: 2166-2181.

Appendix The variation of axial strain δe, and curvature δ t κ are given by t

δe =

t ′ X′ t ′ Y t ′ δX + δY 1+e 1+e

δtk =

t

(44)

t ′ t ′′ XδY X ′′ − 2t k t Y ′ t ′ t Y ′′ + 2t k t X ′ t ′ t Y ′ δ t X ′′ δ Y − δ X + − (1 + e)2 (1 + e)2 (1 + e)2 (1 + e)2

(45)

By using Eq.(19), one can obtain δθ = δ t i1 ·t i2 =

t

t ′ Y′ X δtX ′ − δtY ′ 2 (1 + e) (1 + e)2

(46)

For different integers i (i = 0, . . . , p), the only difference of [Φ(j,i) ] is their associated NURBS basis (j,i) functions. According to the definition of the entries Φrs of [Φ(j,i) ], (j,i)

′ = t T11 Ri,p ,

(j,i)

=

Φ11 Φ21

(j,i)

Φ12

′ = t T12 Ri,p

t ′′ Y + 2t k t X ′ ′ Y′ ′′ R Ri,p , − i,p (1 + e)2 (1 + e)2 t

(47) (j,i)

Φ22

t

=

t ′ X ′′ − 2t k t Y ′ ′ X R R′′ − i,p 2 (1 + e) (1 + e)2 i,p

(48)

Here, one should notice that the index i of NURBS basis functions Ri,p and their derivatives are local ′ ′′ ′′′ numbers. Moreover, Ri,p is defined in the parametric space identified by ξ, while Ri,p , Ri,p , and Ri,p are the first, second, and third derivatives of NURBS basis functions with respective to the arc-length s, respectively. According to the chain rule, one can yield dRi,p (ξ) dξ dRi,p (ξ) = ds dξ ds   2 d dR (ξ) dξ d2 Ri,p (ξ) dξ dRi,p (ξ) d2 ξ i,p ′′ Ri,p (ξ) ≡ = + ds dξ ds dξ 2 ds dξ ds2

′ Ri,p (ξ) ≡

where dξ/ds and d2 ξ/ds2 can be obtained by geometry mapping (24). (j,i) Similarly to [Φ(j,i) ], the entries Frs of [F (j,i) ] are given by   t ′ Y F0 0 F1 [F (j,i) ] = , F0 = Ri,p , F1 = R′ , 0 F0 F2 (1 + e)2 i,p

(49)

t

F2 = −

X′ R′ (1 + e)2 i,p

(50)