NURBS-based isogeometric vibration analysis of generally laminated deep curved beams with variable curvature

Accepted Manuscript NURBS-based isogeometric vibration analysis of generally laminated deep curved beams with variable curvature Anh-Tuan Luu, Nam-Il Kim, Jaehong Lee PII: DOI: Reference:

S0263-8223(14)00404-8 http://dx.doi.org/10.1016/j.compstruct.2014.08.014 COST 5847

To appear in:

Composite Structures

Please cite this article as: Luu, A-T., Kim, N-I., Lee, J., NURBS-based isogeometric vibration analysis of generally laminated deep curved beams with variable curvature, Composite Structures (2014), doi: http://dx.doi.org/10.1016/ j.compstruct.2014.08.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

NURBS-based isogeometric vibration analysis of generally laminated deep curved beams with variable curvature Anh-Tuan Luu,1, ∗ Nam-Il Kim,1, † and Jaehong Lee1, ‡ 1

Department of Architectural Engineering, Sejong University, 98 Kunja Dong, Kwangjin Ku, Seoul 143-747, Republic of Korea

In this paper, the NURBS-based isogeometric analysis is developed for the free vibration analysis of the generally laminated Timoshenko-type of deep curved beams with arbitrary curvature. The equivalent modulus of elasticity is utilized to account for the material couplings of the beam and the deepness term is exactly integrated into ABD parameters of the curved beam. The geometry and the curvature of the laminated free-form curved beams are modeled by NURBS so that the gap between the analysis of the laminated curved beam with constant curvature and the analysis of that with variable curvature is bridged. All the effects of the axis extensibility, the shear deformation and the rotary inertia are taken into account. Results of the non-dimensional frequencies for the laminated elliptic and circular curved beams with different lay-ups are compared with other available results in order to demonstrate the validity of the present isogeometric model. The free vibration analysis of the laminated parabolic curved beams and the laminated elliptic ring are presented. Particularly, the laminated Tschirnhausen’s cubic curved beam is considered to study the dynamic responses as an example of the laminated free-form curved beam.

Keywords: Isogeometric analysis; Laminated composite; Free vibration; Curved beam; Variable curvature.

1. Introduction The laminated composite materials are finding increased usage in a wide variety of structural applications in the aerospace, civil construction, marine and offshore industries during the past four decades. The curved beam made with the laminated composite materials are generally used as structural components of light-weight heavy load bearing elements because of the high strength-to-weight and stiffness-to-weight ratios, the ability of being different strengths in different directions and the nature of being tailored to satisfy the design requirements of strength and stiffness in practical designs. For a curved beam structure that may be subjected to dynamic loads, the determination of the natural frequencies and corresponding mode shapes is critical in the design process. It is usually the first step in a dynamic analysis since a great deal may be deduced concerning the structural behavior and integrity from the knowledge of its natural frequencies and mode shapes. Up to the present, considerable research efforts have been made for the improved free vibration analysis of curved beams made with composite materials. Qatu [1] and Qatu and Elsharkawy [2] applied the classical laminated curved

∗ Graduate

student professor ‡ Professor, Corresponding author. Tel.:+82-2-3408-3287, fax:+82-2-3408-3331, E-mail address: [email protected] † Research

2 beam theory to the vibration problems of thin arches by using the Ritz method for the simply supported and arbitrary boundaries, respectively. For a moderately deep laminated arch with highly anisotropic materials, the natural frequencies may be overestimated since the effects of shear deformation and rotary inertia are neglected. To take into account these effects, Qatu [3] developed the first order shear deformation theory for moderately deep arches including the deepness term in the fundamental equations. Ahmed [4,5] studied the effect of transverse shear deformation for honeycomb sandwich curved beams using the finite element analysis. Raveendranath [6] derived coupled polynomials displacement fields using element equilibrium equations incorporate geometrical and material couplings in a consistent manner and proposed an efficient two-node Timoshenko-type laminated curved beam element. Malekzadeh [7] presented the two-dimensional theory of elasticity for the free vibration analysis of thick laminated deep circular arches. The layerwise finite element method was used to obtain the through-the-thickness discretized forms of the equations of motion and exact solution was obtained for the simply supported antisymmetric cross-ply laminated arches. Malekzadeh et al. [8] introduced an efficient hybrid layerwise and differential quadrature method based on the two-dimensional elasticity theory. To approximate the displacement components in the radial direction, the layerwise theory was used and the resulting equations of motion were discretized using the differential quadrature method. They demonstrated the fast rates of convergence of the method and its high accuracy with low computational efforts. The equations of motion were derived by Khdeir and Reddy [9] from the Hamilton’s principle for the free and forced vibration of the cross-ply laminated thin, moderately thick, and thick shallow arches. The exact natural frequencies were determined for various end conditions using the state space concept. Recently, Hajianmaleki and Qatu [10] developed a modified first order shear deformation theory accounting for deepness, laminate coupling, shear deformation, and rotary inertia for static and free vibration analyses of composite curved beams. They presented the exact solution for simply supported boundary condition as well as numerical solutions using the general differential quadrature method for other boundary conditions. However, the above mentioned studies are limited for the vibration analysis of laminated curved beams with constant curvature. Curved beams with variable curvature used in many engineering applications show the complex structural behavior since deformations of beam depend on the coupled equations between tangential, radial displacements and rotation by curvature effects. A number of investigations have been carried out to study the free vibration behavior of these types of curved beams made of isotropic materials. Oh et al. [11] derived the governing equations to determine natural frequencies and mode shapes for the shear deformable non-circular curved beam. They used the determinant search method combined with the Regular-Falsi method to calculate the natural frequency. Oh et al. [12] also performed the free vibration analysis of curved beams with variable cross-section considering the axial extensibility and the rotary inertia effects, but neglected the shear deformation. Some researchers used polynomials or the power series expansion to approximate displacement field. Huang et al. [13] and Tseng et al. [14] developed the dynamic stiffness method based on the Timoshenko curved beam theory which considers all the effects of axis extensibility, shear deformation, and rotary inertia to study the free vibration of the non-circular curved beams. They decomposed the arch into several subdomains and the series solution of each subdomain was formulated in terms of polynomials, and then solved inplane vibration of laminated arches. Nieh et al. [15] studied the vibrational behvior of elliptic arches by using the subdomains concept to develop the displacement field and Rossi et al. [16] approximated the tangential displacement by using polynomials to solve the in-plane vibration problem of the cantilevered non-circular arches of non-uniform cross-section with a tip mass. Moreover, some used the suitable trial functions to approximate the displacement

3 field. For example, Romanelli and Laura [17] used trial functions of tangential and radial displacements satisfying boundary conditions to solve the fundamental frequency of the hinged arches and Wang and Moore [18] utilized assumed displacements satisfying the elliptic arch with clamped boundary condition to solve the lowest extensional natural frequency. Tseng et al. [19] incorporated the dynamic stiffness method and the series solution to acquire an analytical solution to the free vibration of laminated semi-elliptic curved beam. It is well known that the finite element method (FEM) has been widely used as one of the most versatile and powerful techniques in the free vibration analysis of curved beams. But it has many inherent disadvantages including, for instance, complex equation-based calculations of the curvature, the cumbersome task of mesh generation; a higher order is not easy task to construct conventionally conformable curved beam elements; a time-consuming procedure for the connectivity of elements; re-meshing in moving boundary problems and so on. Therefore, to overcome these drawbacks, as an alternative to FEM, a family of the so-called Isogeometric analysis has been recently introduced by Hughes et al. [20]. The isogeometric analysis (IGA) has been a powerful computational method that has the ability of integrating the finite element analysis (FEA) and the Computer Aided Design (CAD). The fundamental concept behind IGA is to utilize the basis functions that are able to model exactly the geometry from the CAD point of view for numerical simulations of physical phenomena. It can be basically achieved by using the B-splines or NURBS (Non-Uniform Ration B-Splines) for the geometrical representation and invoking the isoparametric concept to define the unknown field variables. A distinct advantage over FEM is that the mesh refinement is not only simply accomplished by a automatic communication with the CAD geometry tools, but leaves the geometry intact. Another intriguing trait of these functions is that they are typically smooth beyond the classical C 0 -continuity of the standard FEM. Comprehensive knowledge of IGA could be found in the text book of Cottrell et al. [21]. This IGA-based approaches have been used to solve many dynamics problems in a wide range of research areas. Realy [22] and Cottrell et al. [23] presented the free vibration study of Euler-Bernoulli beams and Kirchhoff-Love plates by using exact geometrical isogeometric models. Later, Shojaee et al. [24] and Thai et al. [25] implemented NURBS-based isogeometric formulation to investigate the natural frequencies and mode shapes of laminated composite plates governed by Kirchhoff-Love and Reissner-Mindlin theory, respectively. The free vibration of Timoshenko straight beams with the isogeometric approach was carried out by Lee and Park [26]. Nagy et al. [27] implemented a NURBS-based control mesh to find the optimal shape for the maximum natural frequency of Euler-Bernoulli curved beams. From the previously cited references, one can note that despite extensive research for the vibration analysis of the non-circular curved beams made of isotropic materials, to the authors’ knowledge, Ref. [19] is the only study reported on the free vibration analysis of laminated composite curved beams with variable curvature. Besides, it is well known that NURBS have the ability to model free-form curves accurately while the previous literature is based on the complex separate equations to represent the length and the curvature of curved beams. Therefore, that gives us a strong motivation to obtain the accurate natural frequencies and corresponding mode shapes of laminated composite free-form curved beams using the NURBS-based isogeometric analysis. In this paper, the NURBS-based isogeometric finite element analysis is modelled based on the refined Timoshenkotype of laminated curved beam theory using the equivalent modulus of elasticity and modified ABD parameters which account for the material coupling and the deepness term, respectively. The NURBS basis functions are used instead of traditional Lagrange interpolation functions, therefore, geometries and curvatures of the curved beams can be modelled

4 accurately. Furthermore, the h-, p- and k-refinement strategies are implemented to create the NURBS elements with high orders and high continuities. The verification of the present model is obtained by comparing the current results with previously published works for the laminated elliptic and circular curved beams and excellent agreements are observed. Then, laminated parabolic curved beams, laminated full elliptic rings, and laminated Tschirnhausen’s cubic curved beams with various lay-ups, rise to span length, slenderness and orthotropy ratios are taken into account for the investigation of free vibration behaviors.

2. Laminated curved beam theory 2.1. Kinematic relationship In this study, the laminated composite beam with rectangular cross-section of width b and thickness h is characterized by its middle surface. The beam is composed of layers of orthotropic material oriented at arbitrary angle θ with respect to the longitudinal axis. Fig. 1 shows the geometry and positive sign conventions of the laminated curved beam with variable curvature in a right-handed curvilinear coordinate system (x, z), in which x is the curvilinear coordinate coincides with the centroidal axis and z is the coordinate perpendicular to the center line in-plane of curved beam; κ is the curvature; u and v are the tangential and radial displacements, respectively, of the beam middle surface material point; ψ is the bending rotation of the beam cross section about the out-of-plane axis; N is the axial force, M and Q are the bending moment and the shear force, respectively. In the displacement-based formulation, the assumed displacement field is presented as: u ¯(x, z, t) = u(x, t) + zψ(x, t)

(1a)

w(x, ¯ z, t) = w(x, t)

(1b)

where u ¯ and w ¯ are the generalized displacements along the x and z directions, respectively; t is the time. The middle surface strains and the curvature change are given by o = u + κw 

(2)

γ = w + ψ − κu

(3)

χ = ψ

(4)

where o is the middle surface strain, γ is the shear strain at the neutral axis and χ is the curvature change. The superscript “prime” refers to differentiation with respect to x. The normal strain at an arbitrary point can be found from =

1 (o + zχ) 1 + κz

(5)

5 2.2. Generalized constitutive relations The force and moment resultants are the integrals of the stresses over the beam thickness as follows:  h/2 N = b σdz

(6a)

−h/2

 M = b

−h/2

 Q = b

h/2

σzdz

(6b)

τ dz

(6c)

h/2

−h/2

where σ and τ are the normal and shear stresses, respectively. The above equations may be rewritten as: m  zk  N = b σdz M = b Q = b

k=1 zk−1 m  zk  k=1 zk−1 m  zk  k=1

(7a)

σzdz

(7b)

τ dz

(7c)

zk−1

where m is the number of layers and zk and zk−1 are the z coordinates of the top and bottom of the kth layer as shown in Fig. 1. It is well known that the issue of couplings appear in laminates with arbitrary lay-up, thus the in-plane and out-ofplane vibrations need to be solved together. However, Hajianmaleki and Qatu [10] showed that the purely in-plane free vibration behavior of laminated curved beams can be accurately achieved with any of stacking sequence by using equivalent modulus of elasticity of each lamina. The equivalent modulus of elasticity of each lamina is found based on the following equation: 1 (k)

Ex (k)

where Ex

=

cos4 (θk ) + E11



1 2ν12 − G12 E11



cos2 (θk )sin2 (θk ) +

sin4 (θk ) E22

(8)

is the equivalent modulus; E11 and E22 are elastic moduli along directions parallel and perpendicular to

fiber, respectively; G12 and ν12 are the shear modulus and the Poisson ratio, respectively. The stress-strain relationship for any layer of the laminate can be written as: σ = Ex(k) 

(9a)

¯ (k) f (z)γ τ = Q 55

(9b)

¯ 55 is the shear modulus and f (z) is the shear correction function. Assuming parabolic distribution through where Q the thickness, the shear correction function is given by [28]   2  z 5 1− f (z) = 4 h/2

(10)

Substituting Eqs. (9a) to (9b) into Eqs. (7a) to (7c), utilizing Eq. (5) and carrying out the integration over the thickness piecewise, from layer to layer, yields: ⎫ ⎡ ⎧ ⎪ ⎪ N ⎪ ⎪ ⎬ ⎢ A11 B11 0 ⎨ =⎢ M ⎣ B11 D11 0 ⎪ ⎪ ⎪ ⎪ ⎭ ⎩ Q 0 0 A55

⎤⎧ ⎫ ⎪ ⎪ ⎪ε⎪ ⎥⎨ ⎬ ⎥ χ ⎦⎪ ⎪ ⎪ ⎭ ⎩ ⎪ γ

(11)

6 where the A11 , B11 , D11 and A55 are the stiffness coefficients arising from the piecewise integration:   m 1 + κzk 1  (k) A11 = bEx ln κ 1 + κzk−1 k=1    m 1 + κzk 1 1  (k) B11 = bEx (zk − zk−1 ) − ln κ κ 1 + κzk−1 k=1      2   m 2 1 1 1 1 + κzk 1  (k) 1 2 + zk − + zk−1 D11 = bEx − (zk − zk−1 ) − 2 ln κ 2 κ κ κ κ 1 + κzk−1 k=1   m   4 5  ¯ (k) 3 A55 = bQ55 (zk − zk−1 ) − 2 ln zk3 − zk−1 4 3h

(12a) (12b)

(12c) (12d)

k=1

2.3. Equations of motion The strain energy stored in a beam during elastic deformation is  1 L (N o + M χ + Qγ) dx U = 2 0 Writing the strain energy functional for the kth lamina, and summing for the whole laminate yields   1 L U = A11 2o + 2B11 o χ + D11 χ2 + A55 γ 2 dx 2 0

(13)

(14)

Substituting Eqs. (2), (3), and (4) into the above energy expression yields the strain energy functional in terms of the slope and displacements   1 L 2 2 2 U = A11 (u + κw) + 2B11 (u + κw) ψ  + D11 (ψ  ) + A55 (w + ψ − κu) dx 2 0 The kinetic energy for the entire beam including the rotary inertia is expressed from the study of Qatu [3]     2 b L zk T = ρ (1 + κz) u˙ + z ψ˙ + w˙ 2 dzdx 2 0 zk−1

(15)

(16)

where ρk is the mass density and the over “dot” denotes differentiation with respect to t. Integrating over z and summing for m number of layers yields    1 L  I1 + 2κI2 + κ2 I3 u˙ 2 + I1 w˙ 2 + 2 (I2 + κI3 ) u˙ ψ˙ + I3 ψ˙ 2 dx T = 2 0

(17)

where I1 = I2 = I3 =

m  k=1 m  k=1 m  k=1

bρ(k) (zk − zk−1 )

(18a)

bρ(k)

 1 2 2 zk − zk−1 2

(18b)

bρ(k)

 1 3 3 zk − zk−1 3

(18c)

The work done by the distributed external force components px in the tangential direction and pn in the normal direction, and pm in the rotation as the beam displaces is  L W = (px u + pn w + pm ψ) dx 0

(19)

7 By means of Hamilton’s principle, the above energy expressions can be shown to be consistent with the equations of motion and the kinematical relations. This requires taking the variation of the energy functional (U -T -W ) between two arbitrary instants of time t1 and t2 . 

t2

δ

(U − T − W ) dt = 0

(20)

t1

Hamilton’s principle requires the coefficients of δu, δw, and δψ to vanish independently. The three equations of motion can be obtained by separating the terms multiplied by the variation of u, w, and ψ, respectively, as follows: 

A11 (u + κw) + B11 ψ  + A55 κ (w + ψ − κu) = 







 ¨ + (I2 + κI3 ) ψ¨ I1 + 2κI2 + κ2 I3 u



¨ −A11 κ (u + κw) − B11 κψ + A55 (w + ψ − κu) = I1 w 

¨ + I3 ψ¨ D11 ψ  + B11 (u + κw) − A55 (w + ψ − κu) = (I2 + κI3 ) u

(21) (22) (23)

For harmonic vibrations we assume u (x, t) = uf (x)cos ωt

(24a)

w (x, t) = wf (x)cos ωt

(24b)

ψ (x, t) = ψf (x)cos ωt

(24c)

where ω is the natural frequency; uf , wf , and ψf are the amplitudes of the harmonically varying tangential, radial, and sectional rotation motions, respectively. Introducing Eqs. (24a) to (24c) into Eqs. (21) to (23), the following equations of motion are obtained:    !  "  A11 uf + κwf + B11 ψf + A55 κ wf + ψf − κuf + ω 2 I1 + 2κI2 + κ2 I3 uf + (I2 + κI3 ) ψf = 0     −A11 κ uf + κwf − B11 κψf + A55 wf + ψf − κuf + I1 ω 2 wf = 0     D11 ψf + B11 uf + κwf − A55 wf + ψf − κuf + ω 2 {(I2 + κI3 ) uf + I3 ψf } = 0

(25) (26) (27)

3. An isogeometric formulation for laminated free-form curved beams This section shortly represents main features of NURBS as well as the application of IGA for the free vibration of a free-form curved beam. The NURBS curve and NURBS basis functions are constructed in the following. The details for NURBS-based geometry could be found in the NURBS book of Piegl and Tiller [29].

3.1. NURBS basis functions The NURBS is derived from B-splines which are piecewise polynomial curves composed of linear combinations of B-spline basis functions. The primary components of NURBS basis function is a knot vector H = (η1 , η2 , ..., ηn+p+1 ) which is a set of non-decreasing real numbers in the interval [η1 , ηn+p+1 ]. Here, ηi ∈ R is the ith knot, i is the knot index, i = 1, 2, ..., n+p+1 , p is the polynomial order and n is the number of basis functions. The intervals [η1 , ηn+p+1 ] and [ηi , ηi+1 ) are called a patch and a knot span, respectively [29]. A knot vector is called open and uniform if the first and the last knot have the multiplicity p + 1 and the other knots are in uniform space, respectively. Based on an

8 uniform open knot vector H, B-spline basis functions Bi,p (η) are defined recursively starting with piecewise constants (p = 0) as follows: ⎧ ⎨1

Bi,0 (η) =

if ηi ≤ η < ηi+1

⎩0

(28)

otherwise

For p = 1, 2, 3, ..., they are defined by Bi,p (η) =

η − ηi ηi+p+1 − η Bi,p−1 (η) + Bi+1,p−1 (η) ηi+p − ηi ηi+p+1 − ηi+1

(29)

Fig. 2 shows an example of the quadratic B-spline basis functions for the open and non-uniform knot vectors. If the internal knot has multiplicity k, the basis functions are C p−k -continuous at that knot. In addition, the basis functions are interpolatory at the ends of the interval and at the knot whose multiplicity is p (C 0 -continuous). A piecewise polynomial B-spline curves C(η) are defined by a linear combination of B-spline basis functions and coefficients over the parametric space. The coefficients are points in d-dimensional physical space Rd , referred to as control points. C(η) =

n 

Bi,p (η)Pi

(30)

i=1

The control points Pi ∈ Rd , i = 1, 2, ..., n construct the control polygon. The B-spline curve generated by the knot vector H = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5} and 8 control points is shown in Fig. 2. As mentioned before, “NURBS” is abbreviation for Non-Uniform Rational B-Spline. The term “Non-Uniform” refers to the knot vector which is not uniform. The term “Rational” refers to the basis functions which are piecewise rational polynomials. Rational B-Spline in Rd is the projection onto d-dimensional physical space of the polynomial B-spline defined in (d + 1)-dimensional homogeneous coordinate space. The rational polynomial curve is generated by the projective transformation of the B-spline curve. Homogeneous (d + 1)-dimensional control points are as follows: (Pjw )i = (Pj )i , j = 1, 2, ..., d w )i (Pd+1

= wi

(31) (32)

where wi are known as weights. The non-rational (d + 1)-dimensional B-spline curve is expressed from Eq. (30) as follows: Cw (η) =

n 

Bi,p (η)Pw i

(33)

i=1

Applying the projective transformation to Cw (η) yields the corresponding rational B-spline curve. n #

C(η) =

Bi,p (η) wi Pi

i=1 n #

i=1

= Bi,p (η) wi

n 

Ni,p (η) Pi

(34)

i=1

If the knot vector is non-uniform, NURBS basis functions for the NURBS curve are defined as follows: Bi,p (η) wi Ni,p (η) = # . n Bi,p (η) wi

(35)

i=1

The domain [η1 , ηn+p+1 ] is considered as a patch. Important properties of NURBS basis functions are as follows:

9 • Partition of unity, i.e.

n # i=1

Ni,p (η) = 1.

• Non-negativity, i.e. Ni,p (η)  0. • Local support, i.e. Ni,p (η) are non-zero only in the interval [ηi , ηi+p+1 ]. • Linear independence, i.e.

n # i=1

βi Ni,p (η) = 0 ⇔ βi = 0, i = 1, 2, ..., n.

3.2. NURBS-based finite element formulation The NURBS basis functions constructed become the finite element basis functions for the approximate displacement field. The degrees of freedom are the control variables and located at the control points. The number of finite elements are computed as the number of non-zero knot spans. In development of the NURBS-based finite element model, a typical element of a free-form curve beam which is ˆe ⊂ Ω ˆ and the curvilinear domain Ωe ⊂ Ω, is considered as shown in Fig. 3. The defined in the parametric domain Ω fundamental difference between the NURBS-based FEA and the conventional one is that in the new approach, the NURBS basis functions which are capable of exactly representing the known geometry are chosen to approximate the unknown solution fields. The arrangement of “NURBS coordinate” in Eq. (35) is inverted into the finite element notation for the case of one-dimensional parametric space. It should be noted that the beam is the one-dimensional element in the curvilinear and parametric spaces but its geometry is described in the two-dimensional physical rectangular space (ξ, ζ). The geometry a = (ξ, ζ) of the curved beam is defined based on NURBS functions and control points as shown in Eq. (34) and the solution field uf = (uf , wf , θf ) is described as follows: a = ae (η) =

nen 

Ni (η)Pei

i=1 nen 

uf ≈ uf e (η) =

Ni (η)uei

(36) (37)

i=1

where Pi = (ξi , ζi ) are the control points and the coefficients uei = (uei , wie , θie ) are denoted as control variables, nen is the number of control points for each element. The x-curvilinear coordinate, the Jacobian of transformation, the curvature are calculated, respectively, as follows: $

2  2 dζ dξ x (η) = + dη 0 dη 0 dη 0 η 0 =0 $   2 2 dζ dξ dx = + J (η) = dη dη dη 

κ (η) =

η 0 =η

dξ d2 ζ dη dη 2

−

dζ d2 ξ dη dη 2

J3

(38)

(39) (40)

In the following, substituting Eq. (37) into weak forms of Eqs. (25), (26) and (27), the finite element model of a typical element can be expressed as the standard eigenvalue problem. 

 K − ω2 M u = 0

(41)

10 where K is the element stiffness matrix given by ⎡

K11 K12 K13

⎤

⎢ ⎥ 21 22 23 ⎥ K=⎢ ⎣K K K ⎦ K31 K32 K33

(42)

where 11 Kij = 22 Kij = 33 Kij = 12 Kij 13 Kij

= =



Le+1

Le Le+1



Le Le+1



Le 21 Kji 31 Kji



 A11 Ni Nj + A55 κ2 Ni Nj dx



 A11 κ2 Ni Nj + A55 Ni Nj dx



 D11 Ni Nj + A55 Ni Nj dx



Le+1

=



 A11 κNi Nj − A55 κNi Nj dx



 −A55 κNi Nj + B11 Ni Nj dx



 A55 Ni Nj + B11 κNi Nj dx

Le



Le+1

= Le

23 32 Kij = Kji =



Le+1

Le

(43)

and M is the element mass matrix given by ⎡ ⎢ M=⎢ ⎣

M11

0

M13

0

M22

0

31

M

0

⎤ ⎥ ⎥ ⎦

(44)

33

M

where 11 Mij = 22 Mij



Le+1

Le



 I1 + 2κI2 + κ2 I3 Ni Nj dx

Le+1

=

33 Mij =



Le  Le+1 Le

13 31 Mij = Mji =

I1 Ni Nj dx I3 Ni Nj dx 

(45)

Le+1

Le

(I2 + κI3 ) Ni Nj dx

The assembly follows the standard finite element procedure. The standard Gauss-Legendre quadrature is used for numerical integration over the parametric domain. In the isogeometric analysis, each element in the parametric space is mapped into its corresponding element in the physical space as shown in Fig. 3. In order to transform from the parametric space to the parent element, an affine mapping is used [21]. The parametric coordinates of the quadrature ˆ e = [ηi , ηi+1 ] is transformed from its reference element coordinate η˜ ∈ Ω ˜ e as follows : ˆ e |Ω point η ∈ Ω η=

1 [(ηi+1 − ηi ) η˜ + (ηi+1 + ηi )] 2

Remarks on the present isogeometric model are:

(46)

11 (1) The x-curvilinear coordinate (arc-length) and curvature κ become functions of variable η in the parametric space, as shown in Eqs. (38) and (40), respectively, thus the free-form curved beam is completely modelled by NURBS. An analogue in the traditional FEA does not exist when the Lagrange interpolation functions are implemented to approximate the unknown solution field and a series of complex equations are used to model the arc-length and curvature. (2) The selective reduced integration technique is implemented in this study in consistency with the standard FEA. Unlike the standard C 0 Lagrange curved beam elements widely used in FEA, the shear and membrane locking still happen in the quadratic NURBS curved beam elements as stated in the work of Bouclier et al. [30]. However, in this study, the authors are interested in the use of higher orders of NURBS basis functions where the locking is not a problem (see [26; 31]).

3.3. NURBS-based refinement strategies It is well known in FEA that the accuracy in solutions increases by using refinement strategies. In NURBS-based isogeometric modelling, it is interesting that the basis could be enriched by the knot insertion and the order elevation while leaving the geometry unchanged. These two methods could be efficiently utilized in the isogeometric analysis. The insertion of new knot values corresponds to the h-refinement of the traditional FEA since the number of element is raised. The order elevation has common features with the classical p-refinement strategy in FEA because it increases the polynomial order of the basis. When the order of basis is elevated, the multiplicity of knot values increases correspondingly in order to preserve the continuity across elements. Subsequently if an unique knot is inserted in a non-zero knot span, the continuity across the new elements is reduced by one. This mechanism of mesh refinement in analysis is known as the k-refinement. Figs. 4 to 6 depict these three refinement strategies for the case of the elliptic curve. The initial geometry is generated by two non-zero knot spans of quadratic basis functions as shown in Fig. 4(a). The h-refinement technique is implemented in the initial curve by putting 2 or 6 knot values to create the refined curves in Figs. 4(b) or 4(c), respectively. It is realized that the h-refinement strategy creates NURBS elements with higher continuities. The two C 0 /C 0 -quadratic NURBS elements of initial curve result the four and eight quadratic elements with C 1 -continuity across the new interior element boundaries. For the p-refinement strategy, the order of basis for the initial quadratic elements is elevated to generate the cubic, quartic or quintic elements and the continuity of these new elements are preserved, as shown in Fig. 5(b-d). Finally, the k-refinement is utilized by h-refinement of the presented prefinement curves. The basis is enriched with higher continuity as well as polynomial order. Fig. 6 shows 4 elements of cubic, quartic, and quintic basis functions with C 2 -, C 3 - or C 4 -interior continuity at the new element boundaries, respectively, by using the k-refinement technique. It should be noted that two adjacent elements share 1, 2, 3, 4 or 5 control points if they are C 0 -, C 1 -, C 2 -, C 3 - or C 4 -continuous across their common boundary, respectively.

12 4. Results and discussion 4.1. Convergence study and verification In order to validate the accuracy and applicability of the isogeometric finite element analysis framework developed by this study, numerical results for the semi-elliptic curved beam with clamped-clamped boundary condition are presented and compared with existing results shown in literatures. The major and minor radii of the semi-elliptic curved beam are denoted as α and β, respectively. Therefore the length of span l is equal to 2α. The unrefined geometry of the elliptic curved beam has two quadratic NURBS elements created by five control points. The h-, p-, and k-refinement techniques are implemented for generating quadratic, cubic, quartic, and quintic NURBS elements on refined curves as described in details in Figs. 4 to 6. The [0◦ /90◦ ] cross-ply semi-elliptic curved beam is considered to examine the convergence properties of the present IGA model. The following elastic properties are chosen for the composite material utilized in the laminate stacking % sequences: E11 /E22 = 10, G12 /E22 = G13 /E22 = 0.6, G23 /E22 = 0.5, ν12 = 0.25 and E11 /ρ = 103 denotes the quantitative relation between the elastic modulus E11 and the volume mass density ρ. The slenderness ratio l/r = 200, where r is the radius of gyration of the cross-sectional area. The rates of convergence for the lowest ten non-dimensional frequencies of the curved beam with clamped-clamped end constraint and the ratio of minor radius to % major radius β/α = 0.4 are presented in Table 1. The non-dimensional frequencies are defined as ω ¯ i = ωi l2 ρA/E11 I where A is the cross-sectional area and I = bh3 /12 is the moment of inertia. The present results are compared with analytical solutions achieved by Tseng et al. [19] for the first six non-dimensional frequencies. It is found from Table 1 that the results from this study are in excellent agreement with the analytical solutions. Besides, it is a definite strength point that refinement strategies allow IGA to converge efficiently to the expected results. It can be realized from the first five lines of Table 1 that by using the h-refinement, the accuracy of the present IGA-based solution is increased as the number of elements is increased from 8 to 128. The more efficient technique, k-refinement, with higher order, higher continuity of NURBS gives fast convergence speed. With the aid of the k-refinement, the first ten non-dimensional frequencies converge with 32 quintic elements. It can be noted here that the 32 quintic NURBS elements are generated from 41 control points which correspond to 41 × 3 = 123 DOFs. In Table 2, the non-dimensional frequencies for the first six modes of the clamped-clamped semi-elliptic curved beams with [0◦ /90◦ ] cross-ply and [±45◦ ] angle-ply lay-ups by 32 quintic elements are given with the different orthotropy ratios: E11 /E22 = 10 and 30, and various values of β/α. It is seen from Table 2 that the present IGA results and those given by the dynamic stiffness analysis [19] are very consistent for all cases considered. Next, the laminated circular beam with simply supported boundary condition and with a rectangular cross-section having 1 m length, 0.025 m width, and 0.05 m height is considered. The material properties are: E11 = 138 GPa, E22 = 8.96 GPa, ν12 = 0.3, G12 = 7.1 GPa, G13 = G23 = 3.44 GPa and ρ = 1580 kg/m3 [10]. Table 3 shows % the lowest five non-dimensional frequencies ω ˜ i = ωi L2 ρA/E11 I (L is the arc-length of curved beam) with [45◦ 4 ] and [30◦ 2 /60◦ 2 ] lay-ups which have the bending-twisting and stretching-twisting couplings. The present results are compared with the exact solution by Hajianmaleki and Qatu [10] and results from 3D FEA [10] using the ANSYS finite element code. One can see from Table 3 that by using the present NURBS-based IGA model, accurate results are achieved with any kind of stacking sequences.

13 4.2. Parabolic curved beams The free vibration behavior of the laminated parabolic curved beam is described based on the present NURBSbased IGA framework. The schematic and NURBS-based geometrical modelling of the parabolic curved beam with the length of span l and the rise s are depicted in Fig. 7. Tables 4 and 5 tabulate the lowest ten non-dimensional frequencies ω ¯ with various ratios of beam rise to span length s/l for clamped-clamped and hinged-hinged boundary conditions, respectively. The results are obtained by using 32 quintic elements. The material properties are the same as the previous example described in the convergence test but it is assumed that the orthotropy ratio E11 /E22 = 15 and the slenderness ratio l/r = 30 for the results given in Tables 4 and 5. It is seen that the non-dimensional frequencies decrease significantly as the stacking sequence of the laminated parabolic curved beam changes from the [0◦ /90◦ ]s symmetric cross-ply to the [0◦ /90◦ ] cross-ply to the [±45◦ ] angle-ply for all ten modes listed. Besides, the non-dimensional frequencies of the hinged-hinged parabolic curved beam is smaller than those of the clamped-clamped parabolic curved beam, as expected. The variation of the non-dimensional frequencies for the first two modes of the [0◦ /90◦ ] cross-ply parabolic curved beam with clamped-clamped and hinged-hinged boundary conditions with respect to the rise to span length ratio s/l, the slenderness ratio l/r and the orthotropy ratio E11 /E22 is shown in Figs. 8, 9 and 10, respectively. From these figures, it is realized that when s/l, l/r or E11 /E22 varies, the lowest frequency can change from the symmetric mode to the anti-symmetric one and the crossover frequency, which is defined when the curved beam can vibrate in either one of two modes, is seen to occur. (1) The crossover is found for the clamped-clamped and hinged-hinged curved beams as the value of the s/l ratio is around 0.12 and 0.13, respectively in Fig. 8. (2) The same phenomenon occur only for the curved beam with the hinged-hinged end constraint and the corresponding value of the slenderness ratio is close to 18 in Fig. 9. (3) The crossover frequency is only seen for the clamped-clamped parabolic curved beam when E11 /E22 ≈ 3 in Fig. 10. The first ten mode shapes for the [0◦ /90◦ ] cross-ply laminated parabolic curved beam with the clamped-clamped boundary condition, s/l = 0.2, l/r = 30 and E11 /E22 = 15 are depicted in Fig. 11.

4.3. Elliptic ring In this example, the laminated elliptic ring is considered. The geometrical modelling of the elliptic ring is shown in Fig. 12 for the case of quintic NURBS with four non-zero knot spans. The lowest six non-dimensional frequencies with the [0◦ /90◦ ]s symmetric cross-ply, the [0◦ /90◦ ] cross-ply and [±45◦ ] angle-ply lay-ups are presented in Table 6 for various values of α/β when the slenderness ratio 2α/r = 30 and the orthotropy ratio E11 /E22 = 15. It is noted that there are three rigid body modes for the elliptic rings which are excluded in this study. For the perfect circular ring, the six frequencies are grouped into three pairs when the symmetric and anti-symmetric modes are identical. However, the frequency splitting between these pairs of modes is produced when the α/β ratio increases, as expected. The splitting is clear for the pair of the first and second modes and gradually reduces in case of higher-order modes. All the results presented in Table 6 are also produced by using 32 NURBS elements of the order five. Fig. 13 shows the variation of the first two non-dimensional frequencies for the [0◦ /90◦ ] cross-ply lamination of elliptic rings with different orthotropy ratios with respect to the slenderness ratio when α/β=1.6. From this figure, it can be realized that due to the effect of the shear deformation and the rotary inertia, the non-dimensional frequencies decrease as

14 the slenderness ratio reaches to the range of thick rings, as expected. Besides, when the orthotropy ratio increases, this range is extended. For example, for E11 /E22 = 2, the curved thick beam theory is needed to apply when 2α/r is smaller than 30. This value of 2α/r for E11 /E22 = 50 is 60. The first six mode shapes for the [0◦ /90◦ ] cross-ply elliptic ring with α/β = 1.6, 2α/r = 30 and E11 /E22 = 15 are presented in Fig. 14.

4.4. Free-form curved beam In our final example, a laminated curved beam with the Tschirnhausen’s cubic geometry, as a free-form curved beam, is considered to generate results of the free vibration analysis. The geometry of the Tschirnhausen’s cubic curved beam is defined as 27ay 2 = x2 (x + 9a) with x ∈ [−9a, 0] , y ≥ 0, a > 0, the length of span l = 9a [32]. It is well known in CAD that the free-form curve is modeled accurately by NURBS through the technique of interpolation. Therefore, the laminated free-form curved beam can be analysed for the natural frequencies and corresponding mode shapes by using the present NURB-based isogeometric approach. The NURBS-based geometry of the Tschirnhausen’s cubic curved beam is constructed by interpolating of 41 data points with the aid of the chord length method, as shown in Fig. 15(a). This NURBS curve of the order five with 45 control points and 40 knot spans is presented in Fig. 15(b). Using the present geometric modelling, the first six non-dimensional frequencies ω ¯ with [0◦ /90◦ ] cross-ply and [±45◦ ] angle-ply lay-ups are presented in Table 7. The first six mode shapes for the [0◦ /90◦ ] cross-ply Tschirnhausen’s cubic curved beam with and l/r = 30 and E11 /E22 = 15 are plotted in Fig. 16.

5. Concluding remarks The laminated free-form curved beam model considering all the effects of deepness, laminate coupling, axis extensibility, shear deformation, and rotary inertia was developed by using NURBS-based isogeometric approach in order to investigate the in-plane free vibration behavior of laminated curved beams and rings with variable curvature. The proposed model uses deep formulation along with the equivalent lamina modulus for calculation of ABD parameters. The convergence study of the present IGA has been carried out for the case of the semi-elliptic curved beam. It have been confirmed the efficiency and fast convergence speed of the k-refinement in IGA. The validity of the proposed framework has been confirmed through the comparison of the present results with the other solutions available in the literature. The non-dimensional natural frequencies and mode shapes of the laminate parabolic curved beam, elliptic ring, and Tschirnhausen’s cubic curved beam have been presented. The effects of laminated lay-ups, rise to span length, slenderness and orthotropy ratios on the natural frequencies of the laminated curved beams with variable curvature have been investigated. It was found that the crossover of the lowest frequency can occur when the rise to span length, slenderness or orthotropy ratio changes. These numerical results have been provided as benchmark solutions to supplement the available database in the published literature. By using the proposed NURBS-based IGA, the gap between the analysis of the curved beam with constant curvature and that with variable curvature has been eliminated, thus this method can be applied for investigating other behaviors of the free-form curved beam.

15 Acknowledgments

This research was supported by National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology through NRF2010-0019373 and 2012R1A2A1A01007405, and by Korea Ministry of Knowledge Economy under the National HRD support program for convergence information technology supervised by National IT Industry Promotion Agency through NIPA-2013- H0401-13-1003. The support is gratefully acknowledged.

References [1] Qatu MS. In-plane vibration of slightly curved laminated composite beams. J Sound Vib 1992;159(2):327-338. [2] Qatu MS, Elsharkawy AA. Vibration of laminated composite arches with deep curvature and arbitrary boundaries. Comput Struct 1993;47(2):305-311. [3] Qatu MS. Theories and analyses of thin and moderately thick laminated composite curved beams. Int J Solids Struct 1993;30(20):2743-2756. [4] Ahmed KM. Free vibration of curved sandwich beams by the method of finite elements. J Sound Vib 1971;18(1):61-74. [5] Ahmed KM. Dynamic analysis of sandwich beams. J Sound Vib 1972;21(3):263-276. [6] Raveendranath P, Singh G, Pradhan B. Application of coupled polynomial displacement fields to laminated beam elements. Comput Struct 2000;78(5):661-670. [7] Malekzadeh P. In-plane free vibration analysis of laminated thick circular deep arches. J Reinf Plast Compos 2007;26(18):1943-1951. [8] Malekzadeh P, Setoodeh AR, Barmshouri E. A hybrid layerwise and differential quadrature method for in-plane free vibration of laminated thick circular arches. J Sound Vib 2008(1-2);315:212-225. [9] Khdeir AA, Reddy JN. Free and forced vibration of cross-ply laminated composite shallow arches. Int J Solids Struct 1997;34(10):1217-1234. [10] Hajianmaleki M, Qatu MS. Static and vibration analyses of thick, generally laminated deep curved beams with different boundary conditions. Compos Part B: Eng 2012;43(4):1767-1775. [11] Oh SJ, Lee BK, Lee IW. Natural frequencies of non-circular arches with rotatory inertia and shear deformation. J Sound Vib 1999;219(1):23-33. [12] Oh SJ, Lee BK, Lee IW. Free vibration of non-circular arches with rotatory inertia and shear deformation. Int J Solids Struct 2000;37:4871-4891. [13] Huang CS, Tseng YP, Leissa AW, Nieh KY. An exact solution for in-plane vibrations of an arch having variable curvature and cross section. Int J Mech Sci 1998;40(11):1159-1173. [14] Tseng YP, Huang CS, Lin CJ. Dynamic stiffness analysis for in-plane vibrations of arches with variable curvature. J Sound Vib 1997;207(1):15-31. [15] Nieh KY, Huang CS, Tseng YP. An analytical solution for in-plane free vibration and stability of loaded elliptic arches. Comput Struct 2003;81(13):1311-1327. [16] Rossi RE, Laura PA, Verniere PL. In-plane vibrations of cantilvered non-circular arcs of non-uniform cross-section with a tip mass. J Sound Vib 1989;129(2):201-213. [17] Romanelli E, Laura PA. Fundamental frequencies of non-circular, elastic, hinged arcs. J Sound Vib 1972;24(1):17-22. [18] Wang TM, Moore JA. Lowest natural extensional frequency of clamped elliptic arcs. J Sound Vib 1973;30(1):1-7. [19] Tseng YP, Huang CS, Kao MS. In-plane vibration of laminated curved beams with variable curvature by dynamic stiffness analysis. Compos Struct 2000;50(2):103-114.

16 [20] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Engrg 2005;194(39-41):4135-4195. [21] Cottrell JA, Hughes TJR, Bazilevs Y, Isogeometric Analysis: Toward Integration of CAD and FEA. New York: John Wiley & Sons, 2009. [22] Reali A. An Isogeometric analysis approach for the study of structural vibrations. J Earthq Eng 2006;10(1):1-30. [23] Cottrell JA, Reali A, Bazilevs Y, Hughes TJR. Isogeometric analysis of structural vibrations. Comput Methods Appl Mech Engrg 2006;195(41-43):5257-5296. [24] Shojaee S, Valizadeh N, Izadpanah E, Bui T, Vu TV. Free vibration and bucking analysis of laminated composite plates using the NURBS-based isogeometric finite element method. Compos Struct 2012;94(5):1677-1693. [25] Thai CH, Nguyen-Xuan H, Nguyen-Thanh N, Le TH, Nguyen-Thoi T, Rabczuk T. Static, free vibration, and buckling analysis of laminated composite ReissnerMindlin plates using NURBS-based isogeometric approach. Int J Numer Meth Eng 2012;91(6):571-603. [26] Lee SJ, Park KS. Vibrations of Timoshenko beams with isogeometric approach. Appl Math Model 2013;37(22):9174-9190. [27] Nagy AP, Abdalla MM, G¨ urdal Z. Isogeometric design of elastic arches for maximum fundamental frequency. Struct Multidisc Optim 2011;43(1):135-149. [28] Vinson JR, Sierakowski RL. The Behavior of Structures Composed of Composite Materials. Netherlands: Kluwer Academic Publishers, 1986. [29] Piegl L, Tiller W. The NURBS book, 2nd Edition. New York: Pringer-Verlag, 1997. [30] Bouclier R, ElGuedj T, Combescure A. Locking free isogeometric formulations of curved thick beams. Comput Methods Appl Mech Engrg 2012;245-246:144-162. [31] Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. Isogeometric shell analysis: The Reissner-Mindlin shell. Comput Methods Appl Mech Engrg 2010;199:276-289. [32] Lawrence JD. A Catalog of Special Plane Curves. New York: Dover, 1972.

17 TABLE I Convergence rate of non-dimensional frequencies ω ¯ i = ωi l 2

% ρA/E11 I for the clamped-clamped [0◦ /90◦ ] cross-ply

semi-elliptic curved beam with NURBS-based refinement schemes when β/α = 0.4 and E11 /E22 = 10. p

nel a

Mode 1

2

3

4

5

2

3

4

5

6

7

8

9

10

8

24.2612

34.0980

86.3093

141.7366

246.458

348.315

589.975

706.564

968.850

1018.836

16

22.5021

30.0796

69.9474

89.3347

150.832

185.972

281.285

319.978

453.363

476.727

32

22.4223

29.9103

68.5922

87.0118

141.074

171.163

240.765

274.430

365.678

376.232

64

22.4209

29.9069

68.5655

86.9583

140.857

170.794

239.663

273.035

361.680

373.365

128

22.4208

29.9068

68.5648

86.9567

140.851

170.783

239.633

272.995

361.575

373.285

8

22.8937

30.7788

75.5159

100.4230

194.717

237.125

379.578

432.832

707.531

928.524

16

22.4277

29.9220

68.7185

87.3165

142.597

174.577

251.580

289.954

406.368

414.493

32

22.4209

29.9069

68.5660

86.9592

140.863

170.805

239.712

273.107

361.966

373.589

64

22.4208

29.9068

68.5648

86.9567

140.851

170.783

239.632

272.994

361.573

373.284

128

22.4208

29.9068

68.5648

86.9567

140.851

170.783

239.631

272.993

361.570

373.282

8

22.4726

30.0818

70.3021

90.6123

152.695

199.345

338.509

366.477

561.025

573.499

16

22.4213

29.9079

68.5782

86.9979

141.151

171.498

243.235

278.064

381.609

388.406

32

22.4208

29.9068

68.5648

86.9567

140.851

170.784

239.635

273.000

361.599

373.307

64

22.4208

29.9068

68.5648

86.9567

140.851

170.783

239.631

272.993

361.570

373.282

128

22.4208

29.9068

68.5648

86.9567

140.851

170.783

239.631

272.993

361.570

373.282

8

22.4303

29.9362

68.8925

87.4927

146.895

183.847

263.686

322.192

431.815

568.475

16

22.4208

29.9068

68.5663

86.9614

140.912

170.937

240.539

274.870

369.275

379.844

32

22.4208

29.9068

68.5648

86.9567

140.851

170.783

239.631

272.993

361.570

373.282

22.4208

29.9068

68.5648

86.9567

140.851

170.783

239.631

272.993

361.570

373.282

22.421

29.907

68.568

86.959

140.86

170.789

64 Ref. [19] a Number

of elements used for IGA

18 TABLE II Non-dimensional frequencies ω ¯ i = ωi l 2

% ρA/E11 I for the clamped-clamped semi-elliptic curved beam generated by

32 quintic NURBS elements. Lay-up

E11 /E22

β/α

Mode 1

◦

◦

[0 /90 ]

10

30

[±45◦ ]

30

86.957

5

6

140.85

170.78

Ref. [19]

22.421

29.907

68.568

86.959

140.86

170.79

0.6

This study

16.492

26.874

54.476

78.501

118.30

153.99

Ref. [19]

16.492

26.875

54.479

78.505

118.31

154.00

This study

12.182

23.592

45.007

67.964

99.351

131.79

Ref. [19]

12.182

23.593

45.010

67.968

99.357

131.80

This study

18.047

23.903

54.517

68.992

110.51

134.97

Ref. [19]

18.047

23.903

54.519

68.994

110.51

134.98

This study

13.271

21.531

43.451

62.437

93.453

121.56

Ref. [19]

13.271

21.531

43.453

62.440

93.459

121.57

9.8037

18.935

35.997

54.196

78.832

104.27 104.28

0.6

68.565

4

This study

0.4

29.907

3

0.4

0.8

22.421

2

0.8

This study Ref. [19]

9.8041

18.935

35.999

54.198

78.836

0.4

This study

9.4454

12.504

28.952

36.081

59.655

69.352

Ref. [19]

9.4456

12.504

28.954

36.082

59.659

69.355

This study

6.9498

11.307

23.021

32.981

50.130

64.089

Ref. [19]

6.9501

11.307

23.022

32.982

50.133

64.092

This study

5.1428

9.9417

19.013

28.652

42.074

55.465

Ref. [19]

5.1430

9.9421

19.014

28.654

42.076

55.468

0.6

0.8

19 TABLE III Non-dimensional frequencies ω ˜ i = ω i L2

%

ρA/E11 I for the simply supported circular curved beam generated by 16

quintic NURBS elements. Lay-up ◦

[45 ]4

[30◦ 2 /60◦ 2 ]

κL = 1

Mode

κL = 2

This study

Exact solution [10]

3D FEM [10]

This study

Exact solution [10]

3D FEM [10]

1

2.8301

2.830

2.871

1.65949

1.659

1.682

2

12.483

12.48

12.51

11.120

11.12

11.17

3

27.841

27.84

27.82

26.484

26.48

26.52

4

47.971

47.97

47.72

46.667

46.67

46.48

5

71.904

71.90

71.16

70.667

70.67

70.00

1

2.9514

2.951

3.002

1.7381

1.738

1.765

2

13.004

13.00

13.07

11.670

11.67

11.75

3

28.928

28.93

28.96

27.733

27.73

27.83

4

49.688

49.69

49.50

48.697

48.70

48.52

5

74.232

74.23

73.53

73.457

73.46

72.74

20 TABLE IV Non-dimensional frequencies ω ¯ i = ωi l 2

%

ρA/E11 I for clamped-clamped parabolic laminated curved beams (l/r = 30

and E11 /E22 = 15). s/l = 0.2

Mode

s/l = 0.3

s/l = 0.4

[0◦ /90◦ ]s

[0◦ /90◦ ]

[±45◦ ]

[0◦ /90◦ ]s

[0◦ /90◦ ]

[±45◦ ]

[0◦ /90◦ ]s

[0◦ /90◦ ]

[±45◦ ]

1

21.6413

17.2559

11.8658

17.9134

13.9881

9.88571

14.7627

11.1646

7.74356

2

27.1837

25.3212

12.4792

30.7576

26.2063

14.1695

27.7674

22.5114

15.1213

3

38.5170

33.4206

24.3329

35.6189

34.5662

20.4878

36.1985

34.4173

16.9296

4

52.3080

45.4861

28.2473

46.5890

40.0719

26.3600

40.7361

36.5283

24.0543

5

66.6394

59.9970

37.1108

61.0900

54.1081

31.6964

53.8305

47.0047

26.4072

6

68.3572

61.1156

50.3288

61.5902

55.3851

43.7696

55.6483

50.1966

37.1602

7

83.8598

76.5150

55.0703

75.3235

68.0975

51.2519

66.8320

59.5825

48.0168

8

98.8660

92.0394

64.9821

89.0755

82.1973

56.7541

79.2173

72.1018

48.6140

9

114.096

100.064

79.9128

103.040

94.8080

70.1985

91.9241

85.0655

60.5729

10

126.978

107.658

81.1523

116.608

96.4856

74.3427

104.503

91.4911

67.4113

21 TABLE V Non-dimensional frequencies ω ¯ i = ωi l 2

%

ρA/E11 I for hinged-hinged parabolic laminated curved beams (l/r = 30

and E11 /E22 = 15). s/l = 0.2

Mode

s/l = 0.3

s/l = 0.4

[0◦ /90◦ ]s

[0◦ /90◦ ]

[±45◦ ]

[0◦ /90◦ ]s

[0◦ /90◦ ]

[±45◦ ]

[0◦ /90◦ ]s

[0◦ /90◦ ]

[±45◦ ]

1

18.6686

13.1388

8.47157

14.8223

10.4170

6.50481

11.6921

8.10560

4.95711

2

26.4386

20.4995

11.4371

30.2290

21.9298

14.1504

25.6860

18.0535

12.7640

3

36.2992

27.6468

19.4705

33.4764

28.6611

15.9959

35.6720

29.6074

15.3845

4

51.5991

41.7188

28.2461

45.6463

35.7072

26.3407

39.5262

31.4725

22.3305

5

66.5268

50.8060

32.4902

60.4726

47.4337

27.3166

53.0903

43.9380

24.0610

6

67.8499

58.7063

46.4224

61.5190

51.4845

39.8066

55.6286

44.1419

33.3451

7

83.5797

75.1446

54.9228

74.9771

66.5041

51.2346

66.4112

57.5782

45.1734

8

98.6890

90.0902

61.5565

88.8567

80.6326

53.2767

78.9471

70.7315

48.0106

9

114.002

96.8193

77.1360

102.915

91.4110

67.2803

91.7617

84.1955

57.5910

10

126.744

106.953

81.1523

116.440

95.6637

74.3373

104.383

88.0149

67.3191

22 TABLE VI Non-dimensional frequencies ω ¯ i = ωi (2α)2

%

ρA/E11 I for elliptical rings with various α/β ratios (2α/r = 30 and

E11 /E22 = 15). α/β 1.0

Lay-up 2

3

4

5

6

8.14860

8.14860

19.1996

19.1996

30.8256

30.8256

5.18162

5.18162

13.5300

13.5300

23.5062

23.5062

[±45 ]

3.27953

3.27953

8.90911

8.90911

16.2323

16.2323

[0◦ /90◦ ]s

9.30353

9.42710

21.5404

21.5454

34.2885

34.2907

6.06015

6.15518

15.6207

15.6214

26.8872

26.8900

[±45 ]

3.84563

3.89586

10.3696

10.3739

18.8061

18.8072

[0◦ /90◦ ]s

10.1379

10.5892

23.1331

23.1693

36.7636

36.7873

6.71648

7.07082

17.1155

17.1263

29.3613

29.3915

[±45 ]

4.28504

4.47114

11.4499

11.4836

20.7791

20.7906

[0◦ /90◦ ]s

10.7424

11.6597

24.1878

24.2955

38.5281

38.6069

7.20654

7.93610

18.1582

18.2032

31.1494

31.2502

[±45 ]

4.62816

5.00982

12.2221

12.3295

22.2703

22.3157

[0◦ /90◦ ]s

11.1885

12.6571

24.8659

25.0875

39.7936

39.9592

7.57973

8.75748

18.8705

18.9857

32.4458

32.6585

4.90175

5.51629

12.7556

12.9884

23.3657

23.4971

◦

◦

◦

◦

[0 /90 ]s [0 /90 ] ◦

1.2

◦

◦

[0 /90 ] ◦

1.4

◦

◦

[0 /90 ] ◦

1.6

◦

◦

[0 /90 ] ◦

1.8

◦

◦

[0 /90 ] ◦

[±45 ] 2.0

Mode 1

◦

◦

[0 /90 ]s

11.5256

13.5950

25.2807

25.6551

40.7121

40.9866

[0◦ /90◦ ]

7.87227

9.53972

19.3413

19.5689

33.3986

33.7534

5.12587

5.99424

13.1068

13.5163

24.4222

24.5693

◦

[±45 ]

23 TABLE VII Non-dimensional frequencies ω ¯ i = ωi l 2

%

ρA/E11 I for the Tschirnhausen’s cubic laminated curved beams (l/r = 30

and E11 /E22 = 15). Mode

clamped-clamped [0◦ /90◦ ]

hinged-hinged [±45◦ ]

[0◦ /90◦ ]

[±45◦ ]

l/r=30

l/r=100

l/r=30

l/r=100

l/r=30

l/r=100

l/r=30

l/r=100

1

12.8174

16.6931

8.76372

10.9375

9.13513

10.1497

6.15229

6.59864

2

22.7220

33.3211

13.8278

21.7906

21.8297

26.2754

12.8005

16.7436

3

34.8947

59.0957

20.4581

38.8770

26.8522

48.7339

16.8203

32.1066

4

38.8451

85.3094

28.2420

47.7212

38.1375

77.9197

26.1540

46.3679

5

52.7510

112.693

37.8968

63.0989

50.2795

98.8063

35.7344

53.8060

6

66.8694

124.891

45.0675

84.1978

64.6114

117.363

41.7080

74.7010

24

w

\

x z h/2

0

]

u

1 N ( x)

h/2

zk

0

[

zk-1

Q N

N M Q

FIG. 1 Defining sketch for curve beam with variable curvature

M

x=L L

25

.P

3

.

P2

.P

8

.

P5 P1

. P4

.

.P

.

P6

7

(a) B-spline curve 1

N1,2

N5,2

N3,2

N2,2

N8,2

N4,2

N6,2

N7,2

0.5

0 0,0,0

C0

1

2

3,3

4

5,5,5

C1

C1

C0

C1

C0

(b) Quadratic basis functions FIG. 2 B-spline curve and its corresponding quadratic basis functions with open, non-uniform knot vector H = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}

26

] [

Curvilinear domain

e

x x-11

ˆe 

Ki

P Parametric t i domain d i

Ki 1

K

f

e  -1

Parent element

1

K FIG. 3 Geometrical and affine mappings for integration by Gaussian quadrature on NURBS elements of curved beam

27

.

.

.

.

.

e1

1

e2

0 000 0,0,0

1 1 , 2 2

111 1,1,1

(a) Initial quadratic curve with H = {0, 0, 0, 12 , 12 , 1, 1, 1}

. . .

. .

e1

1

. .

0 000 0,0,0

e3

e2

e4

1 1 , 2 2

1 4

3 4

111 1,1,1

(b) h-refined quadratic curve with H = {0, 0, 0, 14 , 12 , 12 , 34 , 1, 1, 1}

..

.

. ... .

.

1

..

0 000 0,0,0

e2 e3

e1

1 8

1 4

e4

3 8

1 1 , 2 2

e5

e6

5 8

e8

e7

3 4

7 8

111 1,1,1

(c) h-refined quadratic curve with H = {0, 0, 0, 18 , 14 , 38 , 12 , 12 , 58 , 34 , 78 , 1, 1, 1} FIG. 4 Semi-elliptic curve modeled exactly by NURBS: unrefined and h-refined curves with their corresponding basis

28

.

.

.

. . .

. .

.

1

e1

0 000 0,0,0

e2

1 1 , 2 2

111 1,1,1

(a) Initial quadratic curve with H = {0, 0, 0, 12 , 12 , 1, 1, 1}

.

. .

1

e1

0 00,0,0,0 000

e2

1 1 1 , , 2 2 2

1111 1,1,1,1

(b) p-refined cubic curve with H = {0, 0, 0, 0, 12 , 21 , 12 , 1, 1, 1, 1}

. .

. . . . .

1

. .

e1

0 00,0,0,0,0 0000

e2

1 1 1 1 , , , 2 2 2 2

11111 1,1,1,1,1

(c) p-refined quartic curve with H = {0, 0, 0, 0, 0, 12 , 12 , 21 , 12 , 1, 1, 1, 1, 1}

. .

. . . . . . .

1

. .

0 00,0,0,0,0,0 00000

e1

e2

1 1 1 1 1 , , , , 2 2 2 2 2

111111 1,1,1,1,1,1

(d) p-refined quintic curve with H = {0, 0, 0, 0, 0, 0, 12 , 12 , 12 , 12 , 12 , 1, 1, 1, 1, 1, 1} FIG. 5 Semi-elliptic curve modeled exactly by NURBS: unrefined and p-refined curves with their corresponding basis

29

.

.

.

.

.

e1

1

e2

0 000 0,0,0

1 1 , 2 2

111 1,1,1

(a) Initial quadratic curve with H = {0, 0, 0, 12 , 12 , 1, 1, 1}

.

..

. . .

.

1

..

e1

0 00,0,0,0 000

e3

e2

1 1 1 , , 2 2 2

1 4

e4

3 4

1111 1,1,1,1

(b) k-refined cubic curve with H = {0, 0, 0, 0, 14 , 12 , 12 , 12 , 12 , 34 , 1, 1, 1, 1}

..

.

. ... .

.

1

..

e1

0 00,0,0,0,0 0000

e3

e2

1 1 1 1 , , , 2 2 2 2

1 4

e4

3 4

11111 1,1,1,1,1

(c) k-refined quartic curve with H = {0, 0, 0, 0, 0, 14 , 12 , 12 , 12 , 12 , 34 , 1, 1, 1, 1, 1}

. ..

. . ... . .

. ..

1

0 00,0,0,0,0,0 00000

e1

e2

1 4

1 1 1 1 1 , , , , 2 2 2 2 2

e3

e4

3 4

111111 1,1,1,1,1,1

(d) k-refined quintic curve with H = {0, 0, 0, 0, 0, 0, 14 , 12 , 12 , 12 , 12 , 12 , 34 , 1, 1, 1, 1, 1, 1} FIG. 6 Semi-elliptic curve modeled exactly by NURBS: unrefined and k-refined curves with their corresponding basis

30

s l/2

l

(a) Schematic of the parabolic curved beam

.

. .

. .

. (b) Quintic NURBS parabola

FIG. 7 Schematic and NURBS-based geometrical modelling of the parabolic curved beam

31

30 1stsymmetricmode(cc) 1st symmetric mode (c c) 1stantisymmetricmode(cc) 1stsymmetricmode(hh) 1stantisymmetricmode(hh)

Nond dimensionalfrequency,Z

25

20

15

10

5

0 0

0.2

0.4

0.6

0.8

Rise to span length ratio s/l Risetospanlengthratio,s/l FIG. 8 Variation of the non-dimensional frequencies for the [0◦ /90◦ ] cross-ply laminated parabolic curved beams with respect to the rise to span length ratio s/l when l/r = 30 and E11 /E22 = 15.

32

45

Nond dimensionalfrequency,Z

40 35 30 25 20 15 1stsymmetricmode(cc) 1stantisymmetricmode(cc) 1stsymmetricmode(hh) 1stantisymmetricmode(hh)

10 5 0 0

20

40

60

80

100

Slenderness ratio l/r Slendernessratio,l/r FIG. 9 Variation of the non-dimensional frequencies for the [0◦ /90◦ ] cross-ply laminated parabolic curved beams with respect to the slenderness ratio l/r when s/l = 0.2 and E11 /E22 = 15.

33

40 1stsymmetricmode(cc) 1stantisymmetricmode(cc) 1stsymmetricmode(hh) 1stantisymmetricmode(hh)

Nond dimensionalfrequency,Z

35

30

25

20

15

10

5 0

10

20

30

40

50

Orthotropy ratio,E ratio E11/E22 FIG. 10 Variation of the non-dimensional frequencies for the [0◦ /90◦ ] cross-ply laminated parabolic curved beams with respect to the orthotropy ratio E11 /E22 when s/l = 0.2 and l/r = 30.

34

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

FIG. 11 The first ten mode shapes of the [0◦ /90◦ ] cross-ply parabolic curved beam with clamped-clamped end constraint (s/l = 0.2, l/r = 30 and E11 /E22 = 15). Solid and dashed lines represent deformed and undeformed configurations, respectively.

35

. . .

. . . . . . . .. . . . . .

. . .

FIG. 12 Quintic NURBS ellipse with the knot vector H = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4} and 21 control points of 4 C 0 /C 0 NURBS elements

36

15

Non ndimensionaalfrequency,,Z

13

11

9

7

5

( 11/E / 22=2) 2) 1stsymmetricmode d (E 1stantisymmetricmode (E11/E22=2) 1stsymmetricmode (E11/E22=15) 1stantisymmetricmode (E11/E22=15) 1stsymmetricmode (E11/E22=50) 1stantisymmetricmode (E11/E22=50)

3

1 0

20

40

60

80

100

Slenderness ratio, 2D/r Slendernessratio,2 FIG. 13 Variation of the non-dimensional frequencies for the [0◦ /90◦ ] cross-ply laminated elliptic ring with respect to the slenderness ratio l/r with different orthotropy ratios when α/β = 1.6.

37

(a)

(b)

(c)

(d)

(e)

(f)

FIG. 14 The first six mode shapes of the [0◦ /90◦ ] cross-ply elliptic ring (α/β = 1.6, 2α/r = 30 and E11 /E22 = 15)

38

                       (a) Quintic NURBS curve constructed through 41 interpolation points

. . .......... . . . . . . . . ..... . . . . .. . . .. .. . (b) Constructed quintic NURBS curve with 45 control points FIG. 15 Quintic NURBS curve of 40 elements is constructed by interpolation method of 41 data points given by the Tschirnhausen’s cubic curve which is defined as 27ay 2 = x2 (x + 9a) with x ∈ [−9a, 0] , y ≥ 0, a > 0

39

(a)

(b)

(c)

(d)

(e)

(f)

FIG. 16 The first six mode shapes of the [0◦ /90◦ ] cross-ply Tschirnhausen’s cubic curved beam with clamped-clamped boundary condition (l/r = 30 and E11 /E22 = 15)