Generalized warping and distortional analysis of curved beams with isogeometric methods

Computers and Structures 191 (2017) 33–50

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Generalized warping and distortional analysis of curved beams with isogeometric methods Ioannis N. Tsiptsis ⇑, Evangelos J. Sapountzakis Institute of Structural Analysis and Antiseismic Research, School of Civil Engineering, National Technical University of Athens, Zografou Campus, GR-157 80 Athens, Greece

a r t i c l e

i n f o

Article history: Received 8 February 2017 Accepted 12 June 2017

Keywords: Curved beam Generalized warping Distortion Isogeometric analysis Curvature

a b s t r a c t Towards improving conventional beam elements in order to include distortional effects in their analysis, independent parameters have been taken into account in this study. Curved beam’s behavior becomes more complex, even for dead loading, due to the coupling between axial force, bending moments and torque that curvature produces. Thus, the importance of simulating geometry exactly arises in order to approximate accurately the response of the curved beam. For this purpose, the isogeometric tools (bsplines and NURBS), either integrated in the Finite Element Method (FEM) or in a Boundary Element based Method (BEM) called Analog Equation Method (AEM), are employed in this contribution for the static analysis of horizontally curved beams of open or closed (box-shaped) cross sections. Responses of the stress resultants, stresses and displacements to static loading have been studied for various cross sections. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Refined models either straight or curved with shell or solid elements are widely used in structures, such as for example the deck of a bridge with a thin-walled cross section, for stress or strain analysis. Over time, higher order beam theories have been developed to include nonuniform warping and distortional effects in beam elements which exhibit important advantages over more refined approaches [1]. The evaluation of the cross sectional properties, which are finally incorporated in the one-dimensional beam analysis, is associated with the accuracy of the model regarding the cross sectional behavior. Over the past decades, classical beam theories based on specific assumptions fail to describe accurately the structural behavior of beam elements, especially in more complex formulations such as in curved geometries. Among these theories, that of Saint-Venant (SV) still plays a crucial role due to the fact that the analysis reduces to the evaluation of warping and distortional functions over the cross sectional domain. However, this solution is exact for the uniform warping of a beam (warping/distortional deformations are not restrained). Towards improving SV theory, several researchers investigated the so-called SV’s principle (stated in [2]) as well as the SV’s end-effects in order to derive a

⇑ Corresponding author. E-mail addresses: [email protected] (I.N. Tsiptsis), [email protected] (E.J. Sapountzakis). URL: http://users.ntua.gr/cvsapoun/ (E.J. Sapountzakis). http://dx.doi.org/10.1016/j.compstruc.2017.06.007 0045-7949/Ó 2017 Elsevier Ltd. All rights reserved.

more general formulation of beams’ kinematics. In most of these studies, the solution is obtained as the sum of the SV’S solution and the residual displacements corresponding to the end-effects, as it will be later explained. In the majority of past research works, thin-walled cross sections have been studied due to their low self-weight comparing to solid ones and, thus, their use in practice. These cross sections are more susceptible to torsional and distortional effects. Vlasov (1961) in [3] presented the Thin Tube Theory (TTT) and treated different cross section types as special cases of this general theory. Kollbrunner and Basler (1969) in [4] and Heilig (1971) in [5] were later reformulated TTT for multi-cell boxes with arbitrary cross sections. Kristek (1970) in [6] obtained analytical solution for simple practical cases and separated the analysis of transverse distortion from that of torsion with longitudinal warping employing the superposition principle. Wright, Abdel-Samad and Robinson (1968) in [7] studied the distortional warping response of singlecell box girders with longitudinally and transversely stiffened plates employing the beam on elastic foundation (BEF) analogy. Steinle (1970) in [8] tackled the torsional distortion problem and introduced distortional stress resultants in the analysis. Kollbrunner and Hajdin (1975) in [9] dealt with the extension of the beam theory of prismatic folded structures to include the deformation of the cross section for open and closed cross sections including warping. Other research efforts later expanded TTT considering only box-shaped cross sections (single- or multi-cell) and, thus, being not general [10–15]. Schardt (1989, 1994) in [16,17]

34

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

developed an advanced formulation known as Generalized Beam Theory (GBT) which is a generalization of the classical Vlasov beam theory in order to incorporate flexural and torsional distortion. A distinguishing feature of GBT stems from the general character of its cross sectional analysis which enables the determination of cross-section deformation modes as well as their categorization to global, distortional or local ones. Discretization within the frame of GBT cross sectional analyses depends on the topology of nodes (dependent, independent or intermediate nodes for branched or unbranched sections) and cross sectional shapes (open- or closed-shaped). Further developments of GBT avoid some of its cumbersome procedures through eigenvalue cross sectional analysis [18–23]. The classification of the cross sections with respect to their geometry is not needed with this approach and the most important shape modes are obtained (compared to those derived by GBT). These approaches are employed nowadays by several researchers. Camotim, Silvestre and co-researchers expanded the method to cover a variety of cross sections, orthotropic materials, as well as geometrically nonlinear and inelastic problems [24– 26]. Towards solving the problem for arbitrarily shaped homogeneous or composite cross sections, El Fatmi and Ghazouani (2011) in [27] presented a higher order composite beam theory (denoted HOCBT) that starts from the exact expression of SV’s solution and introduces in- and out-of-plane independent warping parameters for symmetric orthotropic cross sections with the ability to extended it for arbitrary ones. However, in-plane warpings are only due to the flexural and axial deformation modes and, thus, it could be stated that this research effort studies Poisson ratio effects rather than distortional ones. Ferradi and Cespedes (2014) in [28] presented the formulation of a 3D beam element solving an eigenvalue problem for the distortional behavior of the cross section (in-plane problem) and computing warping functions separately by using an iterative equilibrium scheme. Genoese, Genoese, Bilotta and Garcea (2014) in [29] developed a beam model with arbitrary cross section taking into account warping and distortion with their evaluation being based on the solution of the 3D elasticity problem for bodies loaded only on the terminal bases and a semi-analytic finite element formulation. Finally, Dikaros and Sapountzakis (2016) in [30] presented a general boundary element formulation for the analysis of composite beams of arbitrary cross section taking into account the influence of generalized cross sectional warping and distortion due to both flexure and torsion. In this proposal, distortional and warping functions are evaluated by the same eigenvalue problem and in order of importance. Regarding horizontally curved beams subjected to vertical or radial loads, they inherently exhibit a more complex behavior comparing to straight formulations due to the fact that the effects of primary and secondary torsion are always coupled to those of bending and cross section distortion either for centered or eccentric loads. Dabrowski (1968) in [31] elaborated Vlasov’s theory and introduced distortional behavior of box girders with a symmetric cross section. His model introduces the distortion angle as the single degree of freedom which measures the magnitude of the cross-sectional distortion. Bazant and Nimeiri (1974) in [32] proposed the skew-ended finite element in order to implement the theory of non-uniform torsion for straight or curved thinwalled cross sections. Oleinik and Heins (1975) in [33], and Heins and Oleinik (1976) in [34] employing Vlasov’s and Dabrowski’s theories studied the structural behavior of curved box girders. Inplane deformations were approximated using a differential equation which was solved employing the finite difference method. In addition to these research efforts, Sakai and Nagai (1981) in [35], and Nakai and Murayama (1981) in [36] presented several results on the design procedures of the intermediate diaphragms for curved girders and noted that these play a very important role in moderating distortional warping of girders. Meanwhile, Martin

and Heins (1978) in [37] expanded Dabrowski’s equation, which predicts the cross-sectional deformations, so that the angular deformations induced at given points along an I-girder curved bridge can be calculated. Zhang and Lyons (1984) in [38,39] employed Dabrowski’s theory combined with Finite element method to develop a multi-cell box element for the analysis of curved bridges. Nakai and Yoo (1988) in [40] presented an extended study on the analysis and design of curved steel bridges. Yabuki and Arizumi (1989) in [41] employing BEF analogy for distortion proposed spacing provisions which can be utilized for steelplated intermediate diaphragms. Razaqpur and Li (1994) in [42] extended their previous theory to curved thin-walled box beams. Petrov and Geradin (1998) in [43] employing the same concept with El Fatmi and Ghazouani [27] for straight beams formulated a theory for curved and pre-twisted beams of arbitrary homogeneous cross sections, covering geometrically nonlinear range as well. Kim and Kim (2002) in [44] developed a theory for thinwalled curved beams of rectangular cross section by extending the theory developed earlier for straight beams taking into account warping and distortional deformations. Park, Choi and Kang (2005) in [45] expanded their previous work [13], which was limited to straight box girder bridges, to curved formulations. They developed a curved box beam element which was employed in order to develop design charts for adequate spacing of the intermediate diaphragms of curved bridges. Flexural and torsional displacement functions have been based on those proposed for doubly symmetric cross section by Kang and Yoo (1994) in [46] while distortional functions have been derived for a mono-symmetric cross section as in [12]. Despite the practical interest of their study, their proposal cannot accommodate elastic constraints and due to other assumptions made lacks of generality. In other research efforts, the vibration problems of thin-walled curved box girder bridges due to moving loads have been investigated. The curved box girder bridges have been numerically modelled using finite elements which take into account the torsional warping, distortion and distortional warping [47–49]. Other recent research efforts as the following ones mainly constitute design guides with new formulae for specific practical cases rather than a generalized theory for the analysis of curved beams. Particularly, in the study of Zhang, Hou, Li and Wang (2015) [50], a curved girder is simplified to straight one by using the M/r method and calculation formulae for determining the required diaphragm spacing are obtained by regression analyses. Yoo, Kang and Kim (2015) in [51] applied the concept of the BEF analogy for the analysis of distortional stresses of horizontally curved box-girders. The proposed procedure is capable of handling simple or continuous single cell box girders (or separated multi-cell box girders) with rigid or deformable interior diaphragms or cross-frames. Towards establishing a more general theory, Arici and Granata (2016) in [52] employed the Hamiltonian Structural Analysis Method for the analysis of straight and curved thin-walled structures on elastic foundation extending the so-called GBT. To the authors’ knowledge, there are no research efforts that introduce a unified distortional and warping eigenvalue analysis of arbitrarily shaped cross sections to the analysis of curved beams. In modern regulations and design specifications, the importance of torsional and distortional effects in stress or strain analysis of structural members is recognized. Particularly, in sub-sections 6.2.7.1 and 6.2.7.2 of EN 1993-2, Eurocode 3: Design of steel structures – Part 2: Steel bridges, regarding torsion, the designer is obliged to keep the distortional stresses under a specific limiting value or follow some general design rules in case of neglecting distortion. These are presented in clauses (1)–(9) of section 6.2.7, regarding torsion, of EN 1993-1-1, Eurocode 3: Design of steel structures – Part 1-1: General rules and rules for buildings. Nevertheless, no guidelines and specific modelling methodologies

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

offered for the aforementioned effects. It should also be noted that most of the provisions of Eurocode 3 regarding torsion are valid only when distortional deformation can be neglected. The same case is when the stability of uniform members is checked as it is mentioned in clause (1) of section 6.3.3 of EN 1993-1-1. Distortional buckling is encountered in EN 1993-1-3, Eurocode 3 – Design of steel structures – Part 1-3: General rules – Supplementary rules for cold-formed member and sheeting, and EN 1993-1-4, Eurocode 3 – Design of steel structures – Part 1-4: General rules – Supplementary rules for stainless steels. It is taken into account through reduction factors or special arrangements in order to prevent distortion. In addition to these, distortional effect is also suggested be taken into account during the design of unreinforced joints (Section 7.5.2.1(7) of EN 1993-1-8, Eurocode 3: Design of steel structures – Part 1-8: Design of joints) to prevent chord distortional failure and the evaluation of nominal stresses from fatigue actions (Section 4(1) of EN 1993-1-9, Eurocode 3: Design of steel structures – Part 1-9: Fatigue). Regarding design of aluminum structures (Eurocode 9: Design of aluminum structures – Part 1-1: General structural rules and Part 1-3: Structures susceptible to fatigue), some general rules account for distortion and distortional buckling without any specific guidance. Finally, special attention is paid for distortion by the following bridge design specifications. The Guide Specifications for Horizontally Curved Highway Bridges by the American Association of State Highway and Transportation Officials – AASHTO (1993) [53] specify the maximum spacing of the intermediate diaphragms through an approximate formula. This provision meets the requirement that the distortional warping normal stress is limited within 10% of the bending normal stress and the transverse bending normal stress is limited to 137.3 MPa or lower. The Hanshin Expressway Public Corporation of Japan provides the Guidelines for the Design of Horizontally Curved Girder Bridges – HEPCJ (1988) [54] specifying the maximum spacing of the intermediate diaphragms in curved box girder with respect to that in straight box girders multiplied by a reduction factor, which is equal to unity for a span length less than 60 m. It should be noted here that the boundary conditions and the cross section shape are not taken into account directly for both specifications. In this study, the static analysis of horizontally curved beams of arbitrary cross section, loading and boundary conditions including generalized cross sectional warping and distortional effects due to both flexure and torsion is presented. The aim of this paper is to propose a new formulation by enriching the beam’s kinematics both with out-of- and in-plane deformation modes and, thus, take into account both cross section’s warping and distortion in the final 1D analysis of curved members, towards developing GBT further for curved geometries while employing independent warping parameters, which are commonly used in Higher Order Beam Theories (HOBT). The approximating methods and schemes proposed by Dikaros and Sapountzakis in [30,55] are employed and extended in this study. Adopting the concept of end-effects and their exponential decay away from the support [27], appropriate residual strains are added to those corresponding to rigid body movements. Further, applying Hooke’s stress-strain law and employing the equilibrium equations of 3D elasticity, a system of partial differential equations can be derived for each material over the 2D cross section’s domain together with the corresponding boundary conditions. Consequently, a coupled two-dimensional boundary value problem is formulated, with or without considering Poisson ratio. Applying a proper discretization scheme for the cross section, the above mentioned problem will lead to the formulation of an eigenvalue problem which the eigenvalues and the corresponding eigenvectors, for a desired number of modes, can be extracted from. The obtained set of modes contains axial, flexural and torsional modes in order of significance without distinction between them. To avoid

35

the additional effort needed in order to recognize the most significant modes, the iterative local equilibrium scheme described in the work of Dikaros and Sapountzakis in [55] is adopted until the error due to residual terms becomes minimal. The above scheme is initialized by employing a pre-assumed vector which corresponds to rigid body movements of the cross section (the socalled central solution). Together with the warping functions calculated first, the corresponding distortional ones are also obtained and recursively modify the warping functions due to their coupling. With all these additional modes, the beams’ kinematics is enriched and capable of describing accurately the displacement and stress distribution in the beam. The functions derived are evaluated employing 2D BEM [56]. The coefficient matrices containing the geometric properties of the cross section can now be calculated, as it will be later explained. Thus, a set of boundary value problems are formulated with respect to the unknown kinematical components (displacements, rotations and independent parameters), the number of which is defined by the user depending on the accuracy of the results. The linear system formulated is solved using Isogeometric tools [57], either integrated in the Finite Element Method (FEM) or in the Analog Equation Method (AEM) [58], which is BEM based. Employing the principal of virtual work the new equilibrium equations are derived. The results obtained from the beam element will be compared to those obtained from finite 3D solutions and other research efforts. Numerical examples are presented to illustrate the efficiency and the accuracy of this formulation. To the authors’ knowledge, the numerical procedures previously mentioned have not been reported in the literature for the analysis of curved beams including distortional effects. However, Isogeometric analysis (IGA) has been employed during the last years to provide numerical solutions for various engineering problems. Particularly, NURBS basis functions have been used for vibration and buckling analyses of orthotropic plates with different boundary conditions and configurations [59] as well as for the simulation of 2D fracture mechanics problems in piezoelectric materials under dynamic and static coupled electromechanical loads [60]. More recently, the critical buckling parameters and natural frequencies of plates with functionally graded material (FGM) and internal cracks or voids have been derived by combining the extended Isogeometric analysis based on NURBS and Reissner-Mindlin plate theory [61]. Additionally, the static bending, free vibration and buckling behaviors of functionally graded plates have been studied with the aid of higher-order NURBS and a novel simple quasi-3D hyperbolic shear deformation theory (S-Q3HSDT) [62]. During the last year, the elasto-plastic large deformations in 3D solids have been simulated employing the Isogeometric analysis based on Bézier extraction of NURBS which results in a coarser mesh, higher-order continuity and higher accuracy comparing to FEM [63]. The essential features and novel aspects of the present formulation compared with previous ones are summarized as follows. i. The developed beam formulation is capable of the static analysis of spatial curved beams of arbitrary closed or open and thin- or thick-cross section with one plane of constant curvature (either small or great) considering warping and distortional effects as well as transverse loading to the plane of curvature (as is usually the case in practice). ii. The cross sectional analysis is based on an iterative equilibrium scheme which results in a numerical procedure with less computational effort and complexity comparing to traditional eigenvalue analysis reported in the literature for similar problems. Particularly, modes attributed to different structural phenomena can be separated directly and make the supervision of the results easier. In addition to this, the

36

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

iii.

iv.

v.

vi.

data post-processing and the iterative procedure become faster due to the fact that warping and distortional functions are calculated separately. The accuracy level of the numerical method proposed can be decided by the user by setting the desirable number of the modes taken into account and, thus, increasing the number of higher modes added in the final solution. The numerical solution of advanced beam theories and its application to the analysis of horizontally curved beams is based on b-splines and NURBS (Isogeometric Analysis), as in [64–66], offering the advantage of integrating computer aided design (CAD) in the analysis. In addition to this, the order of the basis functions can be adjusted by the user in the proposed formulation. The curved beam is subjected to arbitrary external loading including warping and distortional moments and is supported by the most general boundary conditions including elastic support or restraint. The proposed method employs a BEM approach which requires only boundary discretization for the cross sectional analysis with line or parabolic elements instead of area elements of the FEM solutions, which require the whole cross section to be discretized into triangular or quadrilateral area elements.

To conclude with, the proposed beam formulation needs further improvements in order to account for non-linear large deformations, buckling behavior, variable or pre-twisted cross sections and more than one plane of curvature. Additionally, non-aligned beam members and connected curved beams with different cross sections are subject of authors’ further research. 2. Statement of the problem Let us consider a curved prismatic element of length L and radius R (Fig. 1). In the general case, an arbitrarily shaped composite cross section of m ¼ 1; 2; . . . ; M homogenous, isotropic and linearly elastic materials with modulus of elasticity Em , shear modulus Gm and Poisson ratio mm , occupying the region Xm of the yz plane with finite number of inclusions is studied. Let also the boundaries of the regions Xm be denoted by CXm (Fig. 2). This

boundary curve is piecewise smooth, i.e. it may have a finite number of corners. In Fig. 2 CXYZ is the principal bending coordinate system through the cross section’s centroid C, while yC , zC are its coordinates with respect to Sxyz reference coordinate system through the cross section’s shear center S. It holds that yC ¼ y  Y and zC ¼ z  Z. The initial radius of curvature, denoted by R is considered constant and it is parallel to Y axis. The displacement vec ðx; y; zÞ of an arbitrary point of the cross section for a specific tor u material is obtained as the sum of SV solution vector corresponding to the rigid body motion combined with a residual (index R) displacement vector due to end-effects which are responsible for the generation of self-equilibrating stress distributions:

 R ðx; y; zÞ  ðx; y; zÞ ¼ u SV ðx; y; zÞ þ u u

Xm a ðxÞW i ðy; zÞ ¼ uðxÞ þ hY ðxÞZ  hZ ðxÞY þ i¼1 i |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} rigid body mov ement

ð1aÞ

out-of-plane warping

v ðx;y;zÞ ¼ v SV ðx;y;zÞ þ v R ðx;y;zÞ ¼ |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} v ðxÞ  zhx ðxÞ

þ

rigid body mov ement

Xm a ðxÞ;x DY i ðy;zÞ i¼1 i |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} distortion in Y direction

Xm   SV ðx;y;zÞ þ w  R ðx;y;zÞ ¼ wðxÞ þ yhx ðxÞ þ wðx;y;zÞ ¼w a ðxÞ;x DZ i ðy;zÞ i¼1 i |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} rigid body mov ement

distortion in Z direction

ð1b;cÞ

where ðÞ;j is for differentiation with respect to j, i is the number of , v , w  are the higher order cross sectional functions considered, u axial, transverse and radial beam displacement components with respect to the Sxyz system of axes, respectively, Wðy; zÞ is the warping function, DYðy; zÞ and DZðy; zÞ are the distortional functions of the in-plane deformation mode Dðy; zÞ while aðxÞ is a function describing the decay of deformation along beam length. Moreover, v ðxÞ and wðxÞ describe the deflection of the center of twist S, while uðxÞ denotes the ‘‘average” axial displacement of the cross section. hx ðxÞ is the angle of twist due to torsion, while hY ðxÞ and hZ ðxÞ are the angles of rotation due to bending about the centroidal Y, Z axes, respectively. The derivation of rigid body motions is in more detail explained in the work of Kang and Yoo (1994), while sin hx  hx ; cos hx  1 assumption is adopted and higher order terms are neglected in this study. Considering the fact that end-effects decay exponentially away from the support, aðxÞ ¼ ecx where c is a constant to be specified. However, different expressions of this

Fig. 1. Prismatic curved beam under indicative axial-flexural-torsional-distortional loading of an arbitrary homogenous cross section occupying the two dimensional region X.

37

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

½clm ðr2 DYÞm þ cðlm þ km ÞðrDÞm;y  cðlm þ km ÞðW ;y Þm þ c3 lm ðDYÞm ecx ¼ 0

ð4bÞ

½clm ðr2 DZÞm þ cðlm þ km ÞðrDÞm;z  cðlm þ km ÞðW ;z Þm þ c3 lm ðDZÞm ecx ¼ 0

ð4cÞ

Thus, the following equations need to be satisfied

  2l þ km l þ km ðr2 WÞm ¼ c2  m ðWÞm  m ðrDÞm

lm

lm

ð5aÞ

ðr2 DYÞm þ

lm þ km ½ðrDÞm;y þ ðW ;y Þm  ¼ c2 ½ðDYÞm  lm

ð5bÞ

ðr2 DZÞm þ

lm þ km ½ðrDÞm;z þ ðW ;z Þm  ¼ c2 ½ðDZÞm  lm

ð5cÞ

Together with the following boundary conditions

 Fig. 2. Arbitrary composite cross section of m homogenous materials occupying the two dimensional region X.

parameter have also been adopted in other research efforts (i.e. polynomials of various degrees). It should be noted here that, usually, within GBT the first derivatives of the parameters multiplying the corresponding distortional modes are selected as dependent warping parameters. This leads to a quadratic generalized eigenvalue problem, thus, requiring the doubling of the unknowns, which is avoided with the proposed formulation. After establishing the displacement field, the strain components for mth material due to end-effects can be computed as

ðrxn Þm ¼ ðrxy Þm ny þ ðrxz Þm nz ¼ 0 on ðCX Þm ðrxn Þm ¼ ðrxn Þm



ðryn Þm ¼ ðryy Þm ny þ ðryz Þm nz ¼ 0 on ðCX Þm ðryn Þm ¼ ðryn Þm



on ðCX Þm [ ðCX Þn m – n

ð6aÞ

ð6bÞ

ð6cÞ

Employing the relation mm ¼ km =½2ðkm þ lm Þ and expanding the stresses in the above boundary conditions, the following boundary value problem is formulated

ð2aÞ

 ðr2 WÞm ¼ c2 

ðeyy Þm ¼ ðv R ðx; y; zÞ;y Þm ¼ a;x ðDY ;y Þm

ð2bÞ

(

 R ðx; y; zÞ;z Þ ¼ a;x ðDZ ;z Þm ðezz Þm ¼ ðw m

ð2cÞ

2 1 þ mem ðWÞm  ðrDÞm 1  mem 1  mem



ð7aÞ

ðW ;n Þm ¼ c2 ½ðDYÞm ny  ðDZÞm nz  on free surface g m ðW ;n Þm þ g n ðW ;n Þn ¼ c2 ðg m  g n Þ½ðDYÞm ny  ðDZÞm nz  on Interfaces

ð7bÞ

R ðx; y; zÞ;y Þ ðcxy Þm ¼ ðcyx Þm ¼ ðv R ðx; y; zÞ;x Þm þ ðu m ¼ a;xx ðDYÞm þ aðW ;y Þm

ð2dÞ

 R ðx; y; zÞ;x Þ þ ðu  R ðx; y; zÞ;z Þ ðcxz Þm ¼ ðczx Þm ¼ ðw m m ¼ a;xx ðDZÞm þ aðW ;z Þm

ð2eÞ

 R ðx; y; zÞ;y Þ þ ðv R ðx; y; zÞ;z Þ ðcyz Þm ¼ ðczy Þm ¼ ðw m m ¼ a;x ½ðDZ ;y Þm þ ðDY ;z Þm 

ð2fÞ

Employing the well-known stress-strain constitutive relationship for elastic media and isotropic solid, the stress components of the mth material are as follows in matrix form

6 ðrij Þm ¼ 4

on ðCX Þm [ ðCX Þn m – n

ðrzn Þm ¼ ðrzy Þm ny þ ðrzz Þm nz ¼ 0 on ðCX Þm ðrzn Þm ¼ ðrzn Þm

 R ðx; y; zÞ;x Þ ¼ a;x ðWÞm ðexx Þm ¼ ðu m

2

on ðCX Þm [ ðCX Þn m – n

ktr½e þ 2lexx

2lcxy

2lcxz

2lcyx

ktr½e þ 2leyy

2lcyz

2lczx

2lczy

ktr½e þ 2lezz

ðr2 DYÞm þ

1 þ mem ½ðrDÞm;y þ ðW ;y Þm  ¼ c2 ½ðDYÞm  1  mem

ð7cÞ

ðr2 DZÞm þ

1 þ mem ½ðrDÞm;z þ ðW ;z Þm  ¼ c2 ½ðDZÞm  1  mem

ð7dÞ

8 g m fðDY ;n Þm þ ½ðDY ;y Þm ny þ ðDZ ;y Þm nz g > > > > þg  ½ðDY Þ þ ðDZ Þ n ¼ g  ðWÞ n on free surface > > ;y m ;z m y m m y < m

g m fðDY ;n Þm þ ½ðDY ;y Þm ny þ ðDZ ;y Þm nz g þ g m ½ðDY ;y Þm > > > þðDZ ;z Þm ny þ g n fðDY ;n Þn þ ½ðDY ;y Þn ny þ ðDZ ;y Þn nz g þ g n ½ðDY ;y Þn > > > : þðDZ ;z Þn ny ¼ ðg m  g n ÞðWÞm ny on Interfaces ð7eÞ

3 7 5

ð3Þ

m

where i; j ¼ x; y; z, ðtr½eÞm ¼ ðexx þ eyy þ ezz Þm and km ; lm are the two Láme parameters of the mth material. In case poisson ratio mm ¼ 0, k ¼ 0 and Em ¼ 2l. Employing local equilibrium equations of threedimensional elasticity considering body forces to be absent, substituting stress components (Eq. (3)) and the exponential function aðxÞ, the following system of partial differential equations can be derived for the mth material as

½lm ðr2 WÞm þ c2 ðlm þ km ÞðrDÞm þ c2 ð2lm þ km ÞðWÞm ecx ¼ 0 ð4aÞ

8 g m fðDZ ;n Þm þ ½ðDY ;z Þm ny þ ðDZ ;z Þm nz g þ g m ½ðDY ;y Þm > > > > þðDZ Þ n ¼ g  ðWÞ n on free surface > > ;z m z m m z < g m fðDZ ;n Þm þ ½ðDY ;z Þm ny þ ðDZ ;z Þm nz g þ g m ½ðDY ;y Þm > > > þðDZ ;z Þm nz þ g n fðDZ ;n Þn þ ½ðDY ;z Þn ny þ ðDZ ;z Þn nz g > > > :  þg n ½ðDY ;y Þn þ ðDZ ;z Þn nz ¼ ðg m  g n ÞðWÞm nz on Interfaces ð7fÞ where m ¼ mm =ð1  mm Þ is the effective Poisson ratio while g m ¼ lm =lref , g m ¼ km =lref are weighted elastic constants with respect to lref which is the shear modulus of reference material. If a plane stress assumption is employed, mem is substituted by mm . When mm ¼ 0 it holds that g m ¼ Em =Eref , g m ¼ 0, with Eref being the e m

38

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

elastic modulus of reference material, and the aforementioned boundary value problem is simplified. Therefore, employing a proper discretization for the cross section, the above coupled boundary value problem (Eqs. (7)) will lead to the formulation of a generalized eigenvalue problem of the form AF ¼ c2 BF where A; B are known coefficient matrices, c is the eigenvalue and F ¼ ½ W DY DZ T is the eigenvector of the problem. The solution of eigenvalue problem yields a set of eigenvalues together with the corresponding eigenvectors which constitute a basis of cross sectional deformation modes suitable for distortional analysis of beams. As mentioned earlier, the iterative equilibrium scheme described by Ferradi, Cespedes and Arquier (2013) in [28] as well as Dikaros and Sapountzakis (2014) in [30] is employed here until a sufficient number of modes is obtained to represent accurately the non-uniform warping effects and the corresponding distortional ones. In order to initialize the above stated boundary value problem, the rigid body movements of the cross section are employed. These correspond to SV flexural and torsional warping modes. Afterwards, in order to restore equilibrium the secondary warping modes are determined together with their corresponding distortional ones. Following this concept, the iterative procedure is formulated converging to the exact shape of the warping in a section. Each functional vector F iþ1 has to fulfil the orthogonality condition with respect to the functions F i corresponding to the previous set of modes. Knowing that each mode is computed with respect to c2i , it follows that

R c2i

¼R

Xm Xm

After

F 2i dX

ð8Þ

F i F iþ1 dX the

evaluation

of

this

constant,

the

normalized

e F iþ1 ¼ can be established and the functions can be obtained. Within the context of the above considerations and considering up to secondary warping as well as distortional displacements, the enriched kinematics of an arbitrary point of the beam for mth material is given as c2i F iþ1

larly, fx ðxÞ, vx ðxÞ are the independent distortional parameters introduced to describe the nonuniform distribution of primary and secondary distortion due to torsion, while fY ðxÞ, fZ ðxÞ, vY ðxÞ, vZ ðxÞ are the independent distortional parameters introduced to describe the nonuniform distribution of primary and secondary distortion due to flexure. All these parameters are multiplied by the corresponding warping and distortional functions which are components of the Wðy; zÞ and Dðy; zÞ vectors derived by the solution of the coupled boundary value problem stated in Eqs. (7). In Eqs. (9), 16 degrees of freedom have been employed in 3D space. These activate 12 cross sectional deformation modes, namely rigid (4), primary (4) and secondary motions (4), including extension. If tertiary displacements have to be employed for accuracy reasons, the beam’s kinematics is enriched further as follows:

 S ðx; y; zÞ  ðx; y; zÞ ¼ u P ðx; y; zÞ þ u u ¼ uðxÞ þ hY ðxÞZ  hZ ðxÞY þ gx ðxÞ/PS ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} primary

þ gY ðxÞ/PCY ðy; zÞ þ gZ ðxÞ/PCZ ðy; zÞ þ nx ðxÞ/SS ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} secondary

þ

nY ðxÞ/SCY ðy; zÞ

þ nZ ðxÞ/SCZ ðy; zÞ þ xx ðxÞ/TS ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

ð10aÞ

tertiary

v ðx; y; zÞ ¼ v ðxÞ  zhx ðxÞ þ fx ðxÞv PS ðy; zÞ þ fY ðxÞv PCY ðy; zÞ þ fZ ðxÞv PCZ ðy; zÞ

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} primary

þ vx ðxÞv þ vY ðxÞv SCY ðy; zÞ þ vZ ðxÞv SCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} S S ðy; zÞ

secondary

þ wx ðxÞv TS ðy; zÞ þ wY ðxÞv TCY ðy; zÞ þ wZ ðxÞv TCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} tertiary

 y; zÞ ¼ wðxÞ þ yhx ðxÞ wðx; þ fx ðxÞwPS ðy; zÞ þ fY ðxÞwPCY ðy; zÞ þ fZ ðxÞwPCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} primary

 ðx; y; zÞ ¼ u  P ðx; y; zÞ þ u  S ðx; y; zÞ u

þ vx ðxÞwSS ðy; zÞ þ vY ðxÞwSCY ðy; zÞ þ vZ ðxÞwSCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

¼ uðxÞ þ hY ðxÞZ  hZ ðxÞY þ gx ðxÞ/PS ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

secondary

primary

þ gY ðxÞ/PCY ðy; zÞ þ gZ ðxÞ/PCZ ðy; zÞ þ nx ðxÞ/SS ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

ð9aÞ

þ wx ðxÞwTS ðy; zÞ þ wY ðxÞwTCY ðy; zÞ þ wZ ðxÞwTCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} tertiary

ð10b; cÞ

secondary

v ðx; y; zÞ ¼ v ðxÞ  zhx ðxÞ þ fx ðxÞv PS ðy; zÞ þ fY ðxÞv PCY ðy; zÞ þ fZ ðxÞv PCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} primary

þ vx ðxÞv SS ðy; zÞ þ vY ðxÞv SCY ðy; zÞ þ vZ ðxÞv SCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} secondary

 y; zÞ ¼ wðxÞ þ yhx ðxÞ wðx; þ fx ðxÞwPS ðy; zÞ þ fY ðxÞwPCY ðy; zÞ þ fZ ðxÞwPCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} primary

þ vx ðxÞwSS ðy; zÞ þ vY ðxÞwSCY ðy; zÞ þ vZ ðxÞwSCZ ðy; zÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} secondary

ð9b; cÞ P , u  S , denote the primary and secondary longitudinal diswhere u placements, respectively. gx ðxÞ, nx ðxÞ are the independent warping parameters introduced to describe the nonuniform distribution of primary and secondary torsional warping, while gY ðxÞ, gZ ðxÞ are the independent warping parameters introduced to describe the nonuniform distribution of primary warping due to shear. Simi-

In this case 22 degrees of freedom have been employed in order to describe the beam’s behavior. The additional 6 degrees, namely nY ðxÞ, nZ ðxÞ, xx ðxÞ, wx ðxÞ, wY ðxÞ and wZ ðxÞ, account for 3 tertiary warping and 3 tertiary distortional effects, respectively. These activate 4 additional cross sectional deformation modes including extension. Particularly, the additional displacements in Eqs. (10) are expressed as the sum of 2D functions (/SCY , /SCZ , /TS , v TS , v TCY , v TCZ , wTS , wTCY , wTCZ ) multiplied by the additional parameters (additional degrees of freedom mentioned previously) which define their longitudinal intensity. These are residual displacements due to end-effects which generate self-equilibrating stress distributions [67,68]. The enrichment of the beam’s kinematics can be done automatically by increasing the number of modes, which are an input value for the boundary value problem to be solved. This results in the evaluation of additional cross sectional operators which will be employed in the analysis of the beam model, after establishing the strain components as it will be described in the following. After establishing the displacement field, the linear straindisplacement relations in the system ðx; y; zÞ can be written for the mth material as follows

39

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

    v R  ;x  ðexx Þm ¼ u ;  RY R

ðeyy Þm ¼ v ;y ;



 ;z ðezz Þm ¼ w

ð11aÞ

    R u  ;y ; ðcxz Þm þu  RY R   R  ;z ; ðcyz Þ ¼ w  ;x   ;y þ v ;z þu ¼w m RY

½kl m ¼



1 6 60 6 60 6 ½Aux1m ¼ 6 60 6 6 40 0

0 0 0 1 0 0

0 0 0 0 1 0

0 0 0 z y 0

Z 0 0 0 0 0

Y 0 0 0 0 0

/PS 0 0 0 0 0

/PCY 0 0 0 0 0

/PCZ 0 0 0 0 0

/SS 0 0 0 0 0

0 0 0

0 0 0

v v

0 0 0

v

0 0 0

0 0 0

v v

P S wPS

P CY wPCY

P CZ wPCZ

S S wSS

S CY wSCY

0

0

0

0

0

0 0 0

k11 ¼

0 6 60 6 6 60 6 6 1 ½Aux2m ¼ 6 6 R eðRÞ 6 6 6 6 60 4 0

 1R eðRÞ 0 Rz eðRÞ 0 0 0 0 0

0 0

0 0

0 0

0

0

0 0

Z eðRÞ R

1 

ð12fÞ

k12 ¼ k21 ¼

Y eðRÞ R

k22 ¼

1 ½Aux1Tm ½Cm ½Aux2m eðRÞ dX

Xm

1 ½Aux2Tm ½Cm ½Aux1m eðRÞ d

Xm

1 ½Aux2Tm ½Cm ½Aux2m eðRÞ dX

ð13Þ

X;

M X R m¼1

From Eq. (12e) after integrating by parts, it holds that

NQM ¼ k22 u;x þ k21 u

ð14Þ

where NQM is the vector of the stress resultants at the beam’s ends. Moreover, the external work can be derived as follows

0 0

v PS

0 0

0

0

0

0

/PCY R

/PCZ R

/SS R

eðRÞ

Xm

M X R

/PS R

eðRÞ

1 ½Aux1Tm ½Cm ½Aux1m eðRÞ dX;

m¼1

3

0 0

Xm M X R m¼1

7 7 7 7 7 7eðRÞ S 7 v CZ 7 7 wSCZ 5 0

0 0

M X R m¼1

ð12aÞ

2

k22

beam for constant radius of curvature, [C]m is the elasticity matrix employed to derive stresses r, dU is the virtual strain energy and [k11, k12, k21 and k22]m are 16  16 coefficient matrices containing the geometric properties of the cross section for the mth material. These are calculated as follows

ð11bÞ

R where RY is set as eðRÞ in the following and introduces the thickness-curvature effect of the curved beam. Employing the expressions of the displacement components (Eqs. (9)), the strains and stresses can be computed. Applying the principle of virtual work or any other variational principle following standard arguments in the calculus of variations, the governing differential equations for the beam in terms of the kinematical components can be derived. Thus, the local stiffness matrix ½kl  of the spatial curved beam can be evaluated after solving a system of linear equations. Finally, the matrix form of stiffness matrix is derived as follows

2

k12

k21

where [Aux1]m, [Aux2]m are auxiliary matrices to express strains e in 1 dXdx ¼ dV is the differential volume of the curved matrix form, eðRÞ

v ;x þ

ðcxy Þm ¼

k11



eðRÞ

R

eðRÞ

v PCY

v PCZ

eðRÞ

R P CY;y

v SS

eðRÞ

R P CZ;y

R

eðRÞ

v SCY

eðRÞ

R S CY;y

v SCZ

3 eðRÞ

R S CZ;y

v PS;y

v

v

v SS;y

v

v

wPS;z

wPCY;z

wPCZ;z

wSS;z

wSCY;z

wSCZ;z 0

eðRÞ . . . . . . 0

0

0

0

0

þ/PS;y

þ/PCY;y

þ/PCZ;y

þ/SS;y

0

0

0

0

0

0

wPS;y

wPCY;y

wPCZ;y

wSS;y

wSCY;y

wSCZ;y

0

0 0

1

0

/PS;z

/PCY;z

/PCZ;z

/SS;z

0

0 0

0

0

0

0

0

0

þv PS;z

þv PCY;z

þv PCZ;z

þv SS;z

þv SCY;z

þv SCZ;z

7 7 7 7 7 7 7 7 7 7 7 7 7 7 5

ð12bÞ

2

em ¼ ½Aux1m u;x þ ½Aux2m u

ð12cÞ

rm ¼ ½Cm ½Aux1m u;x þ ½Cm ½Aux2m u

ð12dÞ

Z dU ¼

M Z LX

Xm

0 m¼1

ðduT;x ½Aux1Tm

þ

Z ) dU ¼ 0

0

0

0 wPS wPCY wPCZ wSS wSCY wSCZ

0

ð15aÞ M Z X

" ðdu

T

½AuxTm tÞdx

Cm m¼1|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

þ

M Z X

m¼1

Xm

#L ðdu

T

½AuxTm tÞd

X

ð15bÞ 0

where t is the traction vector applied on the lateral surface of the beam including the end cross sections and p is the external load vector of the beam. Combining Eqs. (12e) with (15b), the expression of the variational total potential energy can be evaluated and, thus, the governing differential equations of the problem can be obtained together with the boundary conditions as

ðduT;x k22 u;x þ duT k12 u;x þ duT;x k21 u þ duT k11 uÞdx

) by parts integration Z L ðduT fk22 u;xx þ ½k12  k21 u;x þ k11 ugÞdx ) dU ¼ 0

þ ½duT fk22 u;x þ

0 0

duT p

1 dXdx eðRÞ

L

L k21 ug0

0 0 1 y

dW ¼

duT ½Aux2Tm Þ

 ð½Cm ½Aux1m u;x þ ½Cm ½Aux2m uÞ

3 1 0 0 0 Z Y /PS /PCY /PCZ /SS 0 0 0 0 0 0 6 7 P P P S S S ½Auxm ¼ 4 0 1 0 z 0 0 0 0 0 0 v S v CY v CZ v S v CY v CZ 5

ð12eÞ

k22 u;xx þ ½k12  k21 u;x þ k11 u ¼ p

ð16aÞ

a1 u þ a2 NQM ¼ a3

ð16bÞ

40

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

where ai are diagonal matrices and vector containing known coefficients according to the boundary conditions of the beam (i.e. for clamped end a1 ¼ 1 and a2 ¼ a3 ¼ 0). Employing the expressions of the displacement components in Eqs. (10), the cross sectional operators (Eq. (13)) and the governing differential equations of the curved beam can be obtained in a similar way when tertiary or higher warping and distortional effects are considered. 3. Numerical solutions for curved beams

2

d ui ðxÞ 2

dx

¼ qi ðxÞ

where ui ðxÞ are the different kinematical components and qi ðxÞ are the corresponding fictitious loads (i ¼ 1; . . . ; N with N depending on the number of modes employed). B-splines have been employed to approximate the fictitious load curve. The fundamental solution of each Eq. (17) is a partial solution of the following differential equation 2

d ui ðx; nÞ 2

3.1. A BEM method for the cross sectional analysis The evaluation of the warping and distortional functions is accomplished by solving the problems described by Eqs. (7). Warping functions W and their derivatives are at first computed by solving Eqs. (7a) and (7b). Afterwards, these values are inserted as generalized body forces in Eqs. (7c)–(7f) which are solved as a 2D elasticity problem in order to obtain distortional functions D. The solution of the problem is accomplished employing BEM within the context of the method of subdomains and BEM for Navier operator [56,69]. Afterwards, the values are normalized through the constant c2i , as earlier described, and the procedure is repeated for the desired number of modes. Finally, the functions calculated are employed in order to obtain the cross sectional operation factors given in Eq. (13). These are used as input values together with the elasticity and function matrices to solve the curved beam model with the methods described below. 3.2. Isogeometric tools for the 1D analysis of curved beams either with AEM or FEM

ð17Þ

dx

¼ dðx  nÞ

ð18Þ

where ui ðx; nÞ and its derivatives are given as 

dui ðx; nÞ 1 ¼ sgnr dx 2

K1 ðx; nÞ ¼

ð19aÞ

1 K2 ðx; nÞ ¼ ui ðx; nÞ ¼ r 2

ð19bÞ

with r ¼ x  n; x; n points of the beam. Employing this fundamental solution, the integral representation of the kinematical components is obtained as

Z

L

ui ðnÞ ¼ 0

 L du ðxÞ K2 ðx; nÞqi ðxÞdx  K2 ðx; nÞ i  K1 ðx; nÞui ðxÞ dx 0

ð20Þ

i ðxÞ Eq. (20) implies that if qi ðxÞ and all boundary values (ui ðxÞ, dudx ) at the beam ends 0, L are known, ui ðxÞ can be calculated at any internal point of the beam. Differentiating Eq. (20), the expressions for the derivatives of ui ðxÞ are derived as

dui ðnÞ ¼ dn

Z

L

0

 L du ðxÞ K1 ðx; nÞqi ðxÞdx  K1 ðx; nÞ i dx 0

ð21Þ

In traditional FEM or AEM analysis, the design and geometric models never coincide due to the fact that even though they are both representations of a true object, they rely on different basis functions. This, in turn, produces concerns related to accuracy in the computations, particularly for curved structures. If b-splines or NURBS basis functions are used, their properties are inherited by the FEM or the AEM model. This is particularly important because it allows the circumvention of certain serious difficulties in developing the numerical model (e.g. advanced beams in which, except for bending, shear and torsion, higher order phenomena must be considered making the mesh processing more complicated). Moreover, as the shape functions are smoother, the error affecting its derivatives becomes smaller. This is important especially in FEM models because stress fields are not the primary solution variables, but need to be computed by differentiating displacements through post-processing techniques. Thus, smoother displacement fields ensure more accurate approximation of the stresses. This is not an aspect to be considered in the AEM models due to the fact that fictitious loads, which are the second derivatives of the unknowns, are approximated by shape functions and calculated at first.

Then, qi ðxÞ is either approximated with constant elements or quadratic elements [1]. The introduction of b-splines in the above mentioned expressions can be done by substituting qi ðxÞ with the polynomial representation of a quadratic b-spline with a uniform knot vector N with n 2 [0, 1], which is the parameter space similar to the classic FE subdivision. The first and last knot values are repeated depending on the b-spline degree p and their multiplicity is usually p þ 1. In one dimension, basis functions formed are interpolatory at the ends of the parameter space interval (knot vector with multiplicities). However, nonuniform knot vectors and repeated knots can also be used with NURBS. According to [70], the NURBS basis functions can be expressed in terms of b-splines basis defined by the Cox-De Boor recursive formula

3.2.1. A BEM-based method (AEM) combined to Isogeometric analysis The evaluation of the kinematical components ui (either of those in Eqs. (9) or (10)) is accomplished using the AEM [58]. These have continuous derivatives up to the second order with respect to x at the interval ð0; LÞ and up to the first order at x ¼ 0; L satisfying the boundary value problem described by the coupled governing differential equations of equilibrium in Eq. (16a) along the beam and the boundary conditions in Eq. (16b) at the beam ends x ¼ 0; L. According to this method, for these functions, the following relation is valid

Thus, the basis functions can be derived for the quadratic b-spline [65,66]. The quadratic b-spline curve is defined by



Ni; 0ðnÞ ¼

1 if ni 6 n < ni þ 1 0

otherwise

 ;

p¼0

n  ni Ni; ni þ p  ni ni þ p þ 1  n Ni þ 1; p  1ðnÞ; p  1ðnÞ þ ni þ p þ 1  ni þ 1

ð22aÞ

Ni; pðnÞ ¼

C i ðnÞ ¼

pP1

2 X Ni;2 ðnÞPii

ð22bÞ

ð23Þ

i¼0

where Pii are the control points P0i , P1i and P 2i of the initial control polygon [65]. Substituting basis functions to Eq. (23), the expression of the fictitious load qi ðxÞ is derived as

qi ðxÞ ¼ P0i  2xP0i þ x2 P0i þ 2xP1i  2x2 P1i þ x2 P2i

ð24Þ

41

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

Three equidistant collocation points have been used on the longitudinal axis of the beam at first. Considering the cubic b-spline, the expression for the fictitious load qi ðxÞ is derived as 3

qi ðxÞ ¼

ðl  xÞ

3xðl  xÞ

P 0i þ

3

l

l

3

2

P1i þ

3x2 ðl  xÞ l

3

P2i þ

x3 l

3

P3i

ð25Þ

where P ii are the control points P0i , P 1i , P2i and P3i . Similarly, the expression for the fictitious load qi ðxÞ for the quartic B-spline is derived as

qi ðxÞ ¼

ðL  xÞ4 þ

P0i þ

L4 4x3 ðL  xÞ L4

4xðL  xÞ3

L4 x4 P3i þ 4 P4i L

P1i þ

6x2 ðL  xÞ2 L4

P2i

½E11  ½E12  ½E21  ½E22 



fu1i g fu2i g



 ¼

½a10 þ a20 k21 u0 þ ½a20 k22 u0;x ¼ fa30 g for ½a1L þ a2L k21 u þ ½a L

L 2L k22 u;x

¼ fa3L g for

x¼0 x¼L

ð26Þ

   f0g þ ! ½Efui g ¼ fDg þ fT i g f0g fT 2i g

fD1 g

ð27Þ where ½E is a 64  64 square matrix, fui g, fDg, fT i g are 64  1 vectors, ½E11 , ½E12  are 32  32 known coefficient matrices and fD1 g is known 32  1 coefficient vector with their dimensions depending

ð28a; bÞ

and iterating over the number of boundary conditions, which is 16 for each end in this case. ½kij ; i; j ¼ 1; 2 are the geometric constant matrices given in Eq. (13) and inserted as input in the beam analysis. Thus, the following system is derived

½Efui g ¼ fDg þ ½FfPi g ! fui g ¼ ½E1 fDg þ ½E1 ½FfPi g

where P ii are the control points P0i , P 1i , P2i , P 3i and P 4i . Eqs. (20) and (21) written for the boundary points constitute a system of four simultaneous integral equations, while the boundary conditions in Eq. (16b) are formulated in matrix form giving four more equations. Combining the aforementioned equations, the following system is derived for the case described in Eqs. (9)



on the number of modes employed and all explained in the work of [65]. Particularly in this study, ½E11 , ½E12  are calculated writing Eq. (16b) in the following matrix relations

ð29Þ

With known values at beam ends and applying the integral representations in Eqs. (20) and (21) at the internal collocation points, the following can be calculated

fU i g ¼ ½AfP i g þ ½Cfui g

ð30aÞ

fU 0i g ¼ ½A0 fPi g þ ½C 0 fui g

ð30bÞ

ðU i ; U 0i Þ

where are the vectors containing the values of the different kinematical components and their first derivative at internal collocation points. The coefficients of the square matrices ½A; ½A0  are given by the numerical solution of the integrals given in the work of [65] with respect to the vector of the control points fPi g. The ½C; ½C 0 ; ½C 0 ; ½C 000  matrices are also known coefficient matrices. Substituting Eq. (29) into Eqs. (30a) and (30b), the following are derived

fU i g ¼ ½BfPi g þ fRg

Fig. 3. Doubly symmetric cross sections of example 1 with their tip loads and plan views on the right ((a) Open-shaped and (b) Close-shaped).

ð31aÞ

42

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

fU 0i g ¼ ½B0 fPi g þ fR0 g

ð31bÞ

where

½B ¼ ½½A þ ½C½E1 ½F 0

0

0

1

½B  ¼ ½½A  þ ½C ½E ½F 1

ð32aÞ

fPg ¼ fpg  ð½K 12  K 21 fR0 g þ ½K 11 fRgÞ ð32bÞ

fRg ¼ ½C½E fDg

ð32cÞ

fR0 g ¼ ½C 0 ½E1 fDg

ð32dÞ

Afterwards, the element cross sectional operators are assigned to each discretization element in order to calculate the local stiffness matrix of the beam element employing Eq. (16a) transformed into matrix form as follows

½K l  ¼ ½K 22  þ ½K 12  K 21 ½B0  þ ½K 11 ½B

ð33Þ

where ½K 11 ; ½K 12 ; ½K 21 ; ½K 22  are the geometric constant (16XNdof) X (16XNdof) matrices formulated for all of the discretization elements through an iterative procedure with Ndof being their number. All matrices in Eq. (33) are square with the same dimensions. The stiffness matrix is formulated either with respect to the values of the

Fig. 4. Total translation (cm) contours of (a) a 3D solid FEM model with 7875 quadrilateral elements and (b) the proposed one (AEM with 10 cubic splines or 100 constant elements).

Fig. 5. Normal stress constant elements).

fictitious loads in the case of the original AEM or the control points when b-splines are integrated in the AEM. Additionally, substituting Eqs. (31) into Eq. (16a), the load vector applied is written as follows in matrix form

ð34Þ

where fpg is a load vector containing the load values along the beam either concentrated or uniformly distributed. Finally, the global equation system can be formulated and the unknowns can be evaluated. The result is not the displacement vector as in traditional FEM but it is either the fictitious load, which represents the second derivative of the kinematical components, or the control points. Employing Eqs. (31), the kinematical components and their first derivatives can eventually be obtained. These are employed as input values together with the function matrices in a post-processing procedure in order to derive the total displacements, stresses and stress resultants along the curved geometry. Except for degree elevation of b-splines (Eqs. (25) and (26)) or NURBS, knot insertion can be employed as a refinement procedure of the mesh in order to increase accuracy. This process has been described thoroughly in [64].

3.2.2. Numerical solution with FEM and NURBS The problem described in Section 2 can be set in the equivalent variational formulation through the principal of virtual work and finally derive Eqs. (16). In order to discretize Eqs. (16), the displacement field ui ðxÞ given in Eqs. (9) or (10) can be approximated by means of polynomial interpolating functions of p-degree as extensively is described in [71]. Substituting the displacement approximation for each discretization element (n in total), the equilibrium equations in terms of the nodal displacements of the Finite Element (FE) mesh can be expressed. Thus, ½K l  (Eq. (33)) matrix can be obtained by assembling the contributions from the individual elements. This can be performed numerically by using the Gauss quadrature rule. In order to compute Gauss base points and weight factors an algorithm has been employed according to [72]. The same procedure can be followed in BEM or in the AEM technique, with quadrature nodes being the collocation points, if the analytical solution of the integrals employed needs to be avoided for computational reasons.

rxx (kN/cm2) contours for (a) a 3D solid FEM model with 7875 quadrilateral elements and (b) the one proposed (AEM with 10 cubic splines or 100

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

Instead of traditional interpolating functions employed in FEM, p-degree NURBS interpolating functions are employed for the representation of the displacement field ui ðxÞ following a similar procedure with the one described in the previous paragraph. In this case, curve C given in Eq. (23) has a p-degree NURBS representation defined by n X C i ðnÞ ¼ Ri;p ðnÞPii

ð35Þ

i¼1

where Pii are the control points employed for each kinematical component and Ri;p ðnÞ, which are the NURBS basis functions, can be expressed as follows

Ni;p ðnÞwi Ri;p ðnÞ ¼ Pn i¼1 N i;p ðnÞwi

ð36Þ

where N i;p ðnÞ are given in Eqs. (22) and wi (2 R) are weights related to the ith control point and increase the capabilities of the b-splines interpolation [57]. It should be noted that if all weights are equal, then Ri;p ðnÞ ¼ Ni;p ðnÞ and curve C is a b-spline curve. The geometry of the beam is described by a NURBS structure through an initial control polygon given by the following spatial coordinates: ðxi ; yi ; zi ; wi Þ. One or more meshes are defined for the different kinematical components based on the initial NRUBS structure while knot insertion and degree elevation can be employed in order to refine the initial structure directly on these meshes. In addition to these, Gauss points and the corresponding Jacobians are evaluated in parametrical and physical space together with the values of the basis functions, their derivatives and the radius of curvature directly on the meshes. Thus, post-processing becomes easier and no need for additional loops and calculations of basis functions are demanded. ½K l  matrix can finally be obtained.

43

In case of a curved geometry the following relations need to be employed

J ¼ S0 ¼ R¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 02 þ Y 02

J jX 0 Y 00  X 00 Y 0 j

ð37Þ ð38Þ

where ðÞ0 are the derivatives with respect to n, J is the Jacobian of the transformation, S is the arc length of the beam with one plane of curvature, R is the radius of curvature of the beam and X, Y are the coordinates on the plane of curvature. After establishing the kinematical components and their derivatives, total displacements can directly be plotted on the curved geometry. Finally, stresses and stress resultants can also be derived without the need for excessive post-processing. This numerical solution has been employed by authors in previous formulations of curved beams [66] without considering distortional effects and its advantages over traditional FEM and AEM have been highlighted. 4. Numerical examples In order to validate the proposed formulation of the curved beam element described above, investigate the importance of curvature in distortional analysis and examine the advantages attained by the use of the methods proposed in terms of accuracy and computational effort, computer programs have been written and various curved beam models have been studied. The numerical results have been obtained employing NURBS in FEM and constant elements or b-splines for the representation of the AEM fictitious loads. Then, the results are compared to those obtained by the

Fig. 6. Kinematical components gx ðxÞ, fx ðxÞ and hx ðxÞ of the I-shaped beam examined in example 1 (derived by the proposed formulation with 10 cubic NURBS in FEM).

44

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

Table 1 Kinematical components, Shear stresses and Stress resultants of the I-shaped beam examined in example 1.

w ðcmÞ at x = L hx ðradÞ at x = L hY ðradÞ at x = L

" Pz eccentric lateral loading

kN smax ðcm 2 Þ at x = 1/4L xz kN smax xy ðcm2 Þ at x = 1/4L

M Y ðkN cmÞ at x = 0 M t ðkN cmÞ at x = 0 2

M /P ðkN cm Þ at x = 0

AEM 10 cubic b-splines

FEM solid 7875

12.8713 0.2083 0.3274 120.25

12.5466 0.19702 0.3158 111.02

115.23

109.80

330.81 323.01 3401.49

325.05 318.31 –

S

2

36.46

–

2

16.25

–

M /S ðkN cm Þ at x = 0 S

M PDx ðkN cm Þ at x = 0

application of the Finite Element Method (FEM) employing beam, solid (quadrilateral) or plate elements. The computer software [73] has been used for this purpose. In addition to these, up to 1000 boundary elements depending on the cross section type (cross sectional discretization) have been employed in order to evaluate the cross sectional operators and functions. Various cross sections and beam models with different t=d and d=L ratios (t, d and L are the thickness, the width or height depending on the loading plane and the length of the curved close- or open-shaped cross section, respectively) have been employed in order to estimate the magnitude of distortional effects and the number of modes needed in order to derive accurate results.

4.1. Doubly symmetric cross sections

Fig. 7. Convergence of normalized displacements at the free end of the I-shaped beam (example 1).

In this example, in order to validate the proposed formulation and investigate the importance of distortion and curvature in the analysis of curved beam models as well as the number of modes needed, the static problems of two curved cantilever beams with doubly symmetric cross sections (open or closed) are examined (Fig. 3).

Table 2 Kinematical components, Shear stresses and Stress resultants of the Box-shaped beam examined in example 1.

; Pz eccentric lateral loading

w ðmÞ at x = L hx ðradÞ at x = L hY ðradÞ at x = L kN smax ðm 2 Þ at x = 0 xz kN smax ð Þ at x = 0 xy m2

AEM 10 cubic b-splines

FEM solid 2880/FEM plate 960

0.4266 0.0100 0.0131 21,532

0.4316 0.0112 0.0137 20,539 (solid model)

12,602

11,502 (solid model)

2

139,691 127,100 6930.56

139,824 127,324 –

2

1838.93

–

2

610,306

–

2

1999.04

–

M Y ðkN mÞ at x = 0 M t ðkN mÞ at x = 0 M /P ðkN m Þ at x = 0 S

M /S ðkN m Þ at x = 0 S

M PDx M SDx

ðkN m Þ at x = 0 ðkN m Þ at x = 0

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

Fig. 8. Total translation (m) contours of (a) 3D solid FEM model with 2880 quadrilateral elements and (b) the proposed one (AEM with 10 cubic splines or 80 constant elements).

Fig. 9. Normal stresses (kN/m2) rxx (A) and rzz (B) contours for (a) 3D solid FEM model with 2880 quadrilateral elements and (b) the one proposed (AEM with 10 cubic splines or 80 constant elements).

45

The first beam model under consideration has the I-shaped cross section (E ¼ 73; 000 kN=m2 , G ¼ 28; 000 kN=m2 , R ¼ 0:636 m, m ¼ 0:3, t=d ¼ 0:048, d=L ¼ 0:105) shown in Fig. 3a and an arc length of 1 m (Fig. 3a). It is subjected to a concentrated force of 5 kN at its free end eccentrically applied. The displacement field considered is the one described in Eqs. (9). In Fig. 4 the total translation contours are presented for the 3D FEM model (7875 quadrilateral solid elements) and the one proposed (AEM). It is obvious that the results are in agreement. In Fig. 5 the normal stress rxx contours are displayed for the same cases. It is obvious that the models are in coincidence while the use of cubic b-splines significantly reduces the number of discretization elements for the same level of accuracy. In Fig. 6 the distributions of the kinematical components gx ðxÞ, fx ðxÞ and hx ðxÞ are displayed along the length in order to investigate the order of their magnitude and compare them to each other. FEM with NURBS are employed for this purpose. The exponential decay of the primary warping and distortional parameters is illustrated. However, distortion is insignificant. In Table 1 the values of different kinematical components, shear stresses and higher order Moments are compiled for the proposed beam formulation and compared to the FEM solid model when it is possible. It should be noted here that the primary warping Moment is much larger than the primary distortional Moment due to torsion. These quantities cannot be obtained directly by the FEM solid model. Moreover, bending and twisting Moments are similar. Shear stresses are much lower than normal stress rxx (compare values in Table 1 and Fig. 5) and this directly implies that warping is more important than distortion for this cross section according to Eqs. (11a) which are proportional to the corresponding stresses (Eqs. (12c) and (12d)). In addition to these, secondary warping and distortional Moments are much lower than the primary ones. Thus, there is no need to consider additional higher modes in the analysis. Fig. 7 shows the normalized displacements with respect to FEM solid results at the free end of the beam. The proposed solution with either a cubic (4 collocation points) or a quartic (5 collocation points) gives highly accurate results (errors 0–10%). The convergence rate seems to be higher for the knot insertion refinement procedure than that of the degree elevation. Quadratic b-splines give inaccurate results when the least number of collocation points is employed. When employing knot insertion for the cubic b-spline (24 collocation points) the results almost coincide with the FEM solid solution. Finally, considering the same cross section but with an arc length of 3 m, similar conclusions are drawn. The second cantilever beam model studied has a rectangular box-shaped cross section 5.0  3.5 m with plate thickness 0.30 m (E ¼ 3  107 kN=m2 , G ¼ 1:5  107 kN=m2 , R ¼ 25:465 m, m ¼ 0, t=d ¼ 0:085, d=L ¼ 0:087) and an arc length of 40 m (Fig. 3b). The displacement field considered is the one described in Eqs. (9). In Table 2 the values of different kinematical components, shear stresses and higher order Moments are compiled for the proposed beam formulation subjected to a concentrated force of 5000 kN eccentrically applied at its free end. The results are compared to the 3D FEM model (2880 quadrilateral solid elements) or the equivalent plate model (960 plate elements) when it is possible. In contrast to the previous open shaped cross section, primary distortional Moment due to torsion is much larger than the primary warping Moment which can be considered insignificant. The same case stands for the secondary warping Moment. Secondary distortional Moment is not of importance and, thus, there is no need to consider additional modes in the analysis. In addition to these, shear stresses are larger than the previous case comparing to the normal stress rxx . In Fig. 8 the total translation contours are presented for the FEM solid model and the one proposed (AEM with 10 cubic splines or 80 constant elements). In Fig. 9 the normal

46

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

Table 3 Moment ratios for various curvatures of the Box-shaped cross section examined in example 1. M PDx MY

; Pz eccentric lateral loading

R=1 R = 76.394m R = 50.930m R = 38.197m R = 25.465m

at x = 0

0.178 1.141 1.835 2.581 4.365

M PDx Mt

at x = 0

2.844 5.579 5.276 5.101 4.807

MY Mt

at x = 0

15.952 4.891 2.875 1.976 1.101

and exhibit significant values proportionally to normal stress rxx . Finally, in Table 3 three different ratios of Moments have been compiled for various curvatures while length remains the same. It is important to notice that as the curvature increases the M PDx =M Y ratio becomes larger and, thus, distortional effects become more significant. On the other hand, the M Y =M t ratio becomes smaller as it was expected due to the curved geometry and the MPDx =M t ratio is initially increased with the increase of curvature while afterwards slightly reduces in magnitude due to the fact that distortion is steadily important comparing to all other quantities. It seems that a critical value of curvature exists at which distortional effects become much larger. 4.2. Monosymmetric cross sections In this example, in order to further validate the proposed formulation and investigate the importance of distortion and curved geometry in the analysis of beam models as well as the number of modes needed, the static problems of four curved cantilever beams with monosymmetric cross sections (open or closed) are examined. Their cross sections are displayed in Fig. 10. The first beam model has a C-shaped cross section (E ¼ 73; 000 kN=m2 , G ¼ 28; 000 kN=m2 , R ¼ 0:636 m, m ¼ 0:3, t=d ¼ 0:049, d=L ¼ 0:055) and an arc length of 1 m (similar plan view with Fig. 3a). The second model has a trapezoidal box-shaped cross section (E ¼ 3  107 kN=m2 , G ¼ 1:5  107 kN=m2 , R ¼ 25:465 m, m ¼ 0, t=d ¼ 0:086, d=L ¼ 0:086) and an arc length of 40 m (similar plan view with Fig. 3b). The third one has a box-shaped cross section G ¼ 1:39  107 kN=m2 , R ¼ 100 m, (E ¼ 3:25  107 kN=m2 , m ¼ 0:1667, t=d ¼ 0:1, d=L ¼ 0:065) and an arc length of 33 m (similar plan view with Fig. 3b). The last one has the same plan view and cross section shape as previously but a different material for its upper plate (composite cross section), as shown in Fig. 10d in

Fig. 10. Monosymmetric cross sections of example 2.

stresses rxx and rzz contours are displayed for the same cases. It should be noted here that in this case normal stresses rzz and ryy are of more importance comparing to the previous cross section

grey color (E ¼ 4  107 kN=m2 , G ¼ 2  107 kN=m2 , m ¼ 0). The C-shaped cross section is examined for four different load cases, namely concentrically or eccentrically applied radial (parallel to axis of symmetry) and vertical loads, and results are compared to the FEM solid model (4000 quadrilateral solid elements). The proposed formulation is accomplished with 10 cubic b-splines in the AEM and the displacement field is described by Eqs. (9). In Table 4 the maximum values of the total displacement, the normal stress rxx , the shear stress sxy , the primary warping Moment, the secondary warping Moment and the distortional Moment due to torsion are compiled for the various load cases. It should be noted here that when the load is applied radially either concentrically or eccentrically the warping and distortional

47

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

Fig. 11. Kinematical components gx ðxÞ, fx ðxÞ and hx ðxÞ of the trapezoidal Box-shaped (Fig. 10b) beam examined in example 2 (derived by the proposed formulation with 10 cubic NURBS in FEM). Table 4 Total displacement, stresses and higher order Moments of the C-shaped beam (Fig. 10a) for various load cases. AEM 10 cubic b-splines (FEM solid 4000)

umax total ðmÞ at x = L kN rmax ðm 2 Þ at x = 0 xx max kN sxz ðm2 Þ at x = 0

0.127 (0.128) 17.92 (18.63)

4.85 (5.45) 203.01 (247.1)

5.01 (5.99) 217.31 (260.2)

0.31 (0.284) 20.67 (22.56)

4.23 (3.75)

44.75 (38.68)

48.1 (42.44)

5.25 (4.81)

0

5923.16

6591.31

229.68

2

0

1289.18

1171.38

40.5

2

2.3

32.79

34.51

2.41

2

M /P ðkN m Þ at x = 0 S

M /S ðkN m Þ at x = 0 S

M PDx ðkN m Þ at x = 0

Table 5 Total displacement, stresses and higher order Moments of the C-shaped beam (Fig. 10a) for 12 and 16 modes. 12 modes

16 modes

AEM 10 cubic b-splines (FEM solid 4000) umax total ðmÞ at x = L kN rmax ðm 2 Þ at x = 0 xx kN smax ð Þ at x = 0 xz m2

5.01 (5.99) 5.87 (5.99) 217.31 (260.2) 252.01 (260.2) 48.1 (42.44)

43.8 (42.44)

2

6591.31

7413.57

2

1171.38

1286.02

2

34.51

32.03

M /P ðkN m Þ at x = 0 S

M /S ðkN m Þ at x = 0 S

M PDx ðkN m Þ at x = 0

Moments are much lower than the case of the vertically applied load which results in significant primary and secondary warping effects. This indicates that additional modes need to be employed.

Thus, in Table 5 the values of the same quantities are compiled for the proposed formulation employing either 12 or 16 modes. It is important to note that errors regarding total displacement and stresses are reduced from around 12% (12 modes) to less than 5% (16 modes) comparing to the FEM solid model while higher order warping moments are increased by around 10%. The distortional one due to torsion is insignificant, as it is expected. The second cross section is subjected to a vertical concentrated load 10,000 kN concentrically applied at its free end. In Fig. 11 the kinematical components gx ðxÞ, fx ðxÞ and hx ðxÞ are displayed in order to investigate the order of their magnitude and compare them to each other. FEM with NURBS are employed for this purpose. Comparing to the gx ðxÞ and fx ðxÞ distributions of the I-shaped cross section of the previous example (Fig. 6), the corresponding ones for this case exhibit more similar distributions even though they differ in order of magnitude. In Table 6 the values of different kinematical components, stresses and higher order Moments are compiled for the proposed beam formulation, the

48

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

Table 6 Kinematical components, Shear stresses and Stress resultants of the trapezoidal Box-shaped (Fig. 10b) beam examined in example 2.

; Pz eccentric lateral loading

w ðmÞ at x = L hx ðradÞ at x = L hY ðradÞ at x = L kN rmax ðm 2 Þ at x = 0 xx kN smax ðm 2 Þ at x = 0 xy

MY ðkN mÞ at x = 0 Mt ðkN mÞ at x = 0 2

M/P ðkN m Þ at x = 0

AEM 10 cubic b-splines

FEM solid 2714

AEM 50 quad. GWCBT

Timoshenko FEM beam

0.3744 0.0092 0.0120 40,254

0.3547 0.0087 0.0115 38,230

0.3202 0.0067 0.0104 27,633

0.3238 0.0067 0.0106 28,782

24,135

23,085

16,940

3714

254,316 253,683 342,533

254,647 254,647 –

254,820 254,700 244,940

254,648 254,648 0

S

2

36,612

–

30,717

0

2

366,072

–

–

–

2

4127

–

–

–

M/S ðkN m Þ at x = 0 S

MPDx ðkN m Þ at x = 0 MSDx ðkN m Þ at x = 0

3D FEM model (2714 quadrilateral solid elements) when it is possible, the Generalized Warping Curved Beam (GWCB) formulation proposed in [1] (only warping is considered) and the traditional Timoshenko beam element with FEM. It is worth noting here that, comparing to the box-shaped cross section of the previous example, the M PDx =M Y ratio is much lower while primary warping and distortional Moments are of the same order of magnitude in this case. The reasons for these are the difference in the cross section shape and the larger overall volume of this beam. Warping moments are much lower when distortion is not considered in the formulation due to the fact that the two phenomena are coupled (Eqs. (7)). Shear stress of Timoshenko beam element is significantly less than the other cases. Finally, the third cross section is subjected to a vertical concentrated load 5000 kN concentrically applied at its free end. In Table 7 the values of different kinematical components, stresses and higher order Moments are compiled for the proposed beam formulation either curved or straight and the 3D FEM model (6600 quadrilateral solid elements) for validation when it is possible. Additionally, the corresponding values derived by the analysis of the composite cross section showed in Fig. 10d have been compiled for comparison reasons. It should be noted here that the homogenous curved beam exhibits larger distortional effects due to torsion comparing to the straight formulation for which warping and distortional effects due to bending are of importance. Curvature alters significantly the contribution of beam’s resisting mechanisms and makes it more vulnerable to higher phenomena triggered by arising torsion. However, comparing to the previous box-shaped cross sections, the ratio M PDx =M Y is much lower and the beam resists to loading mainly through bending. Considering the corresponding beam element with composite section, it should be noted that its performance is significantly improved through the use of a material with higher stiffness at its upper plate. This becomes obvious through the reduction in stresses and higher order moments due to torsion by more than 50%. However, the primary warping moment due to bending is quite high but still much less than the distortional one (and more desirable than the corresponding quantities due to torsion which causes more brittle failures comparing to bending).

5. Concluding remarks In this paper, the generalized warping and distortional analysis of curved beams is mainly examined. An iterative equilibrium scheme combined with traditional eigenvalue analysis has been developed in order to derive the cross section’s modes that dominate and, thus, permit the isolation of structural phenomena. Various shapes of open or closed cross sections have been considered. Boundary Element Method has been employed for this purpose. Additionally, Isogeometric tools (b-splines and NURBS) integrated in FEM and AEM are applied for the longitudinal analysis of beams allowing for straightforward model handling (i.e. curvature can be easily changed) and reducing significantly computational cost (especially when additional higher order phenomena need to be accounted for). The presented formulation is based on an enriched kinematic field (Eqs. (9) and (10)) taking into account primary, secondary and tertiary higher order phenomena due to both bending and torsion. The importance of the proposed formulation is highlighted when considering the advantages of curved beam models compared to solid and plate ones. Modelling effort can be significantly reduced (solid models require cumbersome post- and pre-processing even in relatively simple cases) and parametric analyses can be facilitated (construction of multiple solid models is quite cumbersome, especially for complex geometries). Moreover, the FEM beam elements commonly employed in commercial software can be inaccurate. Thus, the main purpose of this beam formulation is to remain simple and with the least number of degrees of freedom needed to describe its behavior accurately. The main conclusions that can be drawn from this investigation are: a. Highly accurate results can be obtained by the use of b-splines in the AEM technique as well as NURBS in Finite Element beam formulations for the static analysis of the proposed curved beam elements. Discretization elements are significantly reduced by the use of b-splines in the AEM and post-processing of the results becomes much easier with NURBS in FEM. Knot insertion has proven to be very efficient when the number of unknowns is increased.

0 (x) 1144.57 (y) 0(x) 1072 (y) –

62,817 (x) 1571 (y) 952.89 (x) 1403 (y)

3408

6198

49

b. In general, open shaped cross sections suffer mainly from warping while closed ones from distortion. The consideration of up to secondary higher order phenomena (Eqs. (9)) is generally accurate. However, in cases of very thin-walled cross sections either open or close shaped, tertiary phenomena might need to be considered. As a rule of thumb in choosing the least number of modes in order to achieve the maximum accuracy when the proposed method is employed, the limitations t=d > 0:05 and d=L > 0:05 can be applied. c. Increase in curvature causes increase in the distortion due to torsion for thin-walled box-shaped cross sections. Cross sections with cantilever plates at both sides undergo less severe distortional effects due to the fact that bending resisting mechanisms become of importance, too. d. Direction and position of loads can play a significant role in the behavior of curved beams with monosymmetric and, thus, asymmetric cross sections due to the fact that the development of higher order warping (in open sections) or distortional (in close shaped) phenomena can be significantly altered. e. Composite cross sections can be easily handled with the proposed BEM approach due to the fact that only boundary discretization is required. Additionally, changes in the arrangement of materials can be easily integrated in the model and significantly improve the resistance of the curved beams against warping and distortion. f. The ratios of distortional and warping Moments to the bending ones can be indicative of the curved beam’s behavior and offers an additional insight into the resisting mechanisms that dominate.

–

0

20,418 (x) 245 (y) 1691.09 (x) 262 (y)

11,338 3357.5 –

–

3130

161,642 26,916 9152 162,115 0 0 162,023 26,969 –

3791 3992 7721 8425

161,723 26,929 18,463

0.3914 0 0.0174 52,122 0.3899 0.0038 0.0169 51,825 0.4001 0.0039 0.0170 52,456

0.2616 0.0019 0.0116 21,151

R = 1 AEM 100 constant R = 100 FEM solid 6600

To summarize, this research effort is a contribution to the structural analysis of spatial curved beam elements of composite arbitrary cross section including distortional effects and integrating Isogeometric tools in the final 1D beam’s analysis. As a future research suggestion, the cross sectional analysis can be further extended to account for the occurrence of multiple cracks in the cross section’s domain by employing meshfree discretization techniques in terms of reproducing kernels (RK) [74–76], extended FEM (XFEM) [77], extended IGA (XIGA) with NURBS [60,78] or Bézier extraction based T-spline extended IGA (BEBT-XIGA) [79]. Acknowledgements

Pz lateral loading

References

2

2

M SDx;y ðkN m Þ at x = 0

S

M PDx;y ðkN m Þ at x = 0

M /S ðkN m Þ at x = 0

2

CY

2

M /P ðkN m Þ at x = 0

S

2

M /P ðkN m Þ at x = 0

M Y ðkN mÞ at x = 0 M t ðkN mÞ at x = 0

kN rmax ðm 2 Þ at x = 0 xx kN smax xy ðm2 Þ at x = 0

w ðmÞ at x = L hx ðradÞ at x = L hY ðradÞ at x = L

This work has been supported by IKY Fellowships of Excellence for Postgraduate Studies in Greece-Siemens Program.

;

R = 100 AEM 10 cubic b-splines

Table 7 Kinematical components, Shear stresses and Stress resultants of the Box-shaped (Fig. 10c and d) beams examined in example 2.

R = 100 AEM 10 cubic b-splines composite

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

[1] Sapountzakis EJ, Tsiptsis IN. Generalized warping analysis of curved beams by BEM. Eng Struct 2015;100:535–49. [2] Love AEH. A treatise on the mathematical theory of elasticity. New York: Dover; 1944. [3] Vlasov V. Thin walled elastic beams. 2nd ed. Washington, DC: National Science Foundation; 1961. [4] Kollbrunner CF, Basler K. Torsion in structures – an engineering approach. New York: Springer; 1969. [5] Heilig R. A contribution to the theory of box girders of arbitrary cross-sectional shape (Translated from the German by C.V. Amerongen). Publication No. 145. Cement and Concrete Association; 1971. [6] Kristek V. Tapered box girders of deformable cross section. J Struct Div, ASCE 1970;96(8):1761–93. [7] Wright RN, Abdel-Samad SR, Robinson AR. BEF analogy for analysis of box girders. J Struct Div, ASCE 1968;94(7):1719–43. [8] Steinle A. Torsion and cross-sectional distortion of the single cell box beam. Beton Stalbetonbau 1970;65:215–22. [9] Kollbrunner C, Hajdin N. Dünnwandige Stäbe 2, Stäbe mit deformierbaren Querschnitten, Nicht-elastisches Verhalten dünnwandiger Stäbe. SpringerVerlag; 1975.

50

I.N. Tsiptsis, E.J. Sapountzakis / Computers and Structures 191 (2017) 33–50

[10] Kermani B, Waldron P. Analysis of continuous box girder bridges including the effects of distortion. Comput Struct 1993;47(3):427–40. [11] Kim YY, Kim JH. Thin-walled closed box beam element for static and dynamic analysis. Int J Numer Meth Eng 1999;45:473–90. [12] Park NH, Lim NH, Kang YJ. A consideration on intermediate diaphragm spacing in steel box girder bridges with a doubly symmetric section. Eng Struct 2003;25:1665–74. [13] Park NH, Choi S, Kang YJ. Exact distortional behavior and practical distortional analysis of multicell box girders using an expanded method. Comput Struct 2005;83:1607–26. [14] Razaqpur AG, Li HG. A finite element with exact shape functions for shear lag analysis in multi-cell box-girders. Comput Struct 1991;39(1–2):155–63. [15] Osadebe NN, Chidolue CA. Torsional-distortional response of thin-walled mono symmetric box girder structures. Int J Eng Res Appl 2012;2(3):814–21. [16] Schardt R. Verallgemeinerte Technische Biegetheorie. Springer-Verlag; 1989. [17] Schardt R. Generalized Beam Theory – an adequate method for coupled stability problems. Thin-Wall Struct 1994;19:161–80. [18] Ranzi G, Luongo A. A new approach for thin-walled member analysis in the framework of GBT. Thin-Wall Struct 2011;49:1404–14. [19] Jönsson J. Distortional theory of thin-walled beams. Thin-Wall Struct 1999;33 (4):269–303. [20] Jönsson J, Andreassen MJ. Distortional eigenmodes and homogeneous solutions for semi-discretized thin-walled beams. Thin-Wall Struct 2011;49:691–707. [21] Andreassen MJ, Jönsson J. Distortional solutions for loaded semi-discretized thin-walled beams. Thin-Wall Struct 2012;50:116–27. [22] Andreassen MJ, Jönsson J. Distortional buckling modes of semi-discretized thin-walled columns. Thin-Wall Struct 2012;51:53–63. [23] Andreassen MJ, Jönsson J. A distortional semi-discretized thin-walled beam element. Thin-Wall Struct 2013;62:142–57. [24] Silvestre N, Camotim D. On the mechanics of distortion in thin-walled open sections. Thin-Wall Struct 2010;48:469–81. [25] Camotim D, Dinis PB. Coupled instabilities with distortional buckling in coldformed steel lipped channel columns. Thin-Wall Struct 2011;49:562–75. [26] Dinis PB, Camotim D. Post-buckling behaviour and strength of cold-formed steel lipped channel columns experiencing distortional/global interaction. Comput Struct 2011;89:422–34. [27] El Fatmi R, Ghazouani N. Higher order composite beam theory built on SaintVenant’s solution. Part-I: theoretical developments. Compos Struct 2011;93:557–66. [28] Ferradi MK, Cespedes X. A new beam element with transversal and warping eigenmodes. Comput Struct 2014;131:12–33. [29] Genoese A, Genoese A, Bilotta A, Garcea G. A generalized model for heterogeneous and anisotropic beams including section distortions. ThinWall Struct 2014;74:85–103. [30] Dikaros IC, Sapountzakis EJ. Distorsional analysis of beams of arbitrary cross section by BEM. In: Proc of the 11th HSTAM international congress on mechanics, Athens, Greece, 27–30 May; 2016. [31] Dabrowski R. Curved thin-walled girders theory and analysis. Cement and Concrete Association; 1968. [32] Bazant ZP, El-Nimeiri M. Stiffness method for curved box girders at initial stress. J Struct Div, ASCE 1974;100(10):2071–90. [33] Oleinik JC, Heins CP. Diaphragms for curved box beam bridges. J Struct Div, ASCE 1975;101(10):2161–78. [34] Heins CP, Oleinik JC. Curved box beam bridge analysis. Comput Struct 1976;6:65–73. [35] Sakai F, Nagai M. A proposal for intermediate diaphragm design in curved steel box girder bridges. Proc Jpn Soc Civ Eng 1981;305:11–22. [36] Nakai H, Murayama Y. Distortional stress analysis and design aid for horizontally curved box girder bridges with diaphragms. Proc Jpn Soc Civ Eng 1981;309:25–39. [37] Martin TL, Heins CP. Analysis of distortional stress in curved I girder bridges. Comput Struct 1978;8:691–701. [38] Zhang SH, Lyons LPR. A thin-walled box beam finite element for curved bridge analysis. Comput Struct 1984;18(6):1035–46. [39] Zhang SH, Lyons LPR. The application of the thin-walled box beam element to multibox bridge analysis. Comput Struct 1984;18(5):795–802. [40] Nakai H, Yoo H. Analysis and design of curved steel bridges. New York: McGraw-Hill; 1988. [41] Yabuki T, Arizumi Y. A provision on intermediate diaphragm spacing in curved steel-plated box-bridge-girders. Struct Eng/Earthq Eng, JSCE 1989;6 (2):207–16. [42] Razaqpur AG, Li HG. Refined analysis of curved thin-walled multicell box girders. Comput Struct 1994;53(1):131–42. [43] Petrov E, Géradin M. Finite element theory for curved and twisted beams based on exact solutions for three-dimensional solids, part 1: beam concept and geometrically exact nonlinear formulation. Comput Methods Appl Mech Eng 1998;165:43–92. [44] Kim YY, Kim Y. A one-dimensional theory of thin-walled curved rectangular box beams under torsion and out-of-plane bending. Int J Numer Meth Eng 2002;53:1675–93. [45] Park NH, Choi YJ, Kang YJ. Spacing of intermediate diaphragms in horizontally curved steel box girder bridges. Finite Elem Anal Des 2005;41(9–10):925–43. [46] Kang YJ, Yoo CH. Thin-walled curved beams, I: formulation of nonlinear equations. J Eng Mech, ASCE 1994;120(10):2072–101.

[47] Huang DZ. Dynamic analysis of steel curved box girder bridges. J Bridge Eng, ASCE 2001;6:506–13. [48] Yang YB, Wu CM, Yau JD. Dynamic response of a horizontally curved beam subjected to vertical and horizontal moving loads. J Sound Vib 2001;242:519–37. [49] Nallasivam K, Dutta A, Talukdar S. Dynamic analysis of horizontally curved thin-walled box-girder bridge due to moving vehicle. Shock Vibrat 2007;14:229–48. [50] Zhang Y, Hou Z, Li Y, Wang Y. Torsional behaviour of curved composite beams in construction stage and diaphragm effects. J Constr Steel Res 2015;108:1–10. [51] Yoo CH, Kang J, Kim K. Stresses due to distortion on horizontally curved tubgirders. Eng Struct 2015;87:70–85. [52] Arici M, Granata MF. Unified theory for analysis of curved thin-walled girders with open and closed cross section through HSA method. Eng Struct 2016;113:299–314. [53] American Association of State Highway and Transportation Officials. Guide specifications for horizontally curved highway bridges, Washington, DC; 1993. [54] Hanshin Expressway Public Corporation. Guidelines for the design of horizontally curved girder bridges, Japan; 1988. [55] Dikaros IC, Sapountzakis EJ. Generalized warping analysis of composite beams of arbitrary cross section by BEM. Part I: theoretical considerations and numerical implementation. J Eng Mech, ASCE 2014;140(9):04014062. [56] Katsikadelis JT. Boundary elements: theory and applications. AmsterdamLondon: Elsevier; 2002. [57] Hughes T, Cottrell J, Bazilevs Y. Isogeometric analysis: toward Integration of CAD and FEA. Wiley; 2009. [58] Katsikadelis JT. The analog equation method. A boundary – only integral equation method for nonlinear static and dynamic problems in general bodies. Theoret Appl Mech 2002;27:13–38. [59] Valizadeh N, Bui TQ, Vu VT, Thai HT, Nguyen MN. Isogeometric simulation for buckling, free and forced vibration of orthotropic plates. Int J Appl Mech 2013;5(2):1350017. [60] Bui TQ. Extended isogeometric dynamic and static fracture analysis for cracks in piezoelectric materials using NURBS. Comput Methods Appl Mech Eng 2015;295:470–509. [61] Yin S, Yu T, Bui TQ, Liu P, Hirose S. Buckling and vibration extended isogeometric analysis of imperfect graded Reissner-Mindlin plates with internal defects using NURBS and level sets. Comput Struct 2016;177:23–38. [62] Liu S, Yu T, Bui TQ, Yin S, Thai D-K, Tanaka S. Analysis of functionally graded plates by a simple locking-free quasi-3D hyperbolic plate isogeometric method. Composites B 2017;120:182–96. [63] Lai W, Yua T, Bui TQ, Wanga Z, Curiel-Sosa JL, Das R, et al. 3-D elasto-plastic large deformations: IGA simulation by Bézier extraction of NURBS. Adv Eng Softw 2017;168:68–82. [64] Sapountzakis EJ, Tsiptsis IN. Quadratic B-splines in the analog equation method for the nonuniform torsional problem of bars. Acta Mech 2014;225:3511–34. [65] Tsiptsis IN, Sapountzakis EJ. Isogeometric analysis for the dynamic problem of curved structures including warping effects. Mech Based Des Struct Mach 2016. http://dx.doi.org/10.1080/15397734.2016.1275974. [66] Sapountzakis EJ, Tsiptsis IN. B-splines in the analog equation method for the generalized beam analysis including warping effects. Comput Struct 2017;80:60–73. [67] Xu XS, Zhong WX, Zhang HW. The Saint-Venant problem and principle in elasticity. Int J Solids Struct 1997;34(22):2815–27. [68] Reagan SW, Pilkey WD. Constrained torsion of prismatic bars. Finite Elem Anal Des 2002;38:909–19. [69] Beer G, Smith I, Duenser Ch. The boundary element method with programming – for engineers and scientists. Wien, New York: Springer; 2008. [70] Piegel L, Tiller W. The NURBS book. Berlin, Heidelberg: Springer; 1997. [71] Oñate E. Structural analysis with the finite element method. Linear statics. Basis and solids, vol. 1. CIMNE, Springer; 2009. [72] Davis PJ, Rabinowitz P. Methods of numerical integration. New York: Academic Press; 1975. [73] FEMAP for Windows. Finite element modeling and post-processing software. Help System Index Version 2010;11:1. [74] Tanaka S, Suzuki H, Sadamoto S, Okazawa S, Yu TT, Bui TQ. Accurate evaluation of mixed-mode intensity factors of cracked shear-deformable plates by an enriched meshfree Galerkin formulation. Arch Appl Mech 2017;87:279–98. [75] Sadamoto S, Ozdemir M, Tanaka S, Taniguchi K, Yu TT, Bui TQ. An effective meshfree reproducing kernel method for buckling analysis of cylindrical shells with and without cutouts. Comput Mech 2017;59:919–32. [76] Tanaka S, Suzuki H, Sadamoto S, Sannomaru S, Yu T, Bui TQ. J-integral evaluation for 2D mixed-mode crack problems employing a meshfree stabilized conforming nodal integration method. Comput Mech 2016;58:185–98. [77] Bui TQ, Zhang CH. Analysis of generalized dynamic intensity factors of cracked magnetoelectroelastic solids by X-FEM. Finite Elem Anal Des 2013;69:19–36. [78] Bui TQ, Hirose S, Zhang C, Rabczuk T, Wu CT, Saitoh T, et al. Extended isogeometric analysis for dynamic fracture in multiphase piezoelectric/ piezomagnetic composites. Mech Mater 2016;97:135–63. [79] Singh SK, Singh IV, Mishra BK, Bhardwaj G, Bui TQ. A simple, efficient and accurate Bézier extraction based T-spline XIGA for crack simulations. Theoret Appl Fract Mech 2017;88:74–96.