Isogeometric analysis of large thin shell structures based on weak coupling of substructures with unstructured T-splines patches

Advances in Engineering Software 135 (2019) 102692

Contents lists available at ScienceDirect

Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft

Isogeometric analysis of large thin shell structures based on weak coupling of substructures with unstructured T-splines patches

T

Liu Zhenyua, Cheng Jinb, , Yang Minglongb, Yuan Peia, Qiu Chana, Gao Weic, Tan Jianrongb ⁎

a

State Key Lab of CAD&CG, Zhejiang University, Hangzhou 310027, PR China State Key Lab of Fluid Power & Mechatronic Systems, Zhejiang University, Hangzhou 310027, PR China c School of Civil & Environmental Engineering, The University of New South Wales, Sydney, NSW 2052, Australia b

ARTICLE INFO

ABSTRACT

Keywords: Isogeometric analysis (IGA) Unstructured T-splines Kirchhoff–Love shell element Substructure Dimensionality reduction method Weak coupling

In this paper, a new isogeometric analysis (IGA) method is proposed for large thin shell structures based on the weak coupling of substructures with unstructured T-splines patches. Firstly, the Kirchhoff–Love shell element is developed based on unstructured T-splines for discretizing the analysis models of large thin shell structures. Subsequently, the substructuring method is adopted to divide the analysis model of a large thin shell structure into several substructures composed of multiple unstructured T-splines patches. The approach for the weak coupling of substructures is introduced in detail. To improve the IGA efficiency of large thin shell structures, a dimensionality reduction method is put forward to reduce the stiffness matrix of an IGA model by firstly condensing the interior control points of its substructures’ matrices and then integrating the condensed substructures into a whole system. The effectiveness and applicability of the proposed IGA method with substructuring are verified by two examples.

1. Introduction Traditional finite element analysis (FEA) is usually time-consuming in the meshing procedure of transforming computer aided design (CAD) models into computer aided engineering (CAE) ones. Isogeometric analysis (IGA) [1] is a method for avoiding this tedious transformation by utilizing the same basis for analysis and modeling. Up to now, IGA method has been applied in many fields such as dynamic response analysis [2–5], computational fluid dynamics [6], shape optimization [7–9], topology optimization [10–12], heat transfer analysis [13,14], and so on. Smooth basis functions enable IGA to exactly discretize analysis models, which makes it suitable for the analyses of shell and plate structures. Non-uniform rational B-splines (NURBS) [15] have been widely applied in IGA [16–19] and a lot of research work has been conducted on the NURBS-based discretization for Reissner–Mindlin [20,21], Kirchhoff–Love [22,23] shells and composite plates [24–27]. For instance, Kang and Youn [28] applied the method of moving asymptotes to the isogeometric shape optimization of trimmed shell structures based on the NURBS basic function. Bandara and Cirak [29] investigated the isogeometric shape optimization of thin shell structures using subdivision surfaces with a gradient-based shape optimization technique. Ding et al. [30] investigated the accurate analysis and ⁎

thickness optimization of tailor rolled blanks by utilizing three-dimensional NURBS element in IGA. Bhardwaj et al. [31] utilized NURBSbased extended isogeometric analysis (XIGA) to simulate the cracked functionally graded plates using first-order shear deformation theory. But it is difficult to deal with various topological shapes with IGA due to the nature of tensor-product mesh [32]. Gaps and overlaps are inevitable when representing the model with multi-patch surfaces [33]. To overcome the deficiencies of NURBS-based design, a generalization of NURBS surfaces, namely, T-splines, was developed [34]. The novel T-splines support many valuable operations such as local refinement and the merging of multiple surfaces into a single gap-free surface [35,36], which makes it more suitable for the IGA of practical structures. Bazilevs et al. [37] explored T-splines as a basis for IGA at first. Scott et al. [38,39] developed the finite element for T-splines based on Bézier extraction, and combined the collocated isogeometric boundary element method and unstructured analysis-suitable T-spline surfaces for the solution of linear elastostatic problems. Casquero et al. [40] investigated the IGA of fully nonlinear thin shells by utilizing analysissuitable T-spline surfaces of arbitrary degree. Singh et al. [41] developed Bézier extraction based T-spline XIGA for the crack simulations which was more efficient than the NURBS-based one. Gu et al. [42,43] proposed an adaptive XIGA based on local refined B-splines and achieved the faster convergence rate than the non-adaptive one with

Corresponding author. E-mail address: [email protected] (J. Cheng).

https://doi.org/10.1016/j.advengsoft.2019.102692 Received 17 December 2018; Received in revised form 25 June 2019; Accepted 1 July 2019 0965-9978/ © 2019 Elsevier Ltd. All rights reserved.

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Fig. 1. An unstructured T-spline surface: (a) global view; (b) partial view.

uniform global refinement. The approach for handling the unstructured T-spline meshes with extraordinary points is a hot research topic in the field of IGA since the properties of such meshes differ from those without extraordinary points. Several scholars have conducted fruitful research work in this field. For example, Scott et al. [44] proposed the definition of unstructured analysis-suitable T-splines at first and coupled it with isogeometric boudary element methods for linear elastostatic problems. In their oppion, the conditions for unstructured T-spline meshes to be analysis suitable are: no T-junction extensions intersect, no one-bay face extension spans an element in the three-ring neighborhood of an extraordinary point, no extraordinary point lies within the three-ring neighborhood of another extraordinary point. Karčiauskas et al. [45] proposed a new method for complementing bi-cubic splines by biquartic splines near irregularities in T-spline mesh, which could achieve increased smoothness and fast convergence. Toshniwala et al. [46] presented a framework for building smooth splines on meshes with extraordinary points and achieved signifcantly higher convergence rates in some numerical examples such as the Scordelis–Lo roof model etc. But there is still a lack of investigation on the accurate discretization of complex thin shell structures. Although higher continuity within and between isogeometric elements can be achieved by utilizing high-order geometric basis functions, the bandwidth of the sparse stiffness matrices will increase with the increase of the orders of geometric basis functions [47], which will lead to the increase of computational cost in IGA. To improve the computational efficiency of IGA, Li et al. proposed isogeometricmeshfree coupling approach which saved the computational cost based on the realization of local refnement [48,49]. Karatarakis et al. [50] proposed an interactive approach that had the amenability to parallelism especially in massively parallel systems like GPUs, which significantly accelerated the computation of mass and stiffness matrices for single patch. For large structures with multiple patches, the substructuring method, also known as the Schur complement method, has been applied to overcome the difficulty of large dimensionality in analysis by subdividing them into reasonably smaller parts, namely substructures [51]. The method condensed a group of finite elements into a super-element represented as a matrix [52]. Substructure analysis utilizes the matrix reduction technique to reduce the system matrices to a smaller set of degrees of freedom (DOFs), and thus reduces the computing time and allows the solution of very large problems with limited computer resources [53]. Since the NURBS or T-splines basis functions commonly have no property of interpolation, the coupling of substructures in IGA is more complex than in traditional FEA. The purpose of this paper is to propose a new IGA method for large thin shell structures based on the weak coupling of substructures with unstructured T-splines patches. A new Kirchhoff–Love shell element is developed at first based on the unstructured T-splines for the IGA of large thin shell structures considering the tremendous potential of the T-splines-based IGA in handling complex geometries. With the introduction of the substructuring method in IGA, a dimensionality reduction method for the stiffness matrix of an IGA model is proposed based on the condensation of substructures’ DOFs to improve the

computational efficiency. The condensed substructures are represented by a small number of control points (namely, external nodes) of the elements at the coupling interfaces. Since the calculation and condensation of the substructures’ stiffness matrices are completely independent of each other, multi-core parallel computation technology can be easily applied to improve efficiency. And finally, the condensed substructures are integrated into a whole system with the significantly reduced size of stiffness matrix based on a variationally consistent weak coupling method. Illustrative examples demonstrate that the proposed thin-shell element based on unstructured T-splines and the dimensionality reduction method for the stiffness matrix of IGA model are valid and applicable in engineering. The rest of the paper is organized as follows. Section 2 develops the Kirchhoff–Love shell element based on the unstructured analysis-suitable T-splines. Section 3 presents the weak coupling method of substructures. Section 4 proposes the dimensionality reduction method for the stiffness matrix of IGA model while Section 5 verifies the effectiveness of the proposed IGA method by two illustrative examples. And finally, Section 6 draws the conclusions. 2. Kirchhoff–Love shell element based on unstructured T-splines 2.1. Fundamentals of unstructured analysis-suitable T-splines Analysis-suitable T-splines have linearly independent basis functions that are necessary for IGA. Li et al. [54] and Scott et al. [39] proposed the definition of analysis-suitable T-splines for general and unstructured T-mesh, respectively. 2.1.1. Unstructured T-mesh Unstructured T-mesh is the control grid for unstructured T-splines with extraordinary points. Fig. 1(a) is an unstructured T-splines surface with the T-mesh shown in physical space while Fig. 1(b) is the partial view with the T-junctions denoted by red dots and the extraordinary points denoted by blue squares. For the odd degrees, every vertex in the T-mesh is the anchor for a control point PA∈R3, where index A denotes the number of global control point. The unstructured T-splines whose Tmesh satisfies several topological conditions are analysis-suitable [39]. By adding all T-junction extensions to the T-mesh and eliminating all elements whose knot interval sums on any side are zero, the T-mesh is formed. 2.1.2. Bézier extraction for unstructured T-splines The basis functions for standard T-splines without extraordinary points can be defined by the Cox–de Boor recursion formula while those for unstructured T-splines with extraordinary points are relatively difficult to obtain. Bézier extraction that extracts the linear operator for mapping the Bernstein polynomial basis on Bézier elements to the Tsplines basis is applicable to both the standard and unstructured Tsplines [38,39,16] 2

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

B1 ˜

N1e ( ˜) N e ( ˜) =

C1e

Nae ( ˜)

=

B ( ˜)

Ce

(n × (p + 1)2)((p + 1)2 × 1)

=

Cae Cne

Nne ( ˜)

˜1

(1)

+

wn

1 a¯

, (9a)

[Rae,1·(a¯2,2 × a¯ 2) + Rae,2·(a¯1 × a¯ 2,2)]

× a¯ 3) + Rae,2 ·(a¯ 3 × a¯1)] 2 a¯

2Rae,12 ·a¯ 3 +

+

[Rae,1·(a¯1,1 × a¯2) + Rae,2 ·(a¯1 × a¯1,1)]

× a¯ 3) + Rae,2·(a¯ 3 × a¯1)]

a¯3·a¯ 2,2 [Rae,1·(a¯ 2 a¯

D3a =

2·a¯3·a¯1,2 [Rae,1·(a¯ 2 a¯

(9b)

[Rae,1·(a¯1,2 × a¯ 2) + Rae,2 ·(a¯1 × a¯1,2)]

× a¯ 3) + Rae,2 ·(a¯ 3 × a¯1)]

(9c)

= denotes the displacement vector of the element's control points. The 1st and 2nd-order derivatives of the basis function of T-splines in the matrices can be obtained from Eq. (2): [u1e,

ue

(2)

¯ and w are the weights of control points: where w

w ¯ =

Rae,22 ·a¯ 3 +

D2a =

e

u 2e, …, une]T

w eC e (B , W

R,e =

w1 0 0 0 0 w2 0 0 e ,w = 0 0 0 0 0 0 wn

1 a¯

a¯ ·a¯ + 3 a¯1,1 [Rae,1·(a¯2

˜2

w eN e ( ˜) w eC eB ( ˜) = = T e ˜ ¯ N ( ) w ((C e )T w¯ )T B ( ˜)

w1 w2

Rae,11·a¯ 3 +

D1a = B (p + 1)2 ˜

(8)

where

B2 ˜

where ˜ = ( , ) denotes the parametric coordinates defined over interval [−1, 1], n is the number of the control points of element e, Ce 2 denotes the extraction operator of element e, Cae = [ca1, e, ca2, e, ...,ca(p, e+ 1) ] is the operator vector corresponding to the ath control point, and Bi ( ˜) denotes the ith Bernstein polynomial. Then the unstructured T-spline basis functions for element e can be written as:

Re ( ˜)

D1a · e1 D1a · e2 D1a ·e3 D2a ·e1 D2a · e2 D2a ·e3 D3a · e1 D3a · e2 D3a ·e3

Db, a =

BW , )

(10)

W2

R,e = w eC e

B,

B, W , + B, W ,

W

W2

+

2W , W , B W3

BW , W2

(3)

(11) where

2.2. Discretization for unstructured T-splines thin-shell

¯ )T B W = ((C e )T w e T W , = ((C ) w ¯ )T B ,

In IGA, the functions from the geometrical description are utilized as basis functions, thus the displacement field of the unstructured Tsplines element e can be described as: n

u˜ e ( ˜) =

Rae ( ˜) uae =

a=1

n

(p + 1)

a=1

b=1

K=

(4)

(

C

+

t3 12

)

C

,

· a¯ 3

a

,

·a3 = a¯

,

·a¯ 3

a¯

,

· a3

D m, a

=

Rae,1·(a¯1·e3)

Rae,2 ·(a¯ 2·e 1) e (Ra,2 · a¯1 + Rae,1·a¯ 2)·e 1

Rae,2 ·(a¯ 2·e2) e (Ra,2 ·a¯1 + Rae,1·a¯ 2)·e2

Rae,2 ·(a¯ 2· e3) e (Ra,2 ·a¯1 + Rae,1· a¯ 2)· e3

(14b) (14c)

(a¯ 22)2

(14d)

F23 = F32 = a¯ 22·a¯ 12

(14e)

F33 =

v )·a¯ 11·a¯ 22 + (1 + v )·(a¯ 12)2 2

(1

with the contravariant components a¯ calculated by:

a¯

=

/a¯

=

/(a¯ ·a¯ ).

(14f) of the surface metric tensor (15)

where δαβ is the Kronecker delta. 2.3. A quarter sphere patch example

(7)

Deb = [Db,1, with

v )·(a¯ 12)2

F13 = F31 = a¯ 11· a¯ 12 F22 =

where Dem = [Dm,1, ...,Dm, a , ...,Dm, n] denotes the matrix for membrane strains with Rae,1·(a¯1· e2)

(14a)

F12 = F21 = v ·a¯ 11·a¯ 22 + (1

(6)

Rae,1·(a¯1· e1)

Et 3 T (Deb) FDeb dA 12(1 v 2)

(14)

F11 = (a¯ 11)2

= Dem ue =

(Dem)T FDem +

where

where Cαβγδ denotes the elasticity tensor, εαβ denotes the membrane strain, καβ denotes the bending strain, a¯ denotes the basis vector of middle surface in the undeformed configuration, aα denotes the basis vector in the deformed configuration, u is the displacement of middle surface, and the subscript ‘, α’ denotes the derivative with respect to ξα. By introducing the finite element approximation of the displacement field into Eq. (5), the strains can be rewritten as:

Deb ue

v2

F11 F12 F13 F = F21 F22 F23 F31 F32 F33

(5)

u, · a3

1

(13)

d

= 2 (a¯ ·u , + a¯ · u, + u, ·u , ),

e

A

where F is the transformation matrix linking the undeformed configuration to the deformed one:

1

where = a¯

t

Et

K e= e

where Rae ( ˜) denotes the basis function of the ath control point of element e, and uae = [uae (x ), uae (y ), uae (z )] denotes the displacement of the element's control point. Substituting Eq. (4) into the internal work formula of Kirchhoff–Love theory [22], we have:

Wint =

(12)

With Eq. (6) and the internal work from Eq. (5), the stiffness matrix of the unstructured T-splines thin-shell element can be written as:

2

cab, e Bb ( ˜) uae

= ((C e )T w ¯ )T B ,

W,

The quarter sphere shown in Fig. 2 is utilized as an example to demonstrate the validity of the developed Kirchhoff–Love shell element

, D b, a , ...,D b, n] denotes the matrix for bending strains 3

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Fig. 2. The quarter sphere subjected to self-weight g.

based on unstructured T-splines. The quarter sphere, the boundary of which are fixed, is subjected to the self-weight. The geometrical and material parameters of the quarter sphere are as follows: the radius of the sphere R = 10 m, the thickness t = 0.05 m, Young's modulus E = 200 GPa while Poisson's ratio μ = 0.3. The maximum displacements of the model with different numbers of control points are calculated in the convergence analysis. In order to verify the efficiency of the proposed unstructured Tsplines thin-shell element in local mesh refinement, convergence analyses are implemented based on both the unstructured T-splines model and the NURBS model. The former is constructed in the software Rhino and adjusted to meet the analysis-suitble conditons [44] while the latter is constructed utilizing the spline basis function with the C2 continuous 3rd-order polynomial. Fig. 3 illustrates the normalized maximum displacements and relative errors to the reference value obtained by two types of models with different number of control points. As can be observed from Fig. 3, the proposed unstructured T-splines model based on the developed Kirchhoff-Love thin-shell element converges more efficiently than the NURBS model. The improvement mainly results from the local refinement and reduction of superfluous control points in the proposed T-splines model. The unstructured T-splines model based on local refinement with 333 control points and the NURBS model with 841 control points obtained by inserting control points and refining the patch globally are shown in Fig. 4(a) and (c), respectively. The Tjunctions are highlighted as red points while two extraordinary points are highlighted as blue squares in the T-splines IGA model, see Fig. 4(a). The contours of the total displacement obtained by both T-splines and NURBS models are also provided in Fig. 4(b) and (d), respectively.

Fig. 4. IGA of quarter sphere: (a) T-splines model (red points indicate T-junctions while blue squares indicate extraordinary points); (b) contours of the total displacement obtained by T-splines model; (c) NURBS model; (d) contours of the total displacement obtained by NURBS model.

method. It was originally proposed to handle the essential boundary conditions for the Poisson problems without Lagrange multipliers and later extended to handle the coupling problems [55,56]. Yin et al. [57] validated the feasibility of a multiscale XIGA using Nitsche's method and concluded it is the most stable among all domain decomposition methods. Guo and Ruess [58] and Ruess et al. [59] utilized an additional stability term to describe the normal part of the interface traction. The following governing equations for elasticity problem can be extended from the principal of virtual work:

(U + W ) + Wnit = 0 whereWnit =

*

+

*

({ n}· u ) d + N (n ·

*

S

u · u d

u )·( u ·n) d

(16)

where U is the strain energy of the system, W is the potential energy of external loads, Wnit denotes the variationally consistent extension of Nitsche's weak form, σ denotes the stress field, n is the outward unit normal vector along the coupling interface Γ*, τS is a free parameter corresponding to the shear components of the traction and τN is the parameter of the additional stability term. The last two terms of Eq. (16) are introduced to avoid the singular stiffness matrix of the system. The operators 〈 • 〉 and { • } are defined as:

3. Weak coupling of substructures

u = u1

3.1. Weak coupling based on Nitsche's method

u2, { n} =

1n1

+ (1

)

2n2

(17)

where β∈[0,1] determines the contribution of two coupled domains, superscripts 1, 2 refer to the coupled domains as shown in Fig. 5. Then the consistency and stability terms of the Kirchhoff–Love patch coupling can be expressed in terms of the bending moments and force stress resultants [58,60]:

Nitsche's method is a variationally consistent weak coupling

Fig. 3. Convergence curves and relative errors of the normalized displacement.

Fig. 5. Two surface domains coupled along the interface Γ*. 4

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al. c Wnit =

*

({m n n }·

s Wnit =h

+

h3 12

*

*

S

(n )

·

ds +

d

)d

(18)

u · u ds +

S

Fig. 6(a). The patch with the finer meshes is specified as Ω1 to ensure that each boundary element has corresponding Gauss points when two edges have the same length, see Fig. 6(b). The sets of boundary meshes of Ω1 and Ω2 are defined as B1 and B2 respectively. For each Ei∈B1, the Gauss points G = {( ˜j(1) , wj , Ei)} np j = 1 are defined on the coupling interface Γ* based on the quadrature rule for line elements, where np is the number of Gauss points. The procedure of calculating the mapped Gauss points g = {( ˜j(2) , ej )} np j = 1 in the case of 'curve–curve' coupling is as follows:

+ 2b m ) n }· u ) d

m | n + (m n t ), s · u3

*

*

({(n

*

*

N (n ·

N (n ·

u )·( u · n) ds )·(

·n) ds

Step 1. Calculate the physical coordinate Xj and tangent vector t for the given ˜j(1) in Ω1. Step 2. Calculate the vertex coordinates xk1, xk2 of ek∈B2 on the coupling interface. Step 3. Let J1 = Xj−xk1 and J2=Xj−xk2; element ek is the target element if (J1 • t)(J2 • t) < 0, i.e. ej = ek; else return to Step 2. Step 4. Calculate the coordinate ˜j(2) of the mapped Gauss points. The

(19)

c s where Wnit and Wnit are the consistency and stability terms of the Kirchhoff–Love patch coupling, respectively, nαβ and mαβ are the symmetric tensors of stress resultants and stress couples, respectively; nα and tα are the covariant components of the normal and tangent vector of the coupling interface; b is the mixed components of the surface's 2nd fundamental form; h is the thickness of the shell; subscript s preceded by a comma denotes the partial differentiation with respect to the arc length along the boundary. Then the additional stiffness contributions Knit can be obtained by Guo and Ruess [58]: c Knit , rs

Knit =

c s Knit , sr + Knit

problem is to find the root of equation f ( ˜j(2) ) = [x ( ˜j(2) )

(20)

It is worth mentioning that the boundary meshes are determined interactively in the modeling software during the coupling process. Therefore, the boundary meshes can be conveniently determined in the special cases of coupling either partially overlapping edges or variable resolution edges. When the coupling edges are only partially overlapping, only the boundary meshes within the overlapped region are selected and named as B1 and B2 respectively according to the overlapped lengths of two edges, see Fig. 7(a) for illustration. When the resolution of coupling edges is variable, the boundary meshes of each edge are grouped into several subsets at first to the effect that the edge with the finer mesh can be determined in each subset. For instance, the boundary meshes of both coupling edges in Fig. 7(b) are divided into two subsets by the dash line at first, then the finer boundary meshes in each subset can be selected as B1 and B2, respectively.

where and its transpose denote the stiffness contributions s resulting from the consistency terms, Knit denotes the contribution of the stability term. Readers can refer to [58] for more details. c Knit , sr

c Knit , rs

3.2. Mapping of Gauss points on the coupling interfaces The Gauss quadrature rule is adopted for the interfacial integral of boundary elements. The general procedure is to firstly define the Gauss (1) points ˜ at the interface of Ω1 according to the quadrature rule, and i

(2) then to search the corresponding points ˜i mapped on Ω2. Then the interfacial integral can be described as: gp

f (R1, R2) d = *

i=1

(1) (2) f R1 ˜i , R2 ˜i

wi

Xj ]· t = 0 ,

where x ( ˜j(2) ) is the physical coordinate corresponding to ˜j(2) . The equation can be easily solved by Newton's method with the iterative ˜(2) ˜(2) ˜(2) f ( ˜j(2) formula ˜j(2) (k + 1) = j (k ) (k ) )/ f ( j (k ) ) and the initial j (0) = 0 .

(21)

3.2.2. ‘Curve–surface’ coupling In this case, the patch containing ‘curve’ is specified as Ω1 to facilitate the determination of the Gauss points on the coupling interface as shown in Fig. 8, where Xj is the physical coordinate of Gauss point G (2) and x ( ˜j , ˜j(2) ) is the physical coordinate of the mapped point g. The mapping method is similar as the case of ‘curve-curve’ coupling except that the mapping points fall on the surface, i.e. the corre(2) sponding physical coordinate x ( ˜j , ˜j(2) ) is determined by parameters

α

where R denotes the unstructured T-splines basis functions defined over interval [−1, 1]; wi denotes the Gaussian weights. The mapping approach of Gauss points on the ‘curve–curve’ and ‘curve–surface’ coupling interfaces are provided in detail hereinafter considering that an appropriate and efficient mapping method is very important in the implementation of weak coupling. 3.2.1. ‘Curve-curve’ coupling In this case, two patches are attached to each other at their edges as (1) shown in Fig. 6, and the coordinate ˜ is reduced to a constant, i.e.

˜(2) and ˜(2) . The numerical solution procedure of calculating the j j (2)

mapped Gauss points g = {( ˜j , ˜j(2) , ej )} np j = 1 in the case of 'curve–surface' coupling is as follows:

j

˜(1) = (c (1), ˜(1) ) and ˜(2) = (c (2), ˜(2) ) , where c(1) and c(2) take the values j j j j {−1, 1} dependent on the coordinate direction. The patch with the shorter edge is specified as Ω1 to ensure that there will always be corresponding Gauss points on the other coupling edge as shown in

Step 1. Calculate the physical coordinate Xj and tangent vector t for (1) the given ( ˜j , ˜j(1) ) in Ω1. Step 2. Calculate the vertex coordinates xk1, xk2, xk3, xk4 of ek∈B2 on

Fig. 6. Determination of patch Ω1: (a) the patch with the shorter edge; (b) the patch with the finer meshes.

Fig. 7. Determination of boundary meshes B1 and B2 in special cases: (a) partially overlapping edges; (b) variable resolution edges. 5

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Fig. 8. Mapping of Gauss points in 'curve–surface' coupling case.

the coupling interface. Step 3. Let J1 = Xj − xk1, J2 = Xj − xk2, J3 = Xj − xk3 and J4 = Xj − xk4; element ek is the target element if ∀Ji ∈ U= {J1,J2,J3,J4} makes ∀Jj ∈ ∁UJi not at the same side of surface defined by (Ji, t), i.e. ej = ek; else return to Step 2. (2) Step 4. Calculate the coordinate ( ˜j , ˜j(2) ) of the mapped Gauss points. The problem is to find the root of equation (2) (2) f ( ˜j , ˜j(2) , s ) = C (s ) x ( ˜j , ˜j(2) ) = 0 , where C(s) = Xj + s • n is a line along the normal vector n, the length of which is determined by Xj and s. Also, Newton's method is adopted to solve the equation with the iterative formula in Eq. (22) and the initial ˜(2) = 0, ˜(2) = 0, s0 = 0.5. j (0) j (0)

˜(2)

j (k + 1)

˜j(2) (k + 1) sk + 1

Fig. 10. Flowchart of the dimensionality reduction method for the stiffness matrix of an IGA model.

substructuring method, a dimensionality reduction method for the stiffness matrix of an IGA model is proposed based on the condensation of the substructures’ DOFs. In order to achieve satisfactory computational efficiency, it is recommended that all the substructures should have the similar number of elements and control points. The flowchart of the proposed dimensionality reduction method for the stiffness matrix of an IGA model is shown in Fig. 10. 4.1. Construction and condensation of substructures’ stiffness matrices

˜(2) =

j (k )

˜j(2) (k ) sk

where J = dC (s )/ ds,

(2) J 1· f ˜j , ˜j(2) , s

The objective of this step is to firstly divide the complex large structure into several substructures, which can be represented by a single patch or multiple patches, and then condense them into a super element. The stiffness matrix Ks for a single patch can be easily calculated by Eq. (15), while the additional stiffness contributions must be added in the case of multiple patches. The domains containing coupling interfaces are defined as boundary domains while the elements inside the boundary domains are defined as boundary elements e to facilitate the coupling of substructures. For example, a model composed of patches P1, P2, P3 shown in Fig. 11 can be divided into two substructures, one comprises patches P1, P2 while another comprises patch . P3. Then the boundary domains between two substructures are The control points of e are regarded as the external nodes while the remaining nodes are regarded as the interior nodes to be condensed. Thus the equilibrium equation Ku = f of the rth substructure can be written as a partitioned matrix and vectors with the appropriate numbering of nodes:

(22)

(2) x / ˜j (k ) ,

is the Jacobi matrix. x / ˜j(2) (k )

For the case that one edge is to couple multiple surfaces, the boundary meshes of the curve are divided into the same number of subsets as that of the surfaces at first, and then the ‘curve–more surfaces’ coupling problem can be handled as several ‘curve–surface’ coupling problems. Take the coupling of ‘curve–two surfaces’ case shown in Fig. 9 as an example, the boundary meshes of the curve are divided into two subsets at first so that the original ‘curve–two surfaces’ coupling problem are transformed into two ‘curve–surface’ coupling problems with the boundary meshes B1 and B2 determined as illustrated in Fig. 9 for each resulting ‘curve–surface’ coupling problem. As far as the Gauss points are calculated, the numerical results of the interfacial integral can be easily obtained. And thus the additional contribution is ready for the coupling of both patches and substructures composed of multiple patches.

Ks(r )· us(r ) =

(r ) Kbb Kbi(r )

ub(r )

Kib(r ) Kii(r )

ui(r )

=

f b(r ) f i(r )

(23)

That is:

4. Dimensionality reduction method for the stiffness matrix of IGA model

(r ) Kbb ·ub(r )

To reduce the computational cost of the IGA based on

+ Kbi(r ) ·ui(r ) = f b(r )

(23a)

Kib(r ) ·ub(r ) + Kii(r )· ui(r ) = f i(r )

(23b)

where the subscripts b and i indicate the external and interior nodes

Fig. 9. Determination of boundary meshes B1 and B2 in the case of 'curve–more surfaces' coupling.

Fig. 11. The boundary domains and the substructure P3. 6

between the substructure containing P1, P2

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

respectively, u denotes the displacement vector, f denotes the force vector, and Kbb, Kbi, Kib, Kii denote the corresponding blocks. According to Eq. (23b), the displacement vector of the interior nodes ui(r ) can be represented as:

ui(r )

=

1 (Kii(r ) ) (f i(r )

Kib(r ) ub(r ) )

ub(1) ui(1) ub(2) Knit · u (2) i

(24)

ui(n)

1 (Kii(r ) ) Kib(r )

Kbi(r )

(r ) f˜ = f b(r )

Kbi(r ) (Kii(r ) ) f i(r )

ub(1)

(25a)

1

(28)

Then Eq. (29) can be obtained from Eq. (25) as:

(25)

(r ) K˜ (r ) = Kbb

ub Kn* 0 · u i 0 0

ub(n)

By substituting Eq. (24) into Eq. (23a) and rewriting it in the form of Ku = f, we have: (r ) K˜ (r )· ub(r ) = f˜

=

(r ) (Kbb

1

Kbi(r ) (Kii(r ) ) Kib(r ) )· ub(r ) + Kn*(r ) ·

ub(2)

= f b(r )

1

Kbi(r ) (Kii(r ) ) f i(r )

ub(n)

(25b)

(29)

(r ) (r ) where K˜ and f˜ are the condensed stiffness matrix and force vector, (r ) respectively. The stiffness matrix K˜ corresponding to the displacement vector of the external nodes ub(r ) can be obtained by Eq. (25a) while the internal stiffness matrix Kii(r ) has been condensed into the stiffness matrix of a single super element. Then the displacement vector of the external nodes ub(r ) can be calculated by Eq. (25). Matrix inversion is the most computational expensive operation for matrices. The computational cost of matrix inversion is closely related to the order of a matrix, which is the number of nodes as far as a (r ) stiffness matrix is concerned. Supposing that there are Nbb external (r ) (r ) K N nodes and ii interior nodes, the order of stiffness matrix bb and Kii(r ) (r ) are Nbb and Nii(r ) respectively while the order of the uncondensed (r ) + Nii(r ) . The proposed stiffness matrix of the rth substructure Ks(r ) is Nbb (r ) + Nii(r ) order method substitutes the inversion operation of the Nbb (r ) uncondensed stiffness matrix Ks involved in the no substructuring IGA with the inversion operation of the Nii(r ) order stiffness matrix of interior nodes Kii(r ) . As a result, the computational cost for calculating the displacement vector can be slightly saved by the proposed method.

where Kn*(r ) is the corresponding rows in Kn* for f b(r ) , and the equation can be further integrated as:

K˜ (1) K˜ (2) 0

ub(1)

0 + Kn* · K˜ (n)

ub(2) ub(n)

=

(1) f˜ (2) f˜ (n ) f˜

(30)

The equation with the internal nodes of each substructure eliminated has the same form as Eq. (23). That is, the coupling information can be expediently added to the condensed stiffness matrix. As a result, the computational efficiency can be greatly improved due to the reduced scale of IGA and the application of parallel calculations for substructures. The displacement vectors of the substructures’ external control points can be easily calculated by Eq. (30), based on which the displacement vectors of interior points and the displacement field of the entire structure can be obtained by Eq. (24). 5. Numerical examples

4.2. Assembly of global stiffness matrix by coupling the condensed substructures

In the implementation of the proposed IGA method, the large thin shell structures need to be divided into reasonable substructures at first. The following rules are suggested for the division of substructures to ensure the performance of the proposed IGA method. Firstly, the large shell structures should be divided according to their local geometrical features. Secondly, every repetitive structural unit in a large shell structure should be divided as a substructure. Finally, the number of control points should be similar for every substructure if possible so as to fully exploit the advantage of the parallel computing technology in the proposed IGA. Two numerial examples are investigated thoroughly in this section for illustrating the substructing process and demonstraing the validity of the proposed IGA method.

The discrete equation of the entire structure can be written as Eq. (26) considering Nitsche's coupling terms between substructures, which can be further rewritten as the partitioned matrix as Eq. (27).

f s(1)

us(1)

Ks(1) Ks(2)

+ Knit · Ks(n)

(1) Kbb Kbi(1)

Kib(1)

us(2)

=

us(n)

f s(2) f s(n)

0

Kii(1) (2) Kbb Kbi(2)

Kib(2) 0

(26)

+ Knit ·

Kii(2)

ub(1)

f b(1)

ui(1) ub(2) ui(2)

f i(1)

5.1. Black tortoise The black tortoise shown in Fig. 12 is one of the four symbols of Chinese constellations, which is represented by a single 3rd-order unstructured T-splines patch and FEA mesh constructed by shell1811 elements. Fig. 12(b) illustrates a partial view of the unstructured Tsplines model with the extraordinary point denoted by the blue square and the T-junctions denoted by the red points. It is subjected to the selfweight with the displacements at all feet being zero. Fig. 13 illustrates the results of convergence analyses obtained by both the IGA based on unstructured T-splines and the traditional FEA. The maximum

f b(2) =

f i(2)

(n ) Kbb Kbi(n)

ub(n)

f b(n)

Kib(n) Kii(n)

ui(n)

f i(n) (27)

According to Eq. (20), matrix Knit contains only the information of each substructure's external nodes, i.e. the values of rows and columns in Knit related to ui(r ) equal ‘0’. Thus Eq. (20) remains valid in Eq. (27), and the corresponding Knit term in Eq. (27) can be rewritten as:

1 Shell 181 is a 4-node element with 6 degrees of freedom at each node which is suitable for analyzing thin to moderately-thick shell structures.

7

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Fig. 12. The black tortoise subjected to the self-weight: (a) 3rd-order unstructured T-splines model; (b) partial view of the T-splines model with the extraordinary point denoted by the blue square and the T-junctions denoted by red points; (c) FEA mesh model constructed by shell181 elements. Table 1 Number of elements, control points and external points of 7 substructures.

Sub_1 Sub_2 Sub_3 Sub_4 Sub_5 Sub_6 Sub_7 Total

Fig. 13. Convergence analyses of the black tortoise obtained by both the IGA based on unstructured T-splines and the traditional FEA.

Number of elements

Number of control points

Number of external points

8227 1409 1432 1480 1424 2684 3635 –

1441 751 777 822 766 1280 1922 –

125 125 125 125 125 136 761 1522

unstructured T-splines rapidly converges at 7759 control points while the traditional FEA reaches the same convergence value of maximum displacement until the number of nodes increases to 20,456. To verify the validity of the proposed dimensionality reduction method for the stiffness matrix, the IGA model with 7759 control points is selected as the analysis model, which is divided into 7 substructures

displacements of the black tortoise with the different numbers of control points/nodes are calculated. As can be observed from Fig. 13, the maximum displacement obtained by the proposed IGA based on

Fig. 14. Substructuring of the black tortoise model.

8

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Fig. 15. The contours of displacement obtained by: (a) IGA without substructuring; (b) IGA with substructuring; (c) FEA.

represented by different colors as shown in Fig. 14. The boundary doare marked by black and white frames while the corremains sponding external points are denoted by yellow hollow squares. In this example, every substructure is constructed by a single patch, and different substructures attach to each other by ‘curve–curve’ coupling. By applying the proposed dimensionality reduction method, the entire stiffness matrix of the black tortoise can be assembled by the stiffness matrices of condensed substructures with the additional contribution of weak coupling. The number of elements, control points and external points of the 7 substructures are listed in Table 1. The total number of external points equals 1522, i.e. the size of the condensed matrix is reduced to 1522/7759 = 19.6% of its original magnitude. Fig. 15 shows the contours of displacement obtained by the IGA without substructuring, the proposed IGA with substructuring, and FEA, respectively. As can be observed from Fig. 15, the displacement distributions as well as the maximum displacements obtained by three methods are similar to each other. Nevertheless, the time for constructing the finite element model in traditional FEA is significantly longer than IGA because of the time-consuming approximation process between CAD and CAE models, which greatly relates to the proficiency of engineers and the complexity of models. As far as the black tortoise is concerned, it will take more than half an hour to construct its finite element model. Table 2 lists the computing times involved in different methods as well as the resulting maximum displacements of the black tortoise. It is worth noting that the computing time in Table 2 is the total time of model construction and analysis (similarly hereinafter). As can be seen from Table 2, the proposed IGA with substructuring based on matrix reduction is the most efficient among three methods. Specifically, the proposed IGA with substructuring based on matrix reduction improves the computational efficiency by a factor of 176.92/ 96.55 = 1.82 compared to the IGA without substructuring. To further analyze the improvement of computational efficiency, the computing times for assembling the global stiffness matrix and solving the equilibrium equation involved in both the IGA without substructuring and the proposed IGA with substructuring are investigated, the results of which are illustrated in Fig. 16. It should be noted that the computing time for assembling the global stiffness matrix in the proposed IGA with substructuring includes the time involved in the sub-steps of assembling and condensing the stiffness matrices of substructures as well as the coupling of substructures’ stiffness matrices. As can be observed from

Fig. 16. Computing times of different steps involved in: (a) IGA without substructuring; (b) IGA with substructuring.

Fig. 16, the computing time for assembling the global stiffness matrix in the proposed IGA with substructuring is much less than that in the IGA without substructuring. Meanwhile, the computing time for solving the equilibrium equation is reduced by nearly 80% with the introduction of substructuring. Therefore, the conclusion can be made that the proposed IGA with substructuring can simultaneously reduce the computing time for assembling the global stiffness matrix and solving the equilibrium equation.

Table 2 Results comparison of the black tortoise obtained by different methods.

Computing time (s) Maximum displacement (mm)

IGA without substructuring

IGA with substructuring

FEA

176.92 1.163

96.55 1.192

>1800 1.172

Fig. 17. The schematic of Cassegrain design.

9

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Fig. 18. Load condition of the Cassegrain antenna.

Fig. 19. Mesh models of the Cassegrain antenna: (a) T-splines model; (b) FEA model.

Specifically, the calculation and condensation of the substructures’ stiffness matrices can be implemented on multiple processors utilizing the parallel computing toolbox of MATLAB. 5.2. Cassegrain antenna The Cassegrain design shown in Fig. 17 is widely utilized in reflecting antennas, particularly in large antennas such as those in satellite ground stations, radio telescopes and communication satellites. The gravity and wind loads greatly influence the surface precision of the Cassegrain antenna since they will change the amplitude and phase distribution of the antenna's aperture and finally affect the far electrical field of the antenna. Therefore, the shape errors Δ between the ideal main reflector and the deformed one (denoted by red dashed curve) need to be investigated. The antenna composed of the main reflector and a bracket is subjected to the self-weight and a wind load of 20 m/s at the working angle of 45°, see Fig. 18 for illustration. Fig. 19 shows its 3rd-order unstructured T-splines patch and the FEA model with shell181 elements. The whole model can be divided into 16 repetitive substructures constructed by several T-splines patches attaching to each other by ‘curve–curve’ or ‘curve–surface’ coupling. It is worth noting that all the truss rods of the bracket are designed as tubes and thus they are simulated with the developed thin-shell elements in both the IGA and

Fig. 20. A substructure of the Cassegrain antenna.

Moreover, multi-core parallel computing technology can be adopted to further improve the efficiency of the proposed IGA with substructuring since the calculation and condensation of the stiffness matrices of substructures are completely independent of each other. 10

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Fig. 21. The condensation and coupling of the stiffness matrices of the antenna's substructures.

Fig. 23. The shape errors along the radius of the reflector at the angles of 0° and 90°.

As listed in Table 3, the maximum displacement calculated by the proposed IGA method achieves a relatively accurate result with only 24,624 control points and the computing time of 385.6 s. The performance is outstanding compared with the traditional FEA. Moreover, the computational efficiency of the proposed IGA can be further enhanced by calculating and condensing the substructures’ stiffness matrices with multiple processors and the parallel computing toolbox of MATLAB.

Fig. 22. Deformation of antenna under the gravity and wind loads obtained by proposed IGA method.

FEA models in order to realistically represent their connections. There are a total of 26 coupling pairs and 2052 control points for each substructure, see Fig. 20 for illustration. The blue and red points in Fig. 20 indicate the external and internal points of a substurcture respectively. Stiffness matrices of all substructures are sparse, which is plotted in color as shown in Fig. 21. For rotation-free Isogeometric thin-shell analysis, the size of the original stiffness matrix of substructure r is (2052 × 3) × (2052 × 3), which is reduced to (407 × 3) × (407 × 3) after the condensation of its interior points. And finally, the size of the entire stiffness matrix coupled by the condensed substructures is reduced to (407 × 3 × 16) × (407 × 3 × 16). That is, the DOFs of the whole assembly are reduced to less than 20% of its original magnitude. The deformation of antenna under the gravity and wind loads obtained by proposed IGA method is illustrated in Fig. 22. As can be seen from Fig. 22, the maximum displacement equals 0.19 mm, which is a reasonable value. The shape errors along the radius of the reflector at the angles of 0° and 90° are illustrated in Fig. 23, the effects of which on the radiation characteristics can be further investigated.

6. Conclusions A Kirchhoff–Love thin-shell element was developed based on the unstructured analysis-suitable T-splines. The stiffness matrix of the developed thin-shell element can be easily calculated with its displacement field approximated by the unstructured T-splines basis integrated with the internal work formula of Kirchhoff–Love theory. The quarter sphere patch example demonstrated the validity of the developed Kirchhoff–Love shell element. A reduction method for the stiffness matrices of the IGA model for large thin shell structures was proposed based on the substructuring method and the condensation of the substructures’ DOFs. A large thin shell structure was firstly divided into several substructures. The sizes of the stiffness matrices were reduced by condensing the interior control points of all resulting substructures. Then the condensed substructures were coupled into a system with the significantly reduced

11

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al.

Table 3 Comparison of the results obtained by different methods.

Computing time (s) Number of control points/nodes Maximum displacement (mm)

IGA without substructuring

IGA with substructuring

FEA

1243.73 30,272 0.20

385.60 24,624 0.19

>3960 85,568 0.32

size by adopting the variationally consistent weak coupling method. The implementation of ‘curve–curve’ and ‘curve–surface’ coupling of substructures was provided in detail. Two examples demonstrated that the proposed IGA could achieve good accuracy with a significantly improved computational efficiency.

259,680 0.24

421,824 0.21

Reissner–Mindlin shell. Comput Methods Appl Mech Eng 2010;199:276–89. [21] Dornisch W, Klinkel S, Simeon B. Isogeometric Reissner–Mindlin shell analysis with exactly calculated director vectors. Comput Methods Appl Mech Eng 2013;253:491–504. [22] Kiendl J, Bletzinger KU, Linhard J, Wuechner R. Isogeometric shell analysis with Kirchhoff–Love elements. Comput Methods Appl Mech Eng 2009;198:3902–14. [23] Nguyen-Thanh N, Valizadeh N, Nguyen MN. An extended isogeometric thin shell analysis based on Kirchhoff–Love theory. Comput Methods Appl Mech Eng 2015;284:265–91. [24] Shojaee S, Valizadeh N, Izadpanah E, Bui TQ, Vu TV. Free vibration and buckling analysis of laminated composite plates using the NURBS-based isogeometric finite element method. Compos Struct 2012;94:1677–93. [25] Valizadeh N, Natarajan S, Gonzalez-Estrada O A, Rabczuk T, Bui TQ, Bordas SP. NURBS-based finite element analysis of functionally graded plates: static bending, vibration, buckling and flutter. Compos Struct 2013;99:309–26. [26] Yin SH, Yu TT, Bui TQ, Xia SF, Hirose S. A cutout isogeometric analysis for thin laminated composite plates using level sets. Compos Struct 2015;127:152–64. [27] Yu TT, Yin SH, Bui TQ, Hirose S. A simple FSDT-based isogeometric analysis for geometrically nonlinear analysis of functionally graded plates. Finite Elem Anal Des 2015;96:1–10. [28] Kang P, Youn SK. Isogeometric shape optimization of trimmed shell structures. Struct Multidiscip O 2016;53(4):825–45. [29] Bandara K, Cirak F. Isogeometric shape optimisation of shell structures using multiresolution subdivision surfaces. Comput Aided Des 2018;95:62–71. [30] Ding CS, Cui XY, Li GY. Accurate analysis and thickness optimization of tailor rolled blanks based on isogeometric analysis. Struct Multidiscip O 2016;54(4):871–87. [31] Bhardwaj G, Singh IV, Mishra BK, Bui TQ. Numerical simulation of functionally graded cracked plates using NURBS based XIGA under different loads and boundary conditions. Compos Struct 2015;126:347–59. [32] Kim HJ, Seo YD, Youn SK. Isogeometric analysis for trimmed CAD surfaces. Comput Methods Appl Mech Eng 2009;198:2982–95. [33] Sederberg TW, Cardon DL, Finnigan GT, North NS, Zheng J, Lyche T. T-spline simplification and local refinement. ACM Trans Graph 2004;23(3):276–83. [34] Sederberg TW, Zheng J, Bakenov A, Nasri A. T-Splines and T-NURCCS. ACM Trans Graph 2003;22(3):477–84. [35] Wang A, Zhao G, Li YD. An influence-knot set based new local refinement algorithm for T-spline surfaces. Expert Syst Appl 2014;41:3915–21. [36] Dimitri R, Lorenzis LD, Scott MA, Wriggers P, Taylor RL, Zavarise G. Isogeometric large deformation frictionless contact using T-splines. Comput Methods Appl Mech Eng 2014;269:394–414. [37] Bazilevs Y, Calo VM, Cottrell JA, Evans JA, Hughes TJR, Lipton S, et al. Isogeometric analysis using T-splines. Comput Methods Appl Mech Eng 2010;199:229–63. [38] Scott MA, Borden MJ, Verhoosel CV, Sederberg TW, Hughes TJR. Isogeometric finite element data structures based on Bézier extraction of T-splines. Int J Numer Meth Engng 2011;88:126–56. [39] Scott MA, Li X, Sederberg TW, Hughes TJR. Local refinement of analysis-suitable tsplines. Comput Methods Appl Mech Eng 2012;213–216:206–22. [40] Casquero H, Liu L, Zhang Y, Reali A, Kiendl J, Gomez H. Arbitrary-degree T-splines for isogeometric analysis of fully nonlinear Kirchhoff–Love shells. Comput Aided Des 2017;82:140–53. [41] 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. Theor Appl Fract Mec 2017;88:74–96. [42] Gu JM, Yu TT, Lich LV, Nguyen TT, Tanaka S, Bui TQ. Multi-inclusions modeling by adaptive XIGA based on LR B-splines and multiple level sets. Finite Elem Anal Des 2018;148:48–66. [43] Gu JM, Yu TT, Lich LV, Nguyen TT, Bui TQ. Adaptive multi-patch isogeometric analysis based on locally refined B-splines. Comput Methods Appl Mech Eng 2018;339:704–38. [44] Scott MA, Simpson RN, Evans JA, Lipton S, Bordas SPA, Hughes TJR, et al. Isogeometric boundary element analysis using unstructured T-splines. Comput Methods Appl Mech Eng 2013;254:197–221. [45] Karčiauskas K, Nguyen T, Peters J. Generalizing bicubic splines for modeling and IGA with irregular layout. Comput Aided Des 2015;70(C):23–35. [46] Toshniwal D, Speleers H, Hughes TJR. Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: geometric design and isogeometric analysis considerations. Comput Methods Appl Mech Eng 2017;327:411–58. [47] Karatarakis A, Karakitsios P, Papadrakakis M. Computation of the isogeometric analysis stiffness matrix on GPU. In: Papadrakakis M, Kojic M, Tuncer I, editors. Proceedings of the Third South–East European Conference on Computational Mechanics an ECCOMAS and IACM Special Interest Conference. 2013. [48] Li WD, Nguyen-Thanh N, Zhou K. Geometrically nonlinear analysis of thin-shell structures based on an isogeometric-meshfree coupling approach. Comput Methods

Conflict of interest The authors declare that they have no conflict of interest. Acknowledgments This work was supported by the N ational Natural Science Foundation of China (Grant Nos. 51490663, 51475418, 51775491) and the Fundamental Research Funds for the Central Universities (Grant No. 2019FZA4003). References [1] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194(39):4135–95. [2] Cottrell JA, Reali A, Bazilevs Y, Hughes TJR. Isogeometric analysis of structural vibrations. Comput Methods Appl Mech Eng 2006;195:5257–96. [3] Xia BZ, Qin Y, Yu DJ, Jiang C. Dynamic response analysis of structure under timevariant interval process model. J Sound Vib 2016;381:121–38. [4] Natarajan S, Wang JC, Song C, Birk C. Isogeometric analysis enhanced by the scaled boundary finite element method. Comput Methods Appl Mech Eng 2015;283:733–62. [5] Do VNV, Hai OT, Thai CH. Dynamic responses of euler-bernoulli beam subjected to moving vehicles using isogeometric approach. Appl Math Model 2017;51:405–28. [6] Bastle B, Brandner M, Egermaier J, Michálková K, Turnerová E. IGA-based solver for turbulence modelling on multipatch geometries. Adv Eng Softw 2017;113:7–18. [7] Wang ZP, Turteltaub S, Abdalla M. Shape optimization and optimal control for transient heat conduction problems using an isogeometric approach. Comput Struct 2017;185:59–74. [8] Herath MT, Natarajan S, Prusty BG, John NS. Isogeometric analysis and genetic algorithm for shape-adaptive composite marine propellers. Comput Methods Appl Mech Engrg 2015;284:835–60. [9] Kostas KV, Ginnis AI, Politis CG, Kaklis PD. Shape-optimization of 2d hydrofoils using an isogeometric BEM solver. Comput Aided Des 2017;82:79–87. [10] Lieu QX, Lee J. A multi-resolution approach for multi-material topology optimization based on isogeometric analysis. Comput Methods Appl Mech Eng 2017;323:272–302. [11] Wu JL, Gao J, Luo Z, Brown T. Robust topology optimization for structures under interval uncertainty. Adv Eng Softw 2016;99:36–48. [12] Lee SW, Yoon M, Cho S. Isogeometric topological shape optimization using dual evolution with boundary integral equation and level sets. Comput Aided Des 2017;82:88–99. [13] An Z, Yu T, Bui TQ, Wang C, Trinh NA. Implementation of isogeometric boundary element method for 2-d steady heat transfer analysis. Adv Eng Softw 2018;116:36–49. [14] Gong YP, Dong CY, Qu XY. An adaptive isogeometric boundary element method for predicting the effective thermal conductivity of steady state heterogeneity. Adv Eng Softw 2018;119:103–15. [15] Rogers DF, Earnshaw RA. State of the art in computer graphics – Visualization and modeling. New York: Springer; 1991. p. 225–69. [16] Borden MJ, Scott MA, Evans JA, Hughes TJR. Isogeometric finite element data structures based on Bézier extraction of NURBS. Int J Numer Meth Engng 2011;87:15–47. [17] Yuan P, Liu ZY, Tan JR. Shape error analysis of functional surface based on isogeometrical approach. Chin J Mech Eng Eng 2017;30:544–52. [18] Gravenkamp H, Natarajan S, Dornisch W. On the use of NURBS-based discretizations in the scaled boundary finite element method for wave propagation problems. Comput Methods Appl Mech Eng 2017;315:867–80. [19] Lai WJ, Yu TT, Bui TQ, Wang ZG, Curiel-Sosa JL, Das R, Hirose S. 3-D elasto-plastic large deformations: IGA simulation by Bézier extraction of NURBS. Adv Eng Softw 2017;108:68–82. [20] Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. Isogeometric shell analysis: the

12

Advances in Engineering Software 135 (2019) 102692

Z. Liu, et al. Appl Mech Eng 2018;336:111–34. [49] Nguyen-Thanh N, Li WD, Zhou K. Static and free-vibration analyses of cracks in thin-shell structures based on an isogeometric-meshfree coupling approach. Comput Mech 2018;62:1287–309. [50] Karatarakis A, Karakitsios P, Papadrakakis M. GPU accelerated computation of the isogeometric analysis stiffness matrix. Comput Methods Appl Mech Eng 2014;269:334–55. [51] Wang EK, Jagfeld M. Substructuring for large structures using finite element program ANSYS on a workstation. Control and Dynamic Systems 1993;59:137–66. [52] Qu ZQ. Model order reduction techniques: with applications in finite element analysis. New York: Springer; 2004. [53] Biondi B, Muscolino G, Sofi A. A substructure approach for the dynamic analysis of train–track–bridge system. Comput Struct 2005;83(28–30):2271–81. [54] Li X, Zheng J, Sederberg TW, Hughes TJR, Scott MA. On linear independence of Tspline blending functions. Comput Aided Geom D 2012;29:63–76.

[55] Nguyen VP, Kerfriden P, Brino M, Bordas SPA, Bonisoli E. Nitsche's method for two and three dimensional NURBS patch coupling. Comput Mech 2014;53:1163–82. [56] Nguyen-Thanh N, Zhou K, Zhuang X, Areias P, Nguyen-Xuan H, Bazilevs Y, Rabczuk T. Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling. Comput Methods Appl Mech Eng 2017;316:1157–78. [57] Yin SH, Yu TT, Bui TQ, Zheng XJ, Gu ST. Static and dynamic fracture analysis in elastic solids using a multiscale extended isogeometric analysis. Eng Fract Mech 2019;207:109–30. [58] Guo YJ, Ruess M. Nitsche's method for a coupling of isogeometric thin shells and blended shell structures. Comput Methods Appl Mech Eng 2015;284:881–905. [59] Ruess M, Schillinger D, Bazilevs Y, Varduhn V, Rank E. Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method. Int J Numer Meth Engng 2013;95:811–46. [60] Koiter WT. On the mathematical foundation of shell theory. Actes Congrès Intern Math 1971;3:123–30.

13