Isogeometric analysis for multi-patch structured Kirchhoff-Love shells

We present an isogeometric method for Kirchhoff-Love shell analysis of shell structures with geometries composed of multiple patches and which possibly possess extraordinary vertices, i.e. vertices with a valency different to four. The proposed isogeometric shell discretisation is based on the one hand on the approximation of the mid-surface by a particular class of multi-patch surfaces, called analysis-suitable~$G^1$ [1], and on the other hand on the use of the globally $C^1$-smooth isogeometric multi-patch spline space [2]. We use our developed technique within an isogeometric Kirchhoff-Love shell formulation [3] to study linear and non-linear shell problems on multi-patch structures. Thereby, the numerical results show the great potential of our method for efficient shell analysis of geometrically complex multi-patch structures which cannot be modeled without the use of extraordinary vertices.


Introduction
Isogeometric analysis [4] is a powerful numerical framework for solving partial differential equations by employing the same (rational) spline function space for representing the geometry and the solution space of the considered partial differential equation. In the last years, it has been widely used as an effective tool for the analysis of shells in various engineering disciplines including aerospace, automotive, maritime, civil and biomedical engineering. Thereby, the high continuity of the employed splines allows to achieve a high accuracy and efficiency in the numerical simulation. Examples of studied shells are the Kirchhoff-Love shell [3] and its updated Lagrangian form [5], the Reissner-Mindlin shell [6][7][8][9][10][11], the solid-shell [12][13][14] and hierarchic shells [15,16] for thin, moderately thick and thick shells, respectively. tive, since any (approximate) G 1 -smooth multi-patch surface can be approximated by an AS-G 1 multi-patch surface [2,57].
In this work, we will develop a novel, simple isogeometric method for the analysis of complex Kirchhoff-love shells composed of several patches with possibly extraordinary vertices. The presented technique will follow the second strategy above, and will rely on three main ingredients, namely first on the representation of the mid-surface of the shell by an AS-G 1 multi-patch surface, then on the usage of the globally exactly C 1 -smooth multi-patch spline space [2] as discretisation space of the shell and finally on the application of the Kirchhoff-Love shell formulation [3] for the analysis. Several linear and non-linear benchmark examples will demonstrate the great potential of the presented method for performing analysis of complex Kirchhoff-Love shells.
The outline of the paper is as follows. Section 2 provides the basics of isogeometric Kirchhoff-Love shell analysis with the focus on the used Kirchhoff-Love shell formulation [3]. In Section 3, we introduce the employed C 1 -smooth multi-patch discretisation, which is based on the one hand on the approximation of the mid-surface of the shell by a particular class of multi-patch surfaces, called AS-G 1 multi-patch geometries [1], and on the other hand on the application of the C 1smooth spline space [2]. Section 4 presents the novel isogeometric method for the analysis of multi-patch structured Kirchhoff-Love shells with several linear and non-linear numerical benchmark examples. Finally, we conclude our work in Section 5.

Kirchhoff-Love shell formulation
Based on the works [3,17,19], we will briefly recall the Kirchhoff-Love shell formulation, which will be used throughout the paper. For the sake of brevity, we will restrict ourselves in this section to a single-patch mid-surface r : [0, 1] 2 → R 3 , but the presented formulation can be simply extended to the employed multi-patch setting introduced in Section 3 by just applying it in each case to the single surface patches. Firstly, Section 2.1 will define the shell coordinate system, followed by Section 2.2 where the shell kinematics will be defined accordingly. Finally, in Section 2.3 the variational form of the shell problem will be given. In the following, we will use Greek indices α, β, γ, δ ∈ {1, 2} and Latin indices i, j ∈ {1, 2, 3}.

Shell coordinate system
The Kirchhoff-Love shell element is defined on the surface r(ξ 1 , ξ 2 ) with parametric coordinates ξ α . By the Kirchhoff Hypothesis [58], which implies no shear in the cross section and orthogonality of orthogonal vectors after deformation, any point in the shell y(ξ 1 , ξ 2 , ξ 3 ) can be represented by a point on the mid-surface and a contribution in normal direction: where a 3 is the unit normal vector to the surface and ξ 3 is the through-thickness coordinate. The deformed and undeformed configurations are denoted by y, r andẙ,r, respectively. The covariant basis of the mid-surface and the normal vector are obtained by the partial derivatives of r with respect to its parametric coordinates, i.e.
In addition, the first and second fundamental forms a αβ , b αβ are defined as follows: where a α,β is the Hessian of the surface and a 3,α is the derivative of the normal vector, which can be obtained by Weingarten's formula a 3,α = −b β α a β with b β α = a αγ b γβ as the mixed curvature tensor [59]. From Eq. (2) it can be observed that second-order derivatives are required for evaluation of the second fundamental form.

Shell kinematic and constitutive relations
The Green-Lagrange strains E αβ at any point in the shell are defined as: where ε αβ , κ αβ are the membrane strains and curvature change, respectively, which are obtained by the first and second fundamental forms (2) of the undeformed and deformed configurations [3,17]: Stresses are represented by the stress resultants n and m, corresponding to membrane forces and moments, respectively. Assuming isotropic linear elastic material, they are obtained by where C αβγδ is the plane stress material tensor. Within this paper, only isotropic linear elastic materials are considered, but the formulation can be easily extended to nonlinear materials as shown in [17].

Variational formulation
The variational formulation for the Kirchhoff-Love shell problem is defined by using the principle of virtual work as in [3,17], and by following the notations from [19]. Denoting the internal and external energies by W int and W ext , respectively, the variations of the internal and external work are defined as Here, we denote by u the displacements, by δu the virtual displacements, by δε the virtual membrane strain and by δκ the virtual curvature change. Furthermore,r represents the undeformed mid-surface of the shell, defined on the parametric domain [0, 1] 2 . Note that the metric basis of the deformed configuration is computed on r, which has the same parametric domain as the undeformed configuration. In addition, dΩ = å αβ dξ 1 dξ 2 is the differential area in the undeformed configuration mapped onto the parametric Following a Galerkin approach, we represent the shell displacements u by a finite sum of basis functions φ i and their weights u i . Then, the discrete residual vector R i is defined by taking the first Gateaux derivative with respect to u i [3], i.e.
where Ω ⊂ R 3 represents the surface domain defined by the mid-surface r. Applying a second linearisation of the variational form with respect to u j , the components K i j of the tangential stiffness matrix, can be found as Lastly, for non-linear simulations, Newton iterations are performed for solution u and increment ∆u by solving K∆u = −R.
Throughout these iterations, the deformed geometry r is updated with the displacement field. Due to the appearance of the curvature variations of the curvature tensor κ αβ in Eqs. (5) and (6), second-order derivatives for the basis functions φ i are required for the variational formulation. This means that the basis functions φ i have to be at least globally C 1 -smooth. Considering a multi-patch structured Kirchhoff-Love shell as in this work, the C 1 -smoothness can be trivially satisfied within a single surface patch but has to be imposed across the patch interfaces, which gives rise to the globally C 1 -smooth isogeometric multi-patch spline functions presented in Section 3.

C 1 -smooth multi-patch discretisation space
We will briefly describe the C 1 -smooth multi-patch space construction [2], which will be used as discretisation space for the isogeometric Kirchhof-Love shell analysis of multi-patch geometries with possibly extraordinary vertices in Section 4. Before, we will also present the employed multi-patch structure for the mid-surface r as well as the particular type of multi-patch surface, called AS-G 1 multi-patch surface [1], which is needed to represent the mid-surface r.
The multi-patch surface r defines a surface domain Ω ⊂ R 3 , which can be represented as the disjoint union of the open quadrilateral surface patches Ω (i) , i ∈ I Ω , of open interface and boundary curves Σ (i) , i ∈ I Σ , determined by the boundary curves of the surface patch parameterisations r ( j) , and of inner and boundary vertices x (i) , i ∈ I χ , determined by the corner points of the surface patch parameterisation r ( j) , i.e.
cf. Figure 1. To distinguish the case of an interface or boundary curve and the case of an inner or boundary vertex, we further divide the index sets I Σ and I χ into where in both cases the symbol Γ denotes the boundary case and the symbol • the interface/inner case. In addition, let ν i be the patch valence of the vertex x (i) , i ∈ I Σ .
x 1 x 2 r (1) r (2) r (3) x (1) x (2) x (3) x (4) x (5) x (6) x (7)  Throughout the paper, considering a curve Σ (i) , i ∈ I Σ , or a vertex x (i) , i ∈ I χ , we always locally (re)parameterise (if needed) the corresponding patch parameterisations r (i j ) in the neighborhood of the curve or vertex into the so-called standard form. This means in case of an interface curve that the patch parameterisations r (i 1 ) and r (i 2 ) are (re)parameterised in such a way that the curve Σ (i) is given by Figure 2 (left), and similarly in case of a boundary curve Σ (i) , i ∈ I Γ Σ , by just taking the curve . In case of an inner vertex x (i) , i ∈ I χ , with patch valence ν i , assuming that the patches and interface curves around the vertex x (i) are labeled in counterclockwise order as Σ (i 1 ) , , Ω (i 2ν i ) , the associated patch parameterisations r (i 2 ) , = 1, . . . , ν i , are (re)parameterized in such a way that the interface curve Σ (i 1 ) is given by and the interface curves Σ (i 2 +1 ) , = 1, . . . , ν i − 1, are given by cf. Figure 2 (right). Similarly, the patch parameterisations r (i 2 ) , = 1, . . . , ν i , around a boundary vertex x (i) , i ∈ I Γ Σ , with patch valence ν i are (re)parameterised, where the curve Σ (i 1 ) and the additional curve Σ (i 2ν i +1 ) are the two boundary curves.

AS-G 1 multi-patch surfaces
In this work, we will model the mid-surface r by a specific G 1 -smooth multi-patch surface, called analysis-suitable G 1 (in short AS-G 1 ) multi-patch surface [1]. AS-G 1 multi-patch surfaces are of great importance, since they are needed to ensure the construction of C 1 -smooth isogeometric multi-patch spline spaces with optimal polynomial reproduction properties. More precisely, in case that a multi-patch surface is not AS-G 1 , the approximation power of the resulting C 1 -smooth spline space over the multi-patch surface can be dramatically reduced, cf. [57].
We recall first that a C 0 -smooth multi-patch surface r is G 1 -smooth if and only if for any two neighboring patches Ω (i 1 ) and cf. [60], where r (i 1 ) and r (i 2 ) are the corresponding patch parameterisations of Ω (i 1 ) and Ω (i 2 ) , respectively, given in standard form. Thereby, the functions α (i,i 1 ) , α (i,i 2 ) and β (i) are uniquely 7 determined up to a common function γ (i) : [0, 1] → R via the patch parameterisations r (i 1 ) and r (i 2 ) . A G 1 -smooth multi-patch surface is then called AS-G 1 if for each interface curve Σ (i) , i ∈ I (i) Σ , the functions α (i,i 1 ) and α (i,i 2 ) can be linear polynomials and if there further exists (nonuniquely determined) linear polynomials β (i,i 1 ) : To uniquely determine the functions α (i,i 1 ) , α (i,i 2 ) , β (i,i 1 ) and β (i,i 2 ) , we select those functions α (i,i 1 ) and α (i,i 2 ) , which are relatively prime and minimize and those functions β (i,i 1 ) and β (i,i 2 ) , which minimize cf. [55]. In case of a boundary curve Σ (i) , i ∈ I Σ , we simply set α (i,i 1 ) = 1 and β (i,i 1 ) = 0. The design of AS-G 1 multi-patch spline surfaces was studied so far for the case of planar domains in [1,54,56,57], and for the case of multi-patch surfaces in [2,57]. In this work, we will employ the methods [2,57] to generate AS-G 1 multi-patch parameterisations of the midsurface r, e.g. by approximating a non-AS-G 1 multi-patch parameterisation of the mid-surface r by an AS-G 1 multi-patch geometry.

The specific C 1 -smooth multi-patch spline space
Let the mid-surface r be an AS-G 1 multi-patch surface. The space of C 1 -smooth isogeometric spline functions over r is defined as or equivalently as and cf [2]. For a function φ ∈ V 1 , we denote the equally valued terms r (i,i 1 ) (ξ, 0) by the functions r (i) 0 (ξ) and r (i) 1 (ξ), respectively, where r (i) 0 represents the trace of the function φ along the curve Σ (i) and r (i) 1 describes a specific transversal derivative of the function φ across the curve Σ (i) , cf. [1]. In case of a boundary curve Σ (i) , i ∈ I Σ , we equivalently define the functions r (i) 0 and r (i) 1 just as r (i) Due to the dependence of the dimension of V 1 on the initial geometry of the single patch parameterisations r (i) , i ∈ I Ω , cf. [61], the entire C 1 -smooth space is, even for simple configurations, hard to study and analyze. Therefore, we consider instead the C 1 -smooth subspace A ⊂ V 1 introduced in [2], which is simpler to generate, whose dimension is independent of the initial geometry of the single patch parameterisations r (i) , and which still possesses optimal approximation properties as numerically verified in [2]. The C 1 -smooth subspace A is given as where φ ∈ C 2 T (x (i) ) means that the function φ is C 2 -smooth at the vertex x (i) with respect to the tangent plane at x (i) , cf. [2].
A possible basis of the space A can be given by the set of functions and Φ x (i) collect the basis functions with respect to the individual patches Ω (i) , curves Σ (i) and vertices x (i) , respectively. Thereby, the single C 1 -smooth basis functions φ Ω (i) j 1 , j 2 , j 1 , j 2 ∈ {2, . . . , n − 3}, are fully determined by the properties • φ x (i) j 1 , j 2 • r (i 2m ) ∈ span{N p,s 1 , 2 : 1 = 0, 1, 2 = 0, . . . , n − 3 + 1 and 1 = 0, . . . , n − 3 + 2 , 2 = 0, 1}, m = 1, . . . , ν i , where δ m j is the Kronecker delta, ∂ m 1 , m 1 , m 2 = 0, 1, 2, m 1 + m 2 ≤ 2 represent the six mixed derivative up to order 2 with respect to the tangent plane at the vertex x (i) with appropriate local coordinates x T 1 , x T 2 , and where σ is an appropriate uniform scaling factor, cf. [2]. All C 1 -smooth basis functions are locally supported with a support fully contained within the patch Ω (i) for a function φ Ω (i) j 1 , j 2 , with a support in the vicinity of the curve Σ (i) for a function φ Σ (i) j 1 , j 2 , and with a support in the vicinity of the vertex x (i) for a function φ Σ (i) j 1 , j 2 . For the detailed technical description of the construction of the single C 1 -smooth basis functions, we refer to [2].

Kirchhoff-Love shell analysis of multi-patch structures
Representing the mid-surface of a multi-patch structured Kirchhoff-Love shell by an AS-G 1 multi-patch surface r, see Section 3.2, and employing the globally C 1 -smooth multi-patch discretisation space [2], recalled in Section 3.3, we can directly apply the Kirchhoff-Love shell formulation [3], introduced in Section 2, to the individual surfaces patches r (i) of the multi-patch mid-surface r. This leads to a novel, simple isogeometric method for the analysis of Kirchhoff-Love shells composed of multiple patches with possibly extraordinary vertices. The performance of this method will be tested in this section on the basis of a series of benchmark problems. While in Sections 4.1-4.3 linear Kirchhoff-Love shell problems will be studied, non-linear ones will be considered in Sections 4.4-4.5. For the geometric and material properties of the considered shells we will denote by L, W and a the lengths, by t the thickness, by E the Young's modulus, and by ν the Poisson's ratio.

Paraboloid shell
We represent the mid-surface of a paraboloid shell, described via the parabolic surface   Figure 4 provides the solutions of the paraboloid multi-patch shells y , = 1, 2, with the convergence plots of the solution at the reference point A, see Figure 3, compared to the number of degrees of freedom for the different degrees p by performing uniform h-refinement. Both for the four-patch and for the six-patch case, convergence can be observed for all degrees p towards the reference value u x 3 = 3.93425 · 10 −5 obtained from the corresponding single-patch geometry using a spline space S p,s h with degree p = 8, regularity s = p − 2 and mesh size h = 0.03125.

Hyperboloid shell
We use the same strategy as in Section 4.1 to describe now the mid-surface of a hyperboloid shell, given by the hyperbolic surface with the parameter domain D = [−L/2, L/2] 2 , by three different AS-G 1 multi-patch surfaces r , = 1, 2, 3, see Figure 7 (second and third row). More precisely, we generate an AS-G 1 threepatch surface r 1 , an AS-G 1 four-patch surface r 2 and an AS-G 1 six-patch surface r 3 by constructing the single surface patch parameterisations r (i) via r (i) = r • r (i) , where r (i) are the patch parameterisations of the planar bilinear three-, four-and six-patch parameterisations r 1 , r 2 and r 3 , respectively, shown in Figure 7 (first row). Again, we obtain that r (i) ∈ (S p,s h ) 3 with p = (2, 2), s = (∞, ∞) and h = 1. The shell is fully clamped on the west side at x 1 = −L/2, which means that all displacements and rotations are fixed weakly [30]. We apply a uniformly distributed vertical load with magnitude p x 3 = −8000 · t [N/m 2 ] and perform linear analysis. We study the vertical displacement at point A (x 1 = L/2, x 2 = 0) by performing numerical tests on the three multi-patch shells y , = 1, 2, 3, possessing the mid-surfaces r , for polynomial degrees p = 3, . . . , 8, and maximal 13 possible regularity s = p − 2 for the C 1 spline space A.
The results for the fully clamped hyperboloid shells with the different patch configurations for the mid-surfaces are presented in Figure 6. The convergence plots show for all studied degrees p convergence towards the reference value u x 3 = −0.93137 · 10 −4 obtained from the corresponding single patch geometry using a spline space S p,s h with degree p = 8, regularity s = p − 2 and mesh size h = 0.015625. From the error plots of the three-, four-and six-patch geometry, convergence towards the reference displacement can be seen. Moreover, for the six-patch geometry lower errors can in general be observed compared to the three-and four-patch geometry.

A
A A

Hyperboloid shell with a hole
We consider as in the previous example in Section 4.2 a hyperboloid shell but now with a hole, see Figure 7 (c)-(d). The AS-G 1 four-patch geometry representing the mid-surface of the shell is constructed via trimming from the hyperbolic surface r, given in (8), by following the strategy from [2, Example 2]. For this, we first trim the parameter domain D = [−L/2, L/2] 2 of the hyperbolic surface by cutting out an approximated disk whose boundary is a B-spline curve of degree 2, see Figure 7 (a). The associated untrimmed parameter domain is then described via a planar AS-G 1 four-patch parameterisation r, consisting of individual patch parameterisations r (i) , i = 1, . . . , 4, with r (i) ∈ S p,s h 2 , p = (2, 2), s = (1, 1) and h = 1 4 , see Figure 7 (b). Via r (i) = r • r (i) , i = 1, . . . , 4, we obtain surface parameterisations r (i) ∈ S p,s

Non-linear square plate
To illustrate the performance of the presented method on a non-linear problem, we will model a multi-patch square plate subject to geometric non-linearities. Although the plate will be initially flat, the geometric non-linearity will require evaluation of the solution on the deformed geometry, hence a surface reconstruction of the solution. For this purpose, we consider an initially planar multi-patch square plate of the domain D = [−L/2, L/2] 2 with respect to two different multipatch mid-surface parameterisations r , = 1, 2, where r 1 and r 2 is the three-and four-patch bilinear domain, respectively, given in Figure 9 (first row). Due to the bilinearity of the individual patches, the multi-patch parameterisations r , = 1, 2, are trivially AS-G 1 .
The plate is a square with sides L = 1 [m], thickness t = 0.01 [m] and material parameters E = 2 · 10 11 [N/m 2 ] and ν = 0.3 [−], following a linear Saint-Venant Kirchhoff material law. All sides of the plate are fixed and a uniformly distributed vertical load of p x 3 = −8 · 10 9 [N/m 2 ] is applied. In the numerical tests on the resulting three-and four-patch shells y , = 1, 2, with the mid-surfaces r , we focus on the vertical displacement at the point A (x 1 = 0, y 2 = 0), see Figure 9 (second row), obtained for uniform h-refinement by choosing the polynomial degrees p = 3, 4, and the maximal possible regularity s = p − 2 for the C 1 spline space A. The results are provided in Figure 10. The solutions show fast convergence with respect to the single-patch reference value u x 3 = −0.461767 obtained by using a spline space S p,s h with degree p = 4, regularity s = p − 2 and mesh size h = 0.0625. Contrary to linear shell problems, higher rates of convergence for increasing polynomial order are not visible, and for the current example the difference in speed of convergence for degree 3 and 4 is indistinguishable.

Post-buckling of an L-shaped domain with holes
As last example, we apply the devoloped method to a geometrically non-linear Kirchhoff-Love shell on the post-buckling response of an L-shaped mid-surface with square holes, see Figure 11. The mid-surface is given by the planar bilinearly parameterized 17-patch domain shown in Figure 11, which is AS-G 1 by construction due to the bilinearity of the single patches. Thereby, the two long edges have dimension L = H, the short edges have dimension W and the sides of the squared holes have dimension a = W/2. , which represent the geometric dimensions of the geometry, the thickness, the Young's modulus and the Poisson's ratio, respectively. As can be seen in Figure 11, the displacements of the left short edge of the geometry are fully fixed (u x 1 = u x 2 = u x 3 = 0) while on the bottom edge, the vertical and out-of plane displacements are fixed (u x 2 = u x 3 = 0), keeping horizontal displacements u x 1 17 Furthermore, Crisfield's arc-length method is employed combined with extended iterations to find and apply buckling mode-shapes automatically, hence to find the equilibrium path without applying initial perturbations [63][64][65]. The displacement is evaluated at the three reference points A, B and C visualized in Figure 11 for polynomial degree p = 4 and regularity s = p − 2 of the C 1 -smooth spline space A.
The results for the (post-)buckling analysis of the 17 patch L-shaped domain with square holes are presented in Figure 12 and Figure 13. As can be seen in these figures, the pre-and post-buckling responses are computed with a number of points per branch. Furthermore, the results show that the present method is computing the smooth solutions properly for this nonlinear buckling analysis problem.

Conclusion
We presented an isogeometric method for the analysis of complex Kirchhoff-Love shells, where the mid-surface of the shells is approximated by a particular class of G 1 -smooth multipatch surfaces, called AS-G 1 [1]. This class of multi-patch surfaces allows the use of the globally C 1 -smooth isogeometric multi-patch spline space with optimal polynomial reproduction properties [2] as discretisation space of the considered shell. Our proposed method further relies on the usage of the Kirchhoff-Love shell formulation [3], which can be directly applied due to the Figure 11: Geometry of an L-shaped mid-surface with square holes with the reference points A, B and C for the numerical tests, cf. Figure 13 and 12. The geometry consists of 17 bilinearly parameterized patches. In the performed numerical simulation, the left-side of the domain will be fixed in all directions and the bottom side of the domain will have a sliding boundary condition which will be free in x 1 -direction and will be constrained in x 2 and x 3 directions. The load will be applied to the bottom side with a fixed load f scaled by a magnification factor λ.  globally C 1 -smoothness of the employed discretisation space. We use the developed isogeometric technique to study several linear and non-linear benchmark problems for geometrically complex multi-patch structures. All numerical results show a good convergence behaviour and demonstrate the great potential of the method for the analysis of complex multi-patch structured Kirchhoff-Love shells. One possible future work could be the extension of our approach to an adaptive method to perform local refinement for the analysis of multi-patch Kirchhoff-Love shells. This would allow a reduction of the needed degrees of freedom for solving the problems with an approximation error of similar magnitude.