Accepted Manuscript Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling N. Nguyen-Thanh, K. Zhou, X. Zhuang, P. Areias, H. Nguyen-Xuan, Y. Bazilevs, T. Rabczuk PII: DOI: Reference: S0045-7825(16)31740-6 http://dx.doi.org/10.1016/j.cma.2016.12.002 CMA 11248 To appear in: Comput. Methods Appl. Mech. Engrg. Please cite this article as: N. Nguyen-Thanh, K. Zhou, X. Zhuang, P. Areias, H. Nguyen-Xuan, Y. Bazilevs, T. Rabczuk, Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling, Comput. Methods Appl. Mech. Engrg. (2016), http://dx.doi.org/10.1016/j.cma.2016.12.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. Highlights (for review) Highlights for Review 1. We present a multi-patch isogeometric large deformation thin shell formulation based on RTH splines. It is an extension of our previous work on RHT-spline shells to large deformations and multiple patches. The coupling is based on Nitsche's method and allows also coupling of a shell to a solid model. 2. Furthermore, we present a stress recovery technique to drive the adaptive hrefinement procedure in isogeometric thin structures. 3. The method is validated for several linear and non-linear benchmark problems including the pinched cylinder and hemispherical shell, a wind turbine rotor accounting for large deformations, a hemispherical shell with a stiffener and a pinched cylinder considering large deformations. *Manuscript Click here to download Manuscript: RHT_MultiplePatches_thinshell_Revised_R4.pdf 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Click here to view linked Referenc Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling N. Nguyen-Thanhc , K. Zhouc , X. Zhuangd , P. Areiase , H. Nguyen-Xuanf , Y. Bazilevsg , T. Rabczuka,b,h,∗ a Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh, Vietnam. of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Vietnam. c School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Ave, Singapore. d Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China. e Physics Department, University of Evora, Portugal. f Duy Tan University, Danang, Vietnam. g Department of Structural Engineering, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0085, USA. h Institute of Structural Mechanics, Bauhaus-University Weimar, Marienstrasse 15, 99423 Weimar, Germany. b Faculty Abstract We present an isogeometric thin shell formulation for multi-patches based on rational splines over hierarchical T-meshes (RHT-splines). Nitsche’s method is employed to efficiently couple the patches. The RHT-splines have the advantages of allowing a computationally feasible local refinement, are free from linear dependence, possess high-order continuity and satisfy the partition of unity and non-negativity. In addition, the C1 continuity of the RHT-splines avoids the rotational degrees of freedom. The good performance of the present method is demonstrated by a number of numerical examples. Keywords: Isogeometric Analysis, NURBS, PHT-splines, Thin shell, Multiple patches, Large deformation 1. Introduction Shell elements are widely used to model thin engineering structures such as car bodies in the automotive industry and lightweight components in aerospace structures. Shell elements that are used to model thin structures numerically are classified according to their thickness and the curvature of the mid-surfaces [1]. Based on the thickness of the elements, shell elements can be categorized into thin-shell elements [2, 3, 4] and thick-shell elements [5, 6, 7, 8]. Thin-shell elements are based on the Kirchhoff-Love theory [9] in which the transverse shear stresses are neglected. Therefore, thin shells require C1 displacement continuity, which is difficult to be achieved for free-form geometries with Lagrange finite element basis functions. The thin-shell theory requires the approximation of the deformation to have second- order square integrable derivatives. Using higher-order continuous discretization spaces in the context of thin shell analysis based on the Kirchhoff-Love ∗ Corresponding authors. Email address: thanhnhon@ntu.edu.sg (Nhon Nguyen-Thanh); kzhou@ntu.edu.sg (Zhou Kun); timon.rabczuk@tdt.edu.vn (Timon Rabczuk) 1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 theory obviates the need for rotational degrees of freedom (DOFs) or the additional discretization of the director field. Isogeometric analysis (IGA) was first introduced by Hughes et.al in [10] in order to more closely couple Computer Aided Design (CAD) with Computer Aided Engineering (CAE). The basic idea is to use CAD basis functions for numerical analysis. While the finite element method (FEM) is mostly used in CAE, common CAD basis functions are Non Uniform Rational B-splines (NURBS). Therefore, most studies in the context of IGA focus on NURBS-based formulations [11, 12, 13, 14, 15, 16, 17, 18]. Notwithstanding, NURBS also have certain shortcomings for numerical analysis. Due to the tensor-product form of NURBS, their control points are required to belong to a structured grid (e.g, in a rectangle in two dimensions). This causes an excessive overhead of control points with increasing refinement. Cottrell et al. [19] introduced a local refinement strategy; however, constraint equations are required, increasing the complexity and effort in implementation. Moreover, refinement still propagates across a given patch. Another disadvantage of NURBS is the C0 continuity across the patch boundaries. In that context, when two NURBS surfaces do not share a common boundary, even C0 continuity is not retained without perturbing at least one surface or introducing alternative techniques. IGA based on NURBS for plate/shells approaches has been presented [20, 21, 22, 23, 24, 25, 26, 27]. A formulation that employs solely the mid-surface position and fulfills the Kirchhoff constraint was first proposed for meshfree methods [28, 29] by using a higher continuous formulation. With the goal of overcoming the shortcomings of NURBS for IGA, in [30] it has been proposed to employ T-splines that allow local refinement. However, although T-splines allow for local adaptive refinement, the implementation of knot insertion under adaptive refinement is more intricate, specially for three-dimensional (3D) problems. Furthermore, Buffa et al. [31] showed that linear dependence of the basis functions may occur for generic T-meshes. The linear dependence issue was fixed in [32] where analysis suitable T-splines were defined. Approaches based on Polycube splines [33], Locally Refined (LR) splines [34] or Polynomial splines over hierarchical T-meshes (PHT-splines) [35, 36, 37] emerged as promising alternatives. The PHT-splines inherit the desirable properties of both NURBS and T-splines, and are capable of joining geometric objects without gaps. They also facilitate local refinement due to the inherent hierarchical structure. Moreover, the PHT-splines fulfill the linear independence requirement of the basis functions. The piecewise polynomial approach allows relatively simple adaptive refinement strategies and integration. Local refinement algorithms are relatively simple while the intricacy of knot insertion with T-splines might be high, particularly in 3D. Since the basis of the PHT-splines is polynomial instead of rational, conic sections used in engineering design cannot be exactly represented within the geometry in the framework of IGA. Therefore, a rational form of the PHTsplines needs to be devised. Recently contributions on PHT-splines to analysis can be found in [38, 39, 40, 41, 42, 43]. In this paper, we exploit a novel basis based on the rational form of polynomial splines over hierarchical T-meshes (RHT-splines) for multi-patch thin shell structures which are coupled efficiently by the Nitsche method [44]. The Nitsche method was first proposed to couple non-matching discretizations by weakly enforcing interface constraints. The method is based on replacing the Lagrange multipliers arising in a dual formulation with their representation using flux terms along the coupling interfaces. Recently, several methods have been explored in the IGA framework for multi-patch coupling [45, 46, 47, 48, 49]. In particular, in [47] the authors made use of the Nitsche approach for thin-shells, which we adopt in the present paper. The RHT-splines possess the key 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 features of PHT-splines and can be used to retain the exact geometry with conic sections. Moreover, the RHT-splines do not require the set-up of an additional mesh, nor additional nodes, and local refinement can be readily implemented through the refinement of geometric models. It employs Kirchhoff-Love theory in pristine form avoiding the need to introduce rotational DOFs due to the C1 continuity of the RHT-splines. We also present a stress recovery technique in IGA to drive the adaptive h-refinement procedure. We employ the Super-convergent Patch Recovery (SPR) technique proposed by Zienkiewicz et al. [50]. The procedure is a least square fitting of the finite element solutions over a local patch of elements at pre-selected points, where the rate of convergence is higher than the global rate. The patch in SPR comprises elements that are assembled around a central corner node. In IGA the recovered stress components are considered at an imaginary surface. An imaginary surface is constructed by the same RHT-spline basis functions which are used for the approximation of the displacement field. This imaginary surface will be fitted to these optimal sampling points in a least square sense for the recovered stress components. The manuscript can be considered as an extension of our work in [41]. The key novelty is the multi-patch coupling based on Nitsche’s method, the extension to large deformations and the h-adaptive refinement driven by the SPR method. The content of the paper is outlined as follows. In section 2, a brief review of the Bézier and Bspline basis functions is given. RHT-splines basis functions are studied in section 3, followed by a short introduction of the thin shell formulation in section 4. Section 5 presents the Nitsche approach for multi-patch coupling. Stress recovery in RHT-splines is presented in section 6. Section 7 shows several numerical examples verifying the effectiveness of the present formulation. The paper is concluded with some discussions and remarks. 2. Bézier and B-splines basis functions 2.1. Bézier curve The Bézier parametric function of degree n is given by n C(ξ ) = ∑ Mi,p (ξ )Pi ; i=1 ξ ∈ [0, 1] , (1) where Pi with i = 1, ..., n are the control points, n is the number of basis functions, and Mi,p (ξ ) are the Bernstein polynomials of degree p which are written as Mi,p (ξ ) = p! ξ i (1 − ξ ) p−i ; i!(p − i)! (p = n − 1) . (2) By applying de Casteljau’s algorithm, new control points Pik are evaluated by the following form  k = 1, ..., n k k−1 k−1 Pi (ξ ) = (1 − ξ )Pi (ξ ) + ξ Pi+1 (ξ ) with (3) i = 0, ..., n − k where Pi0 (ξ ) = Pi (ξ ). 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 2.2. B-Splines The B-spline basis functions Ni,p (ξ ) of order p = 0 (piece-wise constant) are defined recursively on the corresponding knot vector as follows  1 if ξi ≤ ξ ≤ ξi+1 (4) Ni,0 (ξ ) = 0 otherwise For p > 1, we write Ni,p (ξ ) = ξi+p+1 − ξ ξ − ξi Ni+1,p−1 (ξ ) . Ni,p−1 (ξ ) + ξi+p − ξi ξi+p+1 − ξi+1 (5) The B-splines surfaces are defined by the tensor product of the basis functions with parameters Ξ and expressed as n S(ξ , η ) = ∑ Ni (ξ , η )Pi , (6) i=1 where Pi are the control points in a bidirectional control net and Ni (ξ , η ) are the B-spline basis functions. In case of objects including conic sections, such as circles and ellipses, NURBS are utilized to give an exact representation of conic sections. Therefore, NURBS have become a standard tool in CAD modeling. Similar to a B-Spline surface, a NURBS surface is defined by n S (ξ , η ) = ∑ Ri (ξ , η )Pi ; Ri = i=1 Ni (ξ , η ) wi n ∑i=1 Ni (ξ , η ) wi , (7) where wi are weights. 3. Rational Splines over Hierarchical T-mesh (RHT-splines) 3.1. The RHT-splines basis functions The initial RHT-splines curve is given by n C (ξ ) = ∑ Ni (ξ )Pi , (8) i=1 where Ni (ξ ) are the cubic B-spline basis functions and Pi are the control points. The construction of RHT-splines basis functions is initiated from the C1 continuity. The interior knots are repeated once   n−1 n−1 1 1 , , 1, 1, 1, 1 . (9) Ξ = 0, 0, 0, 0, , , ..., n n n n It is seen that each basis function at ξi , ηi has its support being [ξi−1 , ξi+1 ] × [ηi−1 , ηi+1 ]. These points are called the basis vertices. These four basis functions are incorporated with knot vectors [ξi−1 , ξi−1 , ξi , ξi , ξi+1 ]×[ξi−1 , ξi , ξi , ξi+1 , ξi+1 ] and [ηi−1 , ηi−1 , ηi , ηi , ηi+1 ]×[ηi−1 , ηi , ηi , ηi+1 , ηi+1 ], respectively (see fig. 1). The Bézier ordinates representation is used to modify the basis functions. The basis functions are represented by specifying their 16 Bézier ordinates in every cell within 4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 1: Four basis functions associated at ξi , ηi . the support of the basis functions. A set of new Bézier ordinates is then generated by applying de Casteljau’s algorithm. We assume that among all the cells at level k, the cells φik , i = 1, ...,Ck are subdivided. For each i, if the basis function Nik (ξ , η ) does not vanish in some cells of φik , then we subdivide Nik (ξ , η ) into these cells at level k + 1 according to equation (3). Then we reset all the associated Bézier ordinates associated with the new basis vertices to zero. Let the new knot be given by (ξī , ηī ), then the four new basis functions are defined through the knot vector [ξī−1 , ξī−1 , ξī , ξī , ξī+1 ]×[ξī−1 , ξī , ξī , ξī+1 , ξī+1 ] and [ηī−1 , ηī−1 , ηī , ηī , ηī+1 ]×[ηī−1 , ηī , ηī , ηī+1 , ηī+1 ], respectively. There are four additional new basis functions when a new basis vertex is inserted. Therefore, four new control points need to be calculated, while the previous control points stay fixed. The refined RHT-splines surface is given by n S (ξ , η ) = ∑ Pi Nik (ξ , η ) + i=1 4 (ξ , η ) , ∑ Pj N k+1 j (10) j=1 where S(ξ , η ) is the new geometry parameter; Nik (ξ , η ) are the basis functions from the previous level k; N k+1 j (ξ , η ) are the new basis functions at level k + 1. Fig. 2 shows the change of basis functions when a new knot is inserted in two-dimensions (2D). 3.2. RHT-splines space Let T be the T-mesh, and Ω ⊂ R2 be the region occupied by T . Then we write n o S(p, q, α , β , T ) = s(x, y) ∈ Cα,β (Ω)|s(x, y)|φ ∈ P pq for any φ ∈ T , (11) where the space P pq consists of all bi-degree (p, q) polynomials and the space Cα,β consists of all the continuously bivariate functions up to order α in the x-direction and order β in the y-direction. The dimension formula of the spline space S(p, q, α , β , T ) when p ≥ 2α + 1 and q ≥ 2β + 1 has already been proved in [51]. Here, we set p = q = 3 and α = β = 1. The evaluation of the dimension formula is reduced to the following form dim S(3, 3, 1, 1, T ) = 4(V b +V + ) , where V b and V + are boundary and interior crossing vertices, respectively. 5 (12) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 (a) (b) Figure 2: Modification of basis functions in 2D: (a) basis functions at level k, and (b) basis functions at level k + 1. 3.3. RHT-splines surface The RHT-splines surface at level 0 is given by n S (ξ , η ) = ∑ Ri (ξ , η )Pi Ri (ξ , η ) = ; i=1 Ni (ξ , η ) wi n ∑i=1 Ni (ξ , η ) wi , (13) where Ni (ξ , η ) are the bi-cubic NURBS basis functions with C1 -continuity, wi are the weights and Pi are the control points. The RHT-splines surface at level k is given by k nk S (ξ , η ) = ∑ Rki (ξ , η )Pik , (14) i=1 where Rki (ξ , η ) are the same as at level 0, but the basis functions of the RHT-spline Nik (ξ , η ) are over the hierarchical T-mesh at level k. The geometry at level k + 1 is given by Sk+1 (ξ , η ) = nk+1 ∑ Rk+1 (ξ , η )Pik+1 ; Rik+1 (ξ , η ) = i i=1 Nik+1 (ξ , η ) wik+1 nk+1 ∑ i=1 3.4. Calculation of the control points The linear operator is defined as  
∂ S (ξ , η ) L S (ξ , η ) = S (ξ , η ) , = S (ξ , η ) , Sξ η (ξ , η ) ∂S 6 . (15) Nik+1 (ξ , η ) wik+1 (16) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 The geometric information of the RHT-splines surface at (ξ0 , η0 ) with which four basis functions are associated with indices i1 , i2 , i3 and i4 , are expressed as 4 L S (ξ0 , η0 ) = ∑ L Si (ξ0 η0 ) Fi = B · F . (17) i=1 The control points of the ith basis vertex are computed as F = L S (ξ0 , η0 ) · B−1 . The matrix B is given by [35]  (1 − λ ) (1 − µ ) λ (1 − µ )  −α (1 − µ ) α (1 − µ ) B=  −β (1 − λ ) −β λ αβ −αβ where α = 1 ∆u1 +∆u2 , β= 1 ∆v1 +∆v2 , (18)  λ µ µ (1 − λ ) αµ −α µ   , β λ β (1 − λ )  αβ −αβ (19) λ = α ∆u1 , and β = ∆v1 . 4. Thin shell model 4.1. Shell kinematics The displacement of a thin shell can be described by the displacement of its mid-surface, which is a 2D manifold embedded in 3D space. The mapping of shell middle surface is parameterized using curvilinear coordinates ξ 1 , ξ 2 ∈ A ⊂ R2 . The displacement of the points of the shell midsurface is defined by


u ξ 1, ξ 2 = x ξ 1, ξ 2 − X ξ 1, ξ 2 , (20) where X is the position vector of a material point of the shell mid-surface in the reference configuration, and x is the position vector of the same point in the deformed geometry. The position vector of a material point in the current (or deformed) configuration is given by


ϕ ξ 1 , ξ 2 , ξ 3 = x ξ 1 , ξ 2 + ξ 3 g3 ξ 1 , ξ 2 , (21) where ξ 3 denotes the coordinate in thickness direction, −0.5t 6 ξ 3 6 0.5t, with t being the shell thickness. Covariant basis vectors in the current configurations are obtained from the derivatives of the position with respect to the curvilinear coordinates as gα = ∂ϕ ∂ϕ = ϕ ,α + ξ 3 · nα ; g3 = ∂ ξα ∂ξ3 (α = 1, 2) . (22) The normal director at the corresponding reference point on the mid-surface is given by g3 = g1 × g 2 =n. |g1 × g2 | 7 (23) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 The covariant metric coefficients of the surface are defined as gi j = gi · g j ; Gi j = Gi · G j (i = 1, 2, 3) . (24) and the contravariant basis vectors are computed by  −1 . gi j = gi j gi = gi j gi ; The Green-Lagrange strain tensor is given by E=
1 T F F−I , 2 (25) (26) where I is the identity tensor and the deformation gradient is given by F = gi ⊗ Gi ; FT = Gi ⊗ gi . (27) The Green-Lagrange strain tensor is represented as E = Ei j Gi ⊗ G j . (28) The contravariant reference basis vectors is given by Ei j =
1 gi j − G i j . 2 (29) The Green-Lagrange strain is decomposed into a constant part due to membrane action and a linear part due to bending. The in-plane strain coefficients are given by E = Eαβ Gα ⊗ Gβ (α , β = {1, 2}) . (30) The covariant strain coefficients are given by Eαβ = εαβ + ξ 3 καβ . (31) The membrane strains εαβ are given by εαβ =
1 gα · gβ − Gα · Gβ , 2 (32) and the bending strains καβ , measuring the changes in curvature due to bending, are καβ = g3 · gα,β − G3 · Gα,β . (33) The substitution of equations (32) and equation (33) into equation (20) results in εαβ =
1 Gα · u,β − Gα · u,α , 2 8 (34) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 and

 1  καβ = − u,αβ + p u,1 · Gα,β × G2 + u,2 · G1 × Gα,β j¯ G3 · Gα,β + p [u,1 · (G2 × G3 ) + u,2 · (G3 × G1 )] j¯ (35) where j¯ = |G1 × G2 | . 4.2. Governing equations and discretization Applying the principle of virtual work, we obtain the stationarity of the total potential: δ U = δ Uint + δ Uext = 0 (36) S KL δ Uint = δ Uint + δ Uint , (37) in which S is the internal virtual energy for the solid and δ U KL is internal virtual energy for the where δ Uint int Kirchhoff-Love shell. The internal virtual energy for solid is given by S δ Uint = Z Ω S : δ EdΩ , (38) where S denotes the second Piola-Kirchhoff stress tensor, and δ E denotes the variation of strains. The internal virtual energy for the Kirchhoff-Love shell can be expressed as [52] KL δ Uint = Z A Z t/2 −t/2   σ : δ x,α ⊗ gα + δ n ⊗ g3 + ξ 3 δ n,α ⊗ gα det (∇Φ) d ξ 3 dA , (39) where the prefix δ identifies the test function, and σ is the Cauchy stress tensor. The external virtual work is given by δ Uext = Z A b · δ udA + Z Γ t̄ · δ udΓ , (40) where δ u denotes the variation of displacements and strains, b is the body force and t̄ is the prescribed traction. Integrating the Cauchy stress tensor σ over the thickness yields 1 n = ¯ j α 1 mα = ¯ j Z t/2 σ gα det (∇ϕ ) d ξ 3 (41) ξ 3 σ gα det (∇ϕ ) d ξ 3 , (42) −t/2 Z t/2 −t/2 where nα denotes the resultant membrane stress vector and mα the result bending stress vector. 9 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Substituting equations (41) and equation (42) into equation (39) results in Z (nα · δ x,α + mα · δ n̄,α ) dA ZA   αβ αβ αβ = n x,β · δ u,α + m n,β · δ u,α + m x,β · δ nα dA ZA   αβ αβ = n · δ εαβ + m · δ καβ dA . KL δ Uint = (43) A The resultant effective membrane stress and bending stress can be expressed by [53, 54] t 3 αβ γδ κγδ . D 12 (44)    m11 κ11 3 t m̄ =  m22  = D  κ22  . 12 m12 2κ12 (45) nαβ = tDαβ γδ εγδ ; Using Voigt’s notation one can write     n11 ε11 n̄ =  n22  = tD  ε22  ; n12 2ε12 mαβ =  An isotropic material matrix is given by  0 1 υ E  . υ 1 0 D= 1 − υ2 0 0 (1 − υ )/2  An orthotropic material matrix is defined as  E  D= 1 (1−υ12 υ21 ) υ12 E2 (1−υ12 υ21 ) υ21 E1 (1−υ12 υ21 ) E2 (1−υ12 υ21 ) 0 0  0  0 , B12 (46) (47) where E1 and E2 are the Young’s moduli; υ12 and υ21 are the Poisson coefficients; B12 is the shear modulus. Correspondingly, equation (43) can be rewritten in a compact form as KL δ Uint = Z A (n̄ · δ ε + m̄ · δ κ ) dA . The residual force vector is given by   ∂ Uint ∂ Uext Rr = + = Frint + Frext , ∂ ur ∂ ur int where Fext r is the external force and Fr is the internal force:  Z  ∂ε ∂κ int + m̄ : dA . Fr = − n̄ : ∂ ur ∂ ur A 10 (48) (49) (50) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 The second derivatives of the virtual work for the stiffness matrix is  2  ∂ Uint ∂ 2Uext ext int Krs = − , + Krs + = Krs ∂ ur ∂ us ∂ ur ∂ us (51) where the stiffness matrix K ext is the derivative of the external loads. The internal Kirchhoff-Love shell stiffness matrix is given by  Z  ∂ n̄ ∂ ε ∂ 2ε ∂ m̄ ∂ κ ∂ 2κ KL : + n̄ : + : + m̄ : dA , Kint = ∂ us ∂ ur ∂ u r ∂ us ∂ u s ∂ ur ∂ ur ∂ u s A where dA is the area of the midsurface in the reference configuration. For the solid, the stiffness matrix is given by [55]  Z  ∂S ∂E s Kint = : dΩ , Ω ∂ us ∂ ur where dΩ = q (52) (53) gi j d ξ1 d ξ2 d ξ3 . (54) 5. The weak form of the Nitsche approach The governing equation (36) of the principle of virtual work coupled with the Nitsche approach is given as δ Uint + δ UNits + δ Uext = 0 (55) in which δ UNits = − + Z ∗ ZΓ Γ∗ δ {S · n} · {u} da − Z Γ∗ δ {u} · {S · n} da + Z Γ∗ τS δ {u} · {u} da (56) τN (n · δ {u}) ({u} · n)da where n is the normal vector along the coupling interface Γ∗ , S denotes the stress tensor referred to the global coordinate system, τS and τN are the stabilization parameters of the shear and normal components, and where u = u(1) − u(2) ; Sn = β S1 n(1) + (1 − β ) S(2) n(2) β = [0, 1] . (57) In this work we have chosen β = 0.5. The stabilization parameters components are given as τS = CS ν ; h τN = CN λ , h (58) where λ and ν are the Lamé constants; h is element size. The stability and accuracy of Nitsche’s method are essentially influenced by the stabilization parameters CS and CN . The analysis will be unstable if the parameters are chosen to be too small. Otherwise, the Nitsche’s method will be reduced to a penalty method that entails disadvantages if they are too large. These values are 11 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 chosen according to the order of the basis functions (for p=3, the values of CS = CN ≈ 10; for p=4, the values can be chosen to be around 15 [46]). 5.1. Multiple-patch shell coupling For the Kirchhoff-Love patch coupling, the first two consistency terms of the Nitsche’s extension in terms of bending moments and force stress resultants are given [56, 47] Z Z     α α γ α α γ cs N + bγ M · δ (uα ) da δ ΠNits = − δ N + bγ M · uα da − − − ZΓ ZΓ ∗ ∗ Γ∗ δ (Q + Mn,s ) · u3 da − δ Mt · Φn da − Z Γ∗ Γ∗ Z (Q + Mn,s ) · δ (u3 ) da Γ∗ (59) Mt · δ (Φn ) da , in which   N α + bαγ M γ = nβ α + 2bαγ mβ γ nβ ; Q = mαβ |α nβ ;   Mn,s = mαβ nα tβ ; ,s Mt = mαβ nα nβ , (60) where nαβ are the stress resultants; mαβ are bending moments; nα ,tα are the covariant components of the normal and tangent vector; Φn = Φ · n is the rotation along the normal direction of the boundary. The covariant displacement components are given by u = uα gα + u3 n ; uα = u · gα ; u3 = u · n . (61) The term bα represents the mixed components of surface second form, bαγ = gαβ bγβ (62) The constraints of mid-surface displacements and the normal vector with the stability terms are enforced by the following term δ ΠstNits = Z Γ∗ τSt δ u · uds + Z Γ∗ τS t3 δ Φ · Φds 12 t3 + τN t (n · δ u) (u · n) ds + τN (n · δ Φ) (Φ · n) ds . 12 Γ∗ Γ∗ Z Z (63) The terms in brackets are defined as  (1)  (2) + (1 − β ) N α + bαγ M γ N α + bαγ M γ = β N α + bαγ M γ Q + Mn,s = β (Q + Mn,s )(1) + (1 − β ) (Q + Mn,s )(2) (1) Mt = β Mt (2) + (1 − β ) Mt u = u(1) − u(2)     (2) (1) . Φ = n(1) − n(2) − N3 − N3 12 (64) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 where N3 is the unit normal vector of mid-surface in the reference configuration. An example of the coupling interface between the solid and Kirchhoff-Love is illustrated in fig. 3. Coupling the Kirchhoff-Love shell with the solid gives δ Πcs Nits =− Z Γ∗ solid δ (Sn) · uda − Z Γ∗ δ u · (Sn)solid da (65) Figure 3: The coupling patches of the Kirchhoff-Love shell and solid 5.2. Discretizations The additional stiffness from the Nitsche extensions can be rewritten as Nits Nits ST K Nits = −Krs − Ksr + Krs , (66) Nits and K Nits are the stiffness contributions of the two consistency terms; K ST is the where Krs sr rs stiffness component of the stabilization terms. From equation (59), the stiffness matrix for the Kirchhoff-Love shell coupling becomes   Z ∂ N α + bα M γ Z Z γ ∂ (Q + Mn,s ) u3 ∂ uα ∂ Mt ∂ Φn Nits Krs = · da + · da+ · da . (67) ∂ ur ∂ us ∂ ur ∂ us ∂ ur ∂ us Γ∗ Γ∗ Γ∗ The displacements of the Kirchhoff-Love shell are assumed to be linear through the shell thickness u = u + ξ3 (n − N3 ) 13 (68) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 The derivatives of the displacement given by u,r = u,r + ξ3 n,r . (69) The displacement interpolation component is written by u,r = RU . (70) 6. Stress recovery in isogeometric analysis One of the most widely used stress recovery methods is the SPR method developed for the FEM [50]. In the SPR method, the reduced integration Gauss points inside each element are exploited as super-convergent points as the gradient of the approximated field variables, i.e. stress for solid mechanics, at these points possessing higher-order accuracy as compared with other points of the element, e.g. nodes. The main idea of the SPR method is to find the parameters describing the stress field at the nodes, such that the smoothed and continuous stress field σ ∗ can be constructed as σ ∗ = N(ξ , η )σ T , (71) where N(ξ , η ) are the same shape functions as used for the displacement field and σ are the nodal parameters that are different from the primary stress results projected as nodes. The recovered stress fields are assumed to follow the Lagrangian polynomial expansion where the parameters of the polynomial are obtained by minimizing the difference between σ and the recovered stress obtained from the Gauss points belonging to the “patch” of a node [50]. The idea of the SPR method can be extended to the IGA, i.e. by approximating the stress field using the same RHTsplines shape function for displacement approximation as for the stress recovery. Here, we follow the stress recovery method originally proposed in [57] and later applied in [58]. It was proposed that there exists an imaginary “stress surface” that can be approximated by the same RHT-splines function and the only unknowns are the knot parameters. Note that in IGA, the geometry is parameterized by knots. Therefore, we describe n σα∗ = ∑ m ∑ Ri, j (ξ , η )(Ti, j )α i=1 j=1 (α = 1, 2, · · · , 4) , (72) where Ri, j (ξ , η ) denote the RHT-spline shape functions being the same for the displacement field, and Ti, j (ξ , η ) are the stress parameters at the knots to be determined. Note that we can solve for each component of the stress tensor σα∗ , i.e. α rotates between 1 to 4 and (Ti, j )α corresponds to σα∗ . With N = n × m, equation (72) becomes σα∗ = N ∑ Rl (ξ , η )(Tl )α = RT Tα , (73) l=1 where R = R(ξ , η ) is an M by one vector collecting all the components of the RHT-splines shape functions and Tα is an M by one vector collecting all the knot parameters. We now introduce a norm G(Tα ) that computes the square of the difference between the recovered stress σα∗ and 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 primary stress σ̄α as kξ G(Tα ) = kη ∑  2 ∗ σ σ̄ − ∑ α p,q α p,q , (74) p=1 q=1 where kξ and kη are the numbers of the Gauss points in each direction of the patch. The goal is to minimize G(Tα ) with respect to each component of Tα . Now we substitute equation (73) into equation (74) and denote N = kξ × kη , and then equation (73) becomes N ∑ G(Tα ) = kG =1 2  RTkG Tα − σ̄αkG , where RkG = R(ξkG , ηkG ). The stationary condition of F(Tα ) leads to   ∂G ∂G = 0 (for l = 1, 2, ..., N) , = ∂ Tα ∂ (Tl )α (75) (76) and ATα = B , where N A= ∑ kG =1 RTkG RkG (77) N B= ; ∑ kG =1 RTkG σ̄αkG . (78) Therefore, we obtain a system of N equations to solve for the N unknowns of Tα . Note that since the stress parameters are directly involved as unknowns, the procedure needs to be applied to each stress component repeatedly. Therefore, it is more expensive than the stress recovery in the FEM, where the shape functions for stress recovery only depend on node and element coordinates; for more details see [50]. In summary, there are two key differences here as compared to the SPR method, namely • The stress parameters in the FEM are recovered at the nodes while in the IGA the stress parameters are recovered at knots. • The patch in the FEM are elements connected to a node while in the IGA it is defined as the overlapping areas between knots. It should be pointed out that here the Gauss points are used as sampling points to solve for knot parameters. The question of whether the super-convergent points are Gauss points or not for IGA remains an open issue. Recovery-type error estimation procedures have already been used by other researchers for nonlinear problems in adaptive analysis. It was discussed in [59] that the SPR method can perform well, even if the sampling points used for stress recovery are not the super-convergent points. This is worth further investigation on the mathematical proof of finding optimal convergent points. 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 7. Numerical examples In this section, we show the performance of the proposed method through benchmark examples, and an industry problem, considering both linear and non-linear problems. These problems were often used to asses the element’s robustness and accuracy under complex strain states. The results are compared by available reference solutions. The polynomial order p = 3 is used for both RHTplines and NURBS. 7.1. Pinched Cylinder A cylindrical shell with rigid end diaphragm is subjected to a point load at the center of the cylindrical surface. The related parameters are: length of the cylinder L = 600mm; radius R = 300mm; thickness t = 3mm; Young’s modulus E = 3 × 106 N/mm2 ; Poisson’s ratio ν = 0.3; point load P = 1N. The expected deflection under a concentrated load is 1.8248×10−5 [60]. The geometry of cylinder is subdivided into four patches and the meshes are shown in fig. 4. The contour plots of displacements force by coupling eight patches and convergence of displacement at the loading point are shown in fig. 5. A good agreement with the reference solution can be observed. It is also observed that the RHT-splines exhibit the higher accuracy than the NURBS-based approach when the same order approximation is used. These results verify our implementation. Figure 4: Pinched cylinder with diaphragms boundary conditions and the mesh. 16 −5 x 10 1.8 Reference NURBS RHT−spline 1.6 1.4 Displacement 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 1.2 1 0.8 0.6 0.4 0.2 0 (a) 0 0.2 0.4 0.6 0.8 1 1.2 Number of DOFs 1.4 1.6 1.8 2 4 x 10 (b) Figure 5: (a) Contour plots of displacement, force and moment resultants on deformed configuration (scaling factor = 1×107 ). (b) Displacement convergence of the pinched cylinder. 17 7.2. Hemispherical shell Let us consider a hemispherical shell subjected to two opposite point loads F = 2N. The bottom circumferential edge of the hemisphere is free (fig. 6). The parameters used are Young’s modulus E = 6.825×107 N/m2 , Poisson’s ratio ν = 0.3, radius R = 10m and thickness of the shell t = 0.04m. The reference value of the radial displacement under the point loads is 0.0924 [60]. The contour of deformed configuration of displacement component by four patches and the convergence of displacement under the applied loads are shown fig. 7. It is again seen that RHT-splines produce more accurate solution than the NURBS. Figure 6: Hemispherical shell at point loads and the mesh. 0.1 0.09 Reference NURBS RHT−spline 0.08 0.07 Displacement 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0.06 0.05 0.04 0.03 0.02 0.01 0 0 5000 10000 15000 Number of DOFs (a) (b) Figure 7: (a) Contour plot of displacement and deformed configuration (scaling factor = 30) of the hemispherical shell. (b) Displacement convergence of the hemispherical shell. 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 7.3. Application to a wind turbine rotor We consider a wind turbine rotor with gravity loading (g = 9.81m/s2 ) as presented in fig. 8a. For this example, we assume isotropic material with the parameters: Young’s modulus E = 19 × 109 N/m2 , Poisson’s ratio ν = 0.29, and the thickness of the wind turbine rotor is shown in fig. 8b. The geometry of this problem is subdivided into sixteen patches and the mesh is shown in fig. 9. For this example, we use different parameters for the mesh per patch and these are refined from the tip to roof; fig. 10 shows the contour plot of the displacement and convergence of displacement under the gravity loading, while fig. 11 presents the force and moment resultants. The good performance of the method is confirmed by solving this complex geometry problem. Maximum displacement at the tip of the blade using the RHT-splines approach is u = 3.05m, while the reference value is ure f = 3.03m [61, 62] and the error is 0.66%. (a) (b) Figure 8: (a) The wind turbine rotor [48]. (b) Blade thickness distribution. 19 Figure 9: The wind turbine blade is subdivided into 16-patches with coarse mesh and refinement meshes. 3.2 3 Reference RHT−splines 2.8 Displacement 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 2.6 2.4 2.2 2 (a) 1 1.5 2 2.5 3 3.5 4 Number of DOFs 4.5 5 5.5 6 4 x 10 (b) Figure 10: (a) Contour plot displacement (scaling factor = 20) of the wind turbine blade. (b) Displacement convergence of the wind turbine blade. 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 (a) (b) (c) (d) (e) (f) Figure 11: Contour plot of force and moment resultants of the wind turbine blade: (a,b,c) normal force Nxx ,Nyy , Nxy , respectively; (d,e,f) bending moment Mxx ,Myy , Mxy , respectively. 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 7.4. Hemispherical shell with a stiffener Let us consider√ a hemispherical shell including a stiffener (see fig. 12). The parameters are: R1 = 10m, R2 = 5 3m, t1 = 0.1m, t2 = 0.4m, L2 = 0.4m, E = 6.825 × 107 kN/m2 , ν = 0.3 , ρ = 500kg/m3 , g = 10m/s2 , q = 100kN/m2 , α = 30o , β = 10o . The cylindrical stiffener with a square cross-section is simply supported at the bottom surface. The structure is subjected to selfweight and a pressure q is acting on the outer surface of the shell and stiffener. Due to symmetry only a quarter of the system has to be considered. The meshes are illustrated in fig. 13. Contour plots of the displacement of shell and solid coupling are shown in fig. 14. Again, the proposed method is in good agreement with the reference solution. Figure 12: Hemispherical shell with a stiffener. (a) (b) Figure 13: Hemispherical shell with a stiffener. (a) Coarse mesh. (b) Fine mesh 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 (a) (b) Figure 14: Vertical displacement contour plot on the deformed configuration. (a) Reference result [10]. (b) Present method. 7.5. Pinched cylinder with end diaphragms We reconsider the cylindrical shell with rigid end diaphragms in section 7.1, but now under large deformations. The length, radius and thickness of the cylinder are L = 200mm, R = 200mm, t = 1mm, respectively. The material parameters are: Young’s modulus E = 30000N/mm2 and Poisson’s ratio ν = 0.3. The pinching is increased to P = 12000N. Large deformations and the load deflection curve are shown in fig. 15. The results obtained show good agreement with the reference results [63, 29]. 23 (a) (b) 12000 WA (Sze) 10000 UB (Sze) W (RHT−spline) A UB (RHT−spline) 8000 Pinching load P 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 6000 4000 2000 0 −20 (c) 0 20 40 60 Displacements at points A and B 80 100 (d) Figure 15: Deformed pinched cylinder with end diaphragms. (a) At t=10. (b) At t=35. (c) At t=50. (d) The displacement of points A and B. 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 8. Conclusions An isogeometric approach based on RHT-splines is formulated for the Kirchhoff-Love thinshell structures. Owing to the C1 -continuity of RHT-splines, only the mid-surface of the shell needs to be discretized, and the isogeometric approach based on RHT-splines fulfills the Kirchhoff-Love constraint. The proposed formulation requires only displacement degrees of freedom. The method is simple to implement, efficient in computation and can be easily applied to complex surface models by thin-shell patches with the use of Nitsche’s method. The advantages of Nitsche’s method are its variational consistency in comparison to other coupled techniques, and the positive definiteness, symmetry of the stiffness matrix, which can be achieved without additional degrees of freedom. Various benchmark problems are tested using the present RHT-splines thin-shell formulation. The results obtained achieve high reliability in all cases. We also demonstrate the performance of the proposed method for a wind turbine blade subjected to a gravity or aerodynamic load, and an example shell-to-solid coupling. The stress recovery approach in IGA to drive the adaptive h-refinement procedure is also presented. An imaginary surface is constructed by the same RHT-spline basis functions which are used for approximation of the displacement field. The proposed method shows good agreement with reference solutions for linear and non-linear problems. 9. Acknowledgment The first author gratefully acknowledges the financial support of the Singapore Maritime Institute (Grant SMI-2014-MA11). Yuri Bazilevs was partially supported by the NASA (Grant NNX15AW13A). Appendix A. The relation between strains and stresses The derived from the strain energy Uint S= ∂ Uint ∂E (79) where S is the second Piola- Kirchhoff stress tensor, E is the Green-Lagrange strain tensor. Stress and strain tensor are related material tensor D= ∂ S ∂ 2Uint = ∂E ∂ E2 (80) A linear relation between strains and stresses are given by using St.Venant-Kirchhoff material model S=D:E (81) S = Si j Gi ⊗ G j (82) Si j = Di jkl Ekl (83) where Related between the Cauchy and the second Piola-Kirchhoff stress tensor by the deformation gradient 1 σ= F · S · FT (84) (det F) 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 S = det F · F−1 · σ · F−T (85) Appendix B. Derivatives of the membrane and bending components Derivatives of the membrane component from equation (31) εαβ ,r =

∂ εαβ 1 1 1 = gαβ − Gαβ ,r = gαβ ,r = gα,r gβ + gα gβ ,r ∂ ur 2 2 2 gα,r = r,α,r = (R + u),α,r = ∑ Ri,α Ui,r (86) (87) i where Ui,r are the displacements at the control points; Ri,α are RHT-splines basis functions. The derivatives of the bending component is given by

καβ ,r = Gαβ − gαβ ,r = −gαβ,r = − gα,β ,r · n + gα,β · n,r (88) where gα,β ,r = ∑ Ri,αβ Ui,r gα,β = ∑ Ri,αβ Ui i i (89) The derivatives of the normal vector is given by n,α = where ĝ3 =g1 × g2 ; ḡ3 = |g1 × g2 | ; ĝ3,α · ḡ3 − ĝ3 · ḡ3,α ḡ23 ĝ3,α = g1,α × g2 + g1 × g2,α ĝ3 · ĝ3,α ḡ3,α = ḡ3 (90) (91) References [1] M. Bischoff, K.-U. Bletzinger, W.A. Wall, and E. Ramm. Models and Finite Elements for Thin-Walled Structures. John Wiley and Sons, New York, 2004. [2] P.K. and T. Belytschko. Analysis of thin shells by the element free Galerkin method. International Journal of Solids and Structures, 33(20-22):3057–3080, 1996. [3] K.J. Bathe and E. Dvorkin. Our discrete Kirchhoff and isoparametric shell elements for nonlinear analysis - An assessment. Computers and Structures, 16:89–98, 1983. [4] D. Millan, A. Rosolen, and M. Arroyo. Thin shell analysis from scattered points with maximum-entropy approximants. International Journal for Numerical Methods in Engineering, 85(6):723–751, 2011. [5] D. Wang and J.S. Chen. Locking-free stabilized conforming nodal integration for meshfree Mindlin-Reissner plate formulation. Computer Methods in Applied Mechanics and Engineering, 193:1065–1083, 2004. 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [6] J.S. Chen and D. Wang. A constrained reproducing kernel particle formulation for shear deformable shell in cartesian coordinates. International Journal for Numerical Methods in Engineering, 68:151–172, 2006. [7] N. Nguyen-Thanh, T. Rabczuk, H.Nguyen-Xuan, and S. Bordas. An alternative alpha finite element method with stabilized discrete shear gap technique for analysis of Mindlin-Reissner plates. Finite Elements in Analysis and Design, 47(5):519–535, 2011. [8] N. Nguyen-Thanh, T. Rabczuk, H.Nguyen-Xuan, and S. Bordas. A smoothed finite element method for shell analysis. Computer Methods in Applied Mechanics and Engineering, 198(2):165–177, 2008. [9] S. Timoshenko and S. Woinowsky-Krieger. Theory of Plates and Shells (2nd edn). McGrawHill, New York, 1959. [10] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194:4135–4195, 2005. [11] Y. Bazilevs, V.M. Calo, T.J.R Hughes, and Y. Zhang. Isogeometric fluid-structure interaction: theory, algorithms, and computations. Computational Mechanics, 43:3–37, 2008. [12] Y. Bazilevs, M.-C.Hsu, I.Akkerman, S.Wright, K.Takizawa, B.Henicke, T.Spielman, and T.E.Tezduyar. 3D simulation of wind turbine rotors at full scale. Part I: Geometry modeling and aerodynamics. International Journal for Numerical Methods in Engineering, 65(13):207–235, 2011. [13] C.H.Thai, H.Nguyen-Xuan, N.Nguyen-Thanh, T.H. Le, T. Nguyen-Thoi, and T. Rabczuk. Static, free vibration, and buckling analysis of laminated composite Reissner-Mindlin plates using NURBS based isogeometric approach. International Journal for Numerical Methods in Engineering, 2012. [14] D. Schillinger, L. Dede, M.A. Scott, J.A. Evans, M.J. Borden, E. Rank, and T.J.R. Hughes. An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of nurbs, immersed boundary methods, and t-spline cad surfaces. Computer Methods in Applied Mechanics and Engineering, 2012. [15] F.Amiri, D.Millan, Y.Shen, T.Rabczuk, and M.Arroyo. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 69:102–109, 2014. [16] R.Kruse, N.Nguyen-Thanh, L.De Lorenzis, and T.J.R.Hughes. Isogeometric collocation for large deformation elasticity and frictional contact problems. Computer Methods in Applied Mechanics and Engineering, 296:73–112, 2015. [17] J. Kiendl, M. Ambati, L. De Lorenzis, H. Gomez, and A. Reali. Phase-field description of brittle fracture in plates and shells. Computer Methods in Applied Mechanics and Engineering, in press, 2016. 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [18] M. Ambati and L. De Lorenzis. Phase-field modeling of brittle and ductile fracture in shells with isogeometric NURBS-based solid-shell elements. Computer Methods in Applied Mechanics and Engineering, in press, 2016. [19] J.A. Cottrell, T.J.R Hughes, and A. Reali. Studies of refinement and continuity in geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 196:4160– 4183, 2007. [20] D.J. Benson, Y. Bazilevs, and M.C. Hsu T.J.R. Hughes. Isogeometric shell analysis: The Reissner-Mindlin shell. Computer Methods in Applied Mechanics and Engineering, 199:276–289, 2010. [21] J. Kiendl, K.-U. Bletzinger, J. Linhard, and R. Wüchner. Isogeometric shell analysis with Kirchhoff-Love elements. Computer Methods in Applied Mechanics and Engineering, 198:3902–3914, 2009. [22] L. Chen, N. Nguyen-Thanh, H. Nguyen-Xuan, T. Rabczuk, S. Bordas, and G. Limbert. Explicit finite deformation analysis of isogeometric membranes. Computer Methods in Applied Mechanics and Engineering, 277:104–130, 2014. [23] N. Nguyen-Thanh, N. Valizadeh, M.N. Nguyen, H. Nguyen-Xuan, X. Zhuang, P. Areias, G. Zi, Y. Bazilevs, L. De Lorenzis, and T. Rabczuk. An extended isogeometric thin shell analysis based on Kirchhoff-Love theory. Computer Methods in Applied Mechanics and Engineering, 284:265–291, 2015. [24] J. Kiendl, H. Ming-Chen, and A. Reali Michael C.H. Wu and. Isogeometric KirchhoffLove shell formulations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 291:280–303, 2015. [25] J.F. Caseiroa, R.A.F. Valenteand A. Reali, J. Kiendl, F. Auricchio, and R.J. Alves de Sousa. Assumed natural strain NURBS-based solid-shell element for the analysis of large deformation elasto-plastic thin-shell structures. Computer Methods in Applied Mechanics and Engineering, 284:861–880, 2015. [26] T.X. Duong, F. Roohbakhshan, and R.A. Sauer. A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries. Computer Methods in Applied Mechanics and Engineering, 2016. [27] J. Niiranen, J. Kiendl, A.H. Niemi, and A. Reali. Isogeometric analysis for sixth-order boundary value problems of gradient-elastic Kirchhoff plates. Computer Methods in Applied Mechanics and Engineering, 2016. [28] T. Rabczuk and P.M.A. Areias. A meshfree thin shell for arbitrary evolving cracks based on an external enrichment. Computer Methods in Applied Mechanics and Engineering, 16(2):115– 130, 2006. [29] T. Rabczuk, P.M.A. Areias, and T. Belytschko. A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, 72:524–548, 2007. 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [30] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J.A. Evans T.J.R Hughes, S. Lipton, and M.A. Scottand T.W. Sederberg. Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering, 199:229–263, 2010. [31] A. Buffa, D. Cho, and G. Sangalli. Linear independence of the T-spline blending functions associated with some particular T-meshes. Computer Methods in Applied Mechanics and Engineering, 199:437–1445, 2009. [32] X. Li and M. A. Scott. Analysis-suitable t-splines: Characterization, refineability, and approximation. Mathematical Models and Methods in Applied Sciences, 24:1141–1164, 2014. [33] H. Wang, Y. He, X. Li, X. Gu, and H. Qin. Polycube splines. Computer-Aided Design, 40:721–733, 2008. [34] T. Dokken and V. Skytt. Locally refined splines. IV European Congress on Computational Mechanics (ECCM), Paris, France, 2010. [35] J. Deng, F. Chen, X. Li, C. Hu, W. Tong, Z. Yang, and Y. Feng. Polynomial splines over hierarchical T-meshes. Graphical Models, 70:76–86, 2008. [36] X. Li, J. Deng, and F. Chen. Surface modeling with polynomial splines over hierarchical T-meshes. Visual Computer, 23:1027–1033, 2007. [37] X. Li, J. Deng, and F. Chen. Polynomial splines over general T-meshes. Visual Computer, 26:277–286, 2010. [38] N. Nguyen-Thanh, H. Nguyen-Xuan, S. Bordas, and T. Rabczuk. Isogeometric finite element analysis using polynomial splines over hierarchical T-meshes. Materials Science and Engineering, 1, 2010. [39] N. Nguyen-Thanh, H. Nguyen-Xuan, S. Bordas, and T. Rabczuk. Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids. Computer Methods in Applied Mechanics and Engineering, 200(21-22):1892–1908, 2011. [40] L. Tian, F. Chen, and Q. Dua. Adaptive finite element methods for elliptic equations over hierarchical T-meshes. Journal of Computational and Applied Mathematics, 236:878–891, 2011. [41] N. Nguyen-Thanh, J. Kiendl, H. Nguyen-Xuan, R. Wüchner, K.U. Bletzinger, Y. Bazilevs, and T. Rabczuk. Rotation free isogeometric thin shell analysis using PHT-splines. Computer Methods in Applied Mechanics and Engineering, 200(47-48):3410–3424, 2011. [42] P. Wang, J.Xua, J.Deng, and F.Chen. Adaptive isogeometric analysis using rational PHTsplines. Computer-Aided Design, 43:1438–1448, 2011. [43] N. Nguyen-Thanh, J. Muthu, X. Zhuang, and T. Rabczuk. An adaptive three-dimensional RHT-splines formulation in linear elasto-statics and elasto-dynamics. Computational Mechanics, 53:369–385, 2014. 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [44] J. Nitsche. Uber ein variationsprinzip zur losung von dirichlet-problemen bei verwendung von teilraumen, die keinen randbedingungen unterworfen sind. Abhandlungen aus demMathematischen Seminar der Universität Hamburg, 1971. [45] V.P. Nguyen, P. Kerfriden, M. Brino, S.P.A. Bordas, and E. Bonisoli. Nitsche’s method for two and three dimensional nurbs patch coupling. Computational Mechanics, 53:1163–1182, 2014. [46] M.Ruess, D.Schillinger, A.I.Ö zcan, and E.Rank. Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries. Computer Methods in Applied Mechanics and Engineering, 269:46–71, 2014. [47] Y. Guo and M. Ruess. Nitsches method for a coupling of isogeometric thin shells and blended shell structures. Computer Methods in Applied Mechanics and Engineering, 284:881–905, 2015. [48] J. Kiendl, Y. Bazilevs, M.-C. Hsu, R. Wüchner, and K.-U. Bletzinger. The bending strip method for isogeometric analysis of Kirchhoff-Love shell structures comprised of multiple patches. Computer Methods in Applied Mechanics and Engineering, 199:2403–2416, 2010. [49] A. Apostolatos, R. Schmidt, R. Wüchner, and K.U. Bletzinger. A Nitsche-type formulation and comparison of the most common domain decomposition methods in isogeometric analysis. International Journal for Numerical Methods in Engineering, 97:473–504, 2014. [50] O.C. Zienkiewicz and J.Z. Zhu. The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique. International Journal for Numerical Methods in Engineering, 33:1331–1364, 1992. [51] J. Deng, F. Chen, and Y. Feng. Dimensions of spline spaces over T-meshes. Journal of Computational and Applied Mathematics, 194:267–283, 2006. [52] T. Belytschko, W.K. Liu, and B. Moran. Nonlinear Finite Elements for Continua and Structures. Wiley, 2000. [53] J.C. Simo and D.D. Fox. On a stress resultant geometrically exact shell-model 1. formulation and optimal parametrization. Computer Methods in Applied Mechanics and Engineering, 72:267–304, 1989. [54] J.C. Simo, D.D. Fox, and M.S. Rifai. On a stress resultant geometrically exact shell-model 2. the linear-theory - computational aspects. Computer Methods in Applied Mechanics and Engineering, 73:53–92, 1989. [55] C.V. Verhoosel S. Hosseini, J.J.C. Remmers and R.D. Borst. An isogeometric solid-like shell element for nonlinear analysis. International Journal for Numerical Methods in Engineering, 95:238–256, 2013. [56] W. Koiter. On the mathematical foundation of shell theory. Actes Congres Inter. Math., 269:123–130, 1970. 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [57] I. Babuska, T. Strouboulis, C. S. Upadhyay, S. K. Gangaraj, and K. Copps. Validation of a posteriori error estimators by numerical approach. International Journal for Numerical Methods in Engineering, 37:1073–1123, 1994. [58] H. Behrooz, G. Ahmad, and T. Mehdi. An isogeometrical approach to error estimation and stress recovery. European Journal of Mechanics and Solids, 31:101–109, 2012. [59] M. Ainsworth and J. T. Oden. A Posteriori Error Estimation in Finite Element Analysis. John Wiley & Sons, NewYork, 2000. [60] T. Belytschko, H. Stolarski, W.K. Liu, N. Carpenter, and J.S.J. Ong. Stress projection for membrane and shear locking in shell finite elements. Computer Methods in Applied Mechanics and Engineering, 51:221–258, 1985. [61] Y. Bazilevs, M.-C. Hsu, I. Akkerman, S. Wright, K. Takizawa, B. Henicke, T. Spielman, and T. E. Tezduyar. 3D simulation of wind turbine rotors at full scale. Part I: Geometry modeling and aerodynamics. International Journal for Numerical Methods in Engineering, 65:207–235, 2011. [62] Y. Bazilevs, M.-C. Hsu, J. Kiendl, R. Wuchner, and K.-U. Bletzinger. 3D simulation of wind turbine rotors at full scale. Part II: Fluidstructure interaction modeling with composite blades. International Journal for Numerical Methods in Engineering, 65:236–253, 2011. [63] K.Y. Sze, X.H. Liu, and S.H. Lo. Popular benchmark problems for geometric nonlinear analysis of shells. Finite Elements in Analysis and Design, 40:1551–1569, 2004. 31