Coupling of non-conforming trimmed isogeometric Kirchhoff-Love shells via a projected super-penalty approach

Penalty methods have proven to be particularly effective for achieving the required $C^1$-continuity in the context of multi-patch isogeometric Kirchhoff-Love shells. Due to their conceptual simplicity, these algorithms are readily applicable to the displacement and rotational coupling of trimmed, non-conforming surfaces. However, the accuracy of the resulting solution depends heavily on the choice of penalty parameters. Furthermore, the selection of these coefficients is generally problem-dependent and is based on a heuristic approach. Moreover, developing a penalty-like procedure that avoids interface locking while retaining optimal accuracy is still an open question. This work focuses on these challenges. In particular, we devise a penalty-like strategy based on the $L^2$-projection of displacement and rotational coupling terms onto a degree-reduced spline space defined on the corresponding interface. Additionally, the penalty factors are completely defined by the problem setup and are constructed to ensure optimality of the method. To demonstrate this, we asses the performance of the proposed numerical framework on a series of non-trimmed and trimmed multi-patch benchmarks discretized by non-conforming meshes. We systematically observe a significant gain of accuracy per degree-of-freedom and no interface locking phenomena compared to other penalty-like approaches. Lastly, we perform a static shell analysis of a complex engineering structure, namely the blade of a wind turbine.


Introduction
The finite element method (FEM) is a well-established technology for the numerical simulation of a wide variety of engineering applications. It is also well-known that a significant amount of resources is invested into mesh generation and geometry clean-up [20]. To alleviate these burdens, isogeometric analysis (IGA) was introduced in the seminal paper [32]. The main idea of IGA is to improve the interoperability between numerical simulations and Computer Aided Design (CAD) by employing the same mathematical objects used in the geometry description, namely B-splines and non-uniform rational B-splines (NURBS) [46], for the discretization of partial differential equations (PDEs). This shift in paradigm has paved the way for an extensive amount of research, where the reader is referred resulting solution. Recently, these issues have attracted the attention of the research community. In [28], the parameters are scaled by the material properties, the thickness of the shell, the mesh size and by a single, user-defined, problem-independent factor which has been validated on an extensive series of benchmarks. In a similar manner, [45] introduces an additional dependency of the penalty coefficients on various loading and boundary conditions. Concerning the interface locking, a possible remedy based on reduced integration has been proposed in [40]. Nevertheless, the development of a fully parameter-free penalty method which is insensitive to locking phenomena and retains optimal convergence is still in its infancy.
Our contribution seeks to mitigate the aforementioned issues. Inspired by [28] and taking [19] as our starting point, we present an algorithm for achieving displacement and rotational continuity in trimmed shells which is inherently locking-free and where the penalty parameters are automatically defined by the problem setup, namely material properties, geometry of the structure and underlying discretization. Moreover, the aforementioned penalty factors are suitably built to retain the higherorder accuracy of B-splines. The methodology relies on the L 2 -projection of the coupling terms along the associated interface onto a suitable degree-reduced space, where the stable p/p − 2 pairing is employed [13]. We recall that the whole procedure is motivated by the analysis of the underlying perturbed saddle point problems, which gives us insight into the selection of appropriate parameters and into a proper way to eliminate the Lagrange multipliers associated to the constraints. Finally, we highlight that for splines of degree p = 2, 3 the projection is a local operation and therefore it introduces only a small overhead in the total run-time, making it computationally appealing. Then, we verify numerically the robustness of the proposed coupling technique on various nontrimmed and trimmed examples and we compare it with other penalty-like methods. In all cases we observe an optimal convergence behavior, where no boundary locking effects are present. This yields a superior accuracy per degree-of-freedom (dof). Finally, we assess the applicability of our numerical framework to complex engineering structures. To demonstrate this, we perform a static analysis and simplified topology optimization of the DTU 10 MW Reference wind turbine blade [4].
The paper is structured as follows. Section 2 provides to the reader the basic notation related to B-splines whereas Section 3 introduces the fundamentals of trimming. Section 4 explains in details the proposed method, focusing on how to choose the penalty parameters and on how to mitigate locking. In Section 5 the method is validated on a selection of non-trimmed and trimmed numerical benchmarks. Then, our approach is applied to the shell analysis of the DTU 10 MW Reference wind turbine blade, where a simplified topology optimization is performed on the stiffening webs. Finally, some conclusions are drawn in Section 6.

B-splines in a nutshell
In this section, some definitions and fundamentals related to B-splines are reviewed. We refer the reader to [46,20,30], and references therein, for a comprehensive review of B-splines and NURBS and their role in isogeometric analysis. Starting from two integers p, n, a univariate B-spline basis function b i,p of degree p is generated starting from a non-decreasing sequence of real values referred to as knot vector, denoted in the following as Ξ = {ξ 1 , . . . , ξ n+p+1 }. It is worth mentioning that the smoothness of the obtained Bspline basis is C p−k at every knot, where k denotes the multiplicity of the considered knot, while it is C ∞ elsewhere. In the remainder of this work, we consider only splines of maximum continuity, i.e. C p−1 , and degree p ≥ 2. The definition of multivariate B-splines B i,p (η) is achieved in a straight-forward manner using the tensor product of univariate B-splines as: where d denotes the dimension of the parameter space. Additionally, the multi-index i = i 1 , ..., i d denotes the position in the tensor product structure and p = p 1 , ..., p d indicates the vector of polynomial degrees, associated to the corresponding parametric dimension η = η 1 , . . . , η d , respectively. Finally, we denote by Q 0 the Bezier mesh associated to the basis B i,p . Although omitted here, it is straightforward to extend the notation to NURBS, for details see [20]. In the rest of the paper, without loss of generality, the degree vector p will be considered equal in each parametric direction and therefore simplified to a single scalar value p.

Mathematical framework of trimming
In the following section, we summarize the basic mathematical foundation of isogeometric methods defined on trimmed domains, following closely the notation used in [3,17]. For a detailed review of trimming and the current state-of-the-art in IGA we refer to [41] and references therein. Let us define the domain Ω 0 ⊂ R d , where d is the dimension of the physical space of the problem at hand, described by a spline map F : where Ω 0 denotes the parametric domain and we recall that d represents its dimension. In particular, Ω 0 is characterize as a linear combination of multivariate B-spline basis functions and corresponding control points as follows: Then, let us introduce Ω 1 , . . . , Ω N ⊂ R d Lipschitz-regular domains that define the trimming regions to be subtracted from Ω 0 . Consequently, the physical domain reads: where an example is provided in Figure 1 for the case N = 1. We remark that the trimming operation does not change the underlying mathematical description of the original domain. Therefore, elements and associated basis functions are defined with respect to the untrimmed domain Ω 0 . Let us now introduce the B-spline basis of degree p restricted to the corresponding parametric trimmed domain Ω as follows: Similarly, we define the parametric mesh Q as the set of elements such that: where, in the following, we refer to Q as an active cell if Q ∈ Q. Consequently, the definition of physical mesh reads: Finally, we introduce the approximation space formed by multivariate B-splines of degree p restricted to a trimmed domain Ω as follows:

The projected super-penalty method
In this section, we extend the method studied by the authors in [19] for coupling non-conforming Kirchhoff plates to the analysis of trimmed multi-patch Kirchhoff-Love shells. Motivated by the work presented in [13] in the context of isogeometric mortar methods, our strategy leverages the L 2 projection of the coupling terms at the interface, typically defined in terms of the degree p of the solution space related to the corresponding patch, onto a reduced space of B-splines of degree p − 2 defined on the so-called active side of the interface. This procedure mitigates the locking phenomena due to the over-constraint of the solution space in the proximity of the corresponding coupling interface [18]. We remark that our method shares some similarities with the penalty coupling proposed in [40].

A review of differential geometry
Let us review some fundamentals of differential geometry needed to describe the Kirchhoff-Love formulation, following closely the notation in [9]. Recalling the spline geometric mapping of the mid surface F : Ω → Ω, the covariant basis vectors a α are defined as follows: where the comma is used to indicate differentiation with respect to the corresponding curvilinear coordinate. Now, the unit normal vector to the mid surface of the shell a 3 is computed as the normalized cross-product of the in-plane vectors a α : Then, let us introduce the covariant metric coefficients as: Now, the contravariant basis vectors are defined via the following algebraic relationship: where the corresponding covariant and contravariant metric coefficients are linked by the inverse operator: With these coefficients the contravariant basis vectors can be obtained as: We now define the in-plane normal vector n = n α a α to the boundary ∂Ω, where n α denotes its contravariant components. We highlight that n is contained in the tangent plane to the shell. Finally, we define the transformation which maps Cartesian components to curvilinear ones as: where e i represents the standard Euclidean basis.

The weak form of the Kirchhoff-Love shell problem
For the sake of conciseness, in the following we directly work in the discrete setting. The interested reader is referred to [9] for a rigorous derivation. Let us consider as computational domain a manifold Ω ⊂ R 3 with a sufficiently smooth boundary Γ = ∂Ω. Let us split the boundary Γ = ∂Ω into a part associated to Dirichlet-type boundary conditions Γ D = Γ u ∪Γ θ and a part corresponding to Neumanntype boundary conditions Γ N = Γ T ∪ Γ Bnn such that Γ = Γ D ∪ Γ N . It also holds that Γ u ∩ Γ T = ∅ and Γ θ ∩ Γ Bnn = ∅ due to the energetically conjugate nature of applied displacements and transverse shear, and applied rotations and bending moments, respectively. Additionally, let us also introduce the set of corners as χ ⊂ Γ, where this set can be further split into a Neumann part χ N = χ ∩ Γ N and a Dirichlet part χ D = χ ∩ Γ D . Let us also assume a given an applied body forcef ∈ [L 2 (Ω)] d , a prescribed bending momentB nn ∈ L 2 (Γ Bnn ), a prescribed Ersatz forceT ∈ L 2 (Γ T ) as defined in [9] and a given twisting momentS ∈ R for all corners in χ N . With these definitions at hand, the weak formulation of the Kirchhoff-Love shell reads, find u h ∈ V h such that: where the choice of discrete space V h ⊂ S p h d depends in general on the boundary conditions of the problem at hand. The bilinear form a and linear form f can be expanded, respectively, as follows: where we recall that the normal rotation θ n (u) is given as: Then, by leveraging the in-plane projector P = I − a 3 ⊗ a 3 , where I denotes the identity tensor, and the surface gradient ∇ , the membrane and bending strain operators can be defined, respectively, as: where we highlight that the operator sym(·) returns the symmetric part of the input tensor. It is worth noting that the bending operator β requires a global C 1 -continuity of the basis to be welldefined. This is readily achieved within one patch by B-splines of degree p ≥ 2. Next, we can compute the corresponding stress operators by employing a constitutive law. In particular, if we consider a linear elastic model and we analytically integrate through the thickness, we can write: where the fourth-order tensor C for homogeneous materials can be expressed in curvilinear coordinates as: where E and ν represent the Young's modulus and Poisson's ratio, respectively. If we consider composite materials defined as a sequence of orthotropic plies, the bilinear form in (15) must be modified as explained in the following. Let us consider a stacking of plies, numbered by an index n = 1, . . . , P , where P denotes the total number of plies. For each ply we can define the material tensor C n , obtained by transforming the corresponding orthotropic ply tensor from the local ply coordinates to the shell curvilinear reference frame, for further details see [36]. Now, following the classical theory of laminates [49], the homogenized extensional stiffness A, the coupling stiffness B and the bending stiffness D are computed, respectively, as: where t n indicates the thickness of the n-th ply and z n denotes the distance between the centroid of the n-th ply and the mid-plane of the shell, where an example is depicted in Figure 2. Then, the bilinear form associated to a laminate shell reads: where for further details we refer to [36]. Finally, Equation (15) can be summarized in matrix form ζ ply centerline shell mid-plane z n t n Figure 2: Example of a laminate along the thickness direction ζ formed by a non-uniform and nonsymmetric ply sequence. as: where K and f are denoted as the global stiffness matrix and force vector, respectively, and u represents the sought solution coefficients.

The multi-patch setting
Following closely the notation introduced in [13], let us split the computational domain Ω into N non-overlapping subdomains Ω i such that: In CAD terminology, Ω is a B-Rep, i.e. a collection of trimmed surfaces endowed with their topological information. In this work, similarly to [12], we use the so-called face-edge-vertex B-Rep representation. By leveraging the topology, we can then define the interface γ between two adjacent trimmed patches Ω m , Ω n , 1 ≤ m, n ≤ N as a common edge between their faces, see Figure 3 for an example on four patches. Note that two surfaces can share more than one edge. Then, the skeleton Γ is defined as the union of all common interfaces and reads: where L denotes the total number of interfaces and is an ordered index such that 1 ≤ ≤ L.

Remark 1. By a slight abuse of notation, γ can represent both a trimmed or a non-trimmed coupling interface.
Further, let us introduce the cross-points c s as the intersection of at least three shared edges and let us label them with an ordered index c s , s = 1, . . . , S, see again Figure 3 for an illustration.

Remark 2.
It is well-known that CAD softwares provide only an approximation of the true common edge γ which depends on the chosen tolerance. For the sake of simplicity, in our derivation we assume exactness, or equivalently watertightness, of the geometric representation. From a computational standpoint, if the B-Rep is not watertight we perform a closest point projection of the relevant quantities, such as quadrature points and interface knots, onto the coupling edge. For further details we refer to [7]. Now, let us denote by u m the value of the displacement field restricted to Ω m , and similarly u n the value of the primary field on the neighboring subdomain Ω n . Then, for each interface γ the following coupling conditions must be satisfied: which can be rewritten by leveraging the jump operator as:

The projected super-penalty formulation
Following the notation presented in [19], let us introduce for each patch Ω i the following space: Additionally, we denote by V i,h ⊂ Z i,h the finite-dimensional space given by the span of splines associated to subdomain Ω i , where the exact definition of V i,h depends on the set of boundary conditions of the problem at hand. This allows us to introduce the following finite-dimensional space, where we highlight the C 0 -continuity requirement at the cross-points c s . Furthermore, for each interface γ , we introduce the associated knot vector Ξ . The latter is constructed as follows. First, we arbitrarily choose one of the neighboring patches as active. Then, we build Ξ by intersecting the knot lines of the active patch and γ . We highlight that this operation can be performed directly in the parameter space of the active surface, since the B-Rep structure provides a representation of γ in the parameter space of both surfaces. For each surface, we denote the latter representation by γ ( Ω i ), see Figure 5 for an example.

Remark 3.
At this stage, in the spirit of developing a simple and efficient method, we disregard the internal knots of the coupling curve for the construction of Ξ . We highlight that the number of these knots depends on the chosen tolerance in the CAD model, with this number being potentially large. We are aware that this choice could potentially yield a loss of optimality of the method, but for smooth interfaces this effect is negligible. We verify this numerically on two trimmed patches in Figure 4, coupled along a C 1 -continuous quadratic B-spline curve. Indeed, for p = 2, 3 the results are practically indistinguishable, whereas only minor differences are present for the case p = 4. Although outside the scope of this work, finding a simple way to remove this source of sub-optimality constitutes a future research direction.
Then, we build the isogeometric space S p−2 h (γ ) leveraging the p/p − 2 pairing. Assuming Bsplines of maximum smoothness, this space is obtained by removing from Ξ the first and last two knots, where an illustrative example is given in Figure 5 for bivariate B-splines of degree p = 2 and corresponding p−2 = 0 degree-reduced splines defined on the interface knot vector Ξ = [0 1/3 2/3 1].

Remark 4.
The p/p − 2 pairing has been proven to be inf-sup stable in the context of isogeometric mortar methods in [13] and it has been extended to the coupling of non-trimmed Kirchhoff plates in [19]. Although its stability for trimmed geometries has not been rigorously studied, we verify numerically its applicability to the coupling of trimmed Kirchhoff-Love shells.
Consequently, let us define the following space: which is used to characterize the Lagrange multipliers associated to the coupling conditions. We are now ready to define the discretized version of the singularly-perturbed saddle point problem associated to the Kirchhoff-Love shell. Without loss of generality, let us consider homogeneous Neumann-type boundary conditions. Then, the saddle problem reads: where we have introduced the parameters α disp and α rot corresponding to the displacements and normal rotations, respectively. For a rigorous derivation of the singularly-perturbed saddle point formulation in the scope of Kirchhoff plates we refer to [19]. As highlighted in [28,45], these coefficients depend in general on the problem definition, e.g. the material parameters, the thickness of the shell, (a) Initial discretization. The trimmed interface is represented by the dashed red line. Blue dots denote the location of the B-spline curve knots.     Figure 5: Example of the projection setup on a coupling interface for B-splines of degree p = 2. We arbitrarily select the finer mesh (on Ω 1 in this example) to define the projection space S p−2 h (γ ), where the intersections between the parametric coupling curve γ 1 ( Ω 1 ) and the knot lines of Ω 1 are represented by red crosses. Additionally, an intersection mesh at the interface is created only for integration purposes to properly compute the projected penalty terms in Equation (33). The p + 1 integration points are schematically represented by blue dots. Note that we have separated the subdomains for visualization purposes. the applied boundary conditions, the mesh size and discretization degree, where a precise definition of our parameters will be provided in a later section. Let us now eliminate the Lagrange multipliers and rewrite (31) only in terms of the displacement field. In particular, rearranging the second and third equations we obtain: where, with a slight abuse of notation, Π stands for the L 2 -projection, defined on the interface γ , related to the displacements and onto the space S p−2 h (γ ) associated to the normal rotations, respectively. By substituting Equation (32) into the first line of (31) and leveraging the properties of the L 2projection, we obtain: These coupling terms weakly impose the transmission conditions in (27) on the displacements and normal rotations, respectively.

Remark 5.
From a computational standpoint, we rewrite the coupling term associated to the rotations in (33) as defined in [28], where the constraint is recast into two complementary terms. This ensures a non-zero penalty contribution for patches meeting at an arbitrary angle. Then, the L 2 -projection of these terms is performed. For further details, we refer to [28] and references therein.
Now, let us further characterize the aforementioned projection from a computational viewpoint. Let us consider a generic function u ∈ V h defined as the linear combination of basis functions and their corresponding coefficientsû as: Similarly, its projection Π (u) onto the space S p−2 h (γ ) can be written as another linear combination of spline functions and their associated coefficientsũ: The orthogonality of the projection can now be expressed as: which can be rewritten in matrix form by substituting Equations (34) and (35) into Equation (36) as follows: where M denotes the mass matrix associated to the degree-reduced basis and F represents the righthand-side matrix corresponding to the inner product between the basis functions in S p−2 h (γ ) and V h , respectively. In particular, for the projection of the displacement term introduced in Equation (32), F disp is defined as the inner product between the splines in S p−2 h (γ ) d and the jump of the basis functions in V h . Analogously for the rotational term, F rot is assembled as the inner product between the basis functions in S p−2 h (γ ) and the jump of discrete normal rotations in V h . Similarly, we distinguish between the mass matrix M associated to the splines in S p−2 h (γ ) and its vectorial counterpart M corresponding to the functions in S p−2 h (γ ) d . With these definitions at hand, we summarize the computation of the projected terms in Algorithm 1.

Algorithm 1
Computation of the penalty terms in Equation (33).
1: procedure Computation of the penalty terms 2: for each interface γ in Γ do 3: Build the spaces S p−2 h (γ ) and S p−2 h (γ ) d 4: Build the intersection mesh for integration 5:ũ disp ← solve Equation (37)  Lastly, we remark that the solution of Equation (37) is computationally inexpensive for B-splines of degree p = 2, 3 associated to a reduced space of degree p − 2 = 0, 1, respectively, for which the mass matrix is either diagonal or can be lumped.

Selection of penalty parameters
It is well-known that the perturbed problem (31) is variationally consistent only if we select α disp = α rot → ∞ = 1, . . . , L. However, the well-posedness of the underlying problem is insensitive to the choice of the parameters α disp and α rot . Therefore, our method is inherently free from boundary locking, independently of the choice of penalty values, see [19] for further details in the context of isogeometric Kirchhoff plates. This allows us to select α disp and α rot to guarantee the high-order convergence rates achievable by B-splines. Furthermore, in the spirit of developing a parameter-free penalty method, we modify the choice proposed in [28], scaling the displacement and rotation penalty parameters by the physical constants of the underlying problem, the local mesh size, the spline degree and the geometry. For homogeneous isotropic materials they read: where the measure of γ serves as a characteristic length and the exponent β is chosen solely to ensure the optimal convergence of the method. Therefore, it must be a function of the degree p of the underlying discretization. Numerically we have observed that the scaling factor β = p − 1 in (41) is necessary to attain optimal convergence of the method in the H 2 norm, whereas for a scaling of β = p we noticed optimality in the H 2 and H 1 norms. Finally, a factor of β = p + 1 provides optimality in the H 2 , H 1 and L 2 norms. If not stated otherwise, we will use β = p + 1 in all our numerical examples. In case of orthotropic laminates, we adapt the minimum strategy presented in [28], where the minimum local stiffness between adjacent patches Ω m and Ω n is used.
Consequently, the penalty parameters are defined as: Note that all of these parameters are known and depend only on the problem definition, meaning that no user-defined factor is required. Moreover, it is straightforward to check that the penalty terms are dimensionally consistent with respect to their corresponding energy contribution in the weak form (33).

Remark 6.
Clearly, the choice of β influences the condition number of the associated system matrix. This, together with small trimmed elements, can potentially yield ill-conditioned systems of equations and, consequently, loss of accuracy due to numerical round-off errors. In the context of trimmed single-patch shells, a possible remedy based on extended B-splines has been studied in [50]. Furthermore, in the scope of immersed methods, an ad-hoc multigrid preconditioner has been developed in [22]. In this contribution, we employ a direct solver where the stiffness matrix is preconditioned by a simple diagonal scaling. This seems to suffice for the level of accuracy reached in our numerical experiments. We remark that a thorough study of the condition number in the context of trimmed multi-patch Kirchhoff-Love shells is beyond the scope of this paper.

Numerical Examples
In this section we assess the performance of the proposed coupling technique with several numerical examples defined both on trimmed and untrimmed, non-conforming, multi-patch geometries. All the numerical examples presented in the following have been implemented in the open-source and free Octave/Matlab package GeoPDEs [51], where the reparametrization of the trimmed elements for integration purposes is provided by the tool presented in [3]. The analytical shell solutions are taken from the new shell obstacle course studied in [9], where the exact manufactured functions are evaluated in the freely-available Mathematica notebook 1 with 100 digits of precision. Moreover, similarly to [9], for every element we employ 25 × 25 quadrature points to properly capture the highly non-linearity of the quantities of interest. The results of these computations are then imported into GeoPDEs. Also, in all examples taken from [9], we derive from the manufactured solution and apply on the entire boundary ∂Ω non-homogeneous Dirichlet boundary conditions for the displacements and non-homogeneous Neumann boundary conditions for the bending moments. Finally, throughout this section, we compare our choice of penalty factors to a classical approach where the parameters are kept constant: and to the method proposed in [28]: (a) Geometry setup and physical parameters.
x y (b) Initial discretization. Figure 6: Problem setup and initial non-trimmed, non-conforming, multi-patch discretization for the four planar patches example.
where the problem-independent, user-defined parameter δ = 10 3 has been numerically validated on an extensive series of benchmarks.

Four non-trimmed planar patches
The first example is meant to test and verify the implementation of our strategy in a non-trimmed planar setting. The geometrical setup is taken from [19], where this problem is studied in the context of Kirchhoff plates. In particular, the domain Ω = [0, 2]×[0, 2] is subdivided into four non-conforming patches Ω i , i = 1, . . . , 4 coupled along curved interfaces, see Figure 6. To enforce the non-conformity of the latter, the initial interface knots have been shifted by the irrational factor √ 2/100. In our problem definition, we set the Young's modulus E = 10 6 [P a], the thickness of the plate t = 0.005 [m] and the Poisson's ratio ν = 0.3 [−], respectively. Then, to verify the theoretical orders of convergence, we compute the approximation error in the L 2 and H 2 norms with respect to a manufactured smooth solution of the form: The results are summarized in Figure 7, where the convergence of the error measured in the L 2 and H 2 norms, respectively, is plotted against the square root of the number of dofs. We observe that our proposed method attains the expected order of convergence starting from very coarse meshes, whereas interface locking hinders the convergence rates of other penalty methods in the pre-asymptotic regime. As a consequence, we observe a substantial gain of accuracy per degree-of-freedom of the projection strategy, particularly in the L 2 norm. Additionally, we highlight the suboptimal convergence rates achieved by the method proposed in [28], noticeable in the asymptotic regime for p = 4. Further, also for p = 4 and the L 2 norm, we observe the detrimental impact of our choice of penalty parameters on the conditioning of the stiffness matrix and, consequently, on the solution accuracy. This effect can be mitigated by reducing the exponent β in Equation (41), knowing that the method will converge sub-optimally, as depicted in Figure 8. For this reason, we will focus solely on moderate spline degrees p = 2, 3 in the following numerical experiments.

Remark 7.
In order to retain optimal rates of convergence, whenever a cross-point is present in the geometry, we must impose a C 0 -continuity constraint at the cross-point. For further details and a   Comparison of a classic penalty method, the scaled version with respect to the problem parameters proposed in [28] (scaled) and our projection approach (proj).    (a) Geometry setup and physical parameters.  possible implementation, we refer to [19].

Scordelis-Lo roof
In this example we asses the performance of our method on the well-known Scordelis-Lo roof, firstly introduced as part of the shell obstacle course in [8]. The geometrical setup, the chosen parameters and the initial non-conforming multi-patch design, where the roof is split into six subdomains Ω i , i = 1, . . . , 6, are summarized in Figure 9. The structure is supported at both ends of the cylindrical roof by so-called rigid diaphragms, which fix the displacement in the y and z directions, respectively. Moreover, the roof is subjected to a uniform gravity load, directed in the negative z-direction. As studied in [28], we modify the original thickness of the benchmark problem. In particular, we set the Young's modulus, the Poisson's ratio and the thickness of the structure to 4.32 ·  Figure 10b for different penalty methods and also for the single-patch case. We observe that our approach, the method presented in [28] and the single patch case show a similar convergence behavior. However, in case when the penalty parameter is only scaled by the Young's modulus, interface locking phenomena arise, and they are particularly severe for quadratic B-splines. Moreover, in Figure 10a, we compare the time needed to compute and assemble the penalty terms for the aforementioned approaches, where the projection method shows its computational efficiency. This is linked to the fact that, although the projection algorithm requires the solution of an additional system, the corresponding coupling terms involve significantly fewer dofs compared to standard penalty-like methods.

L-beam
This example is meant to demonstrate the applicability of our approach to couple patches at an arbitrary angle, where the corresponding rotational constraint keeps the angle fix during deformation. We consider a beam with an L-section discretized by two non-conforming patches Ω i , i = 1, 2, as    Figure 12a. We remark that on coarse meshes and for the projection method, the rotational constraint is imposed in a less "rigid" way compared to other penalty approaches. This allows to mitigate the effects related to interface locking starting from coarse meshes. Similarly to the previous example, we observe a faster convergence behavior of the vertical displacement under the point load when our approach is employed, see Figure 12b, especially compared to a classical penalty method.

Pure bending of three trimmed planar patches
In this example we consider the computational domain Ω = [0, 2] × [0, 1] split into three trimmed subdomains Ω i , i = 1, 2, 3 as depicted in Figure 13. We remark that in this particular setup, the middle patch is coupled on both sides along trimming interfaces, defined by quadratic spline curves. The applied boundary conditions and loading function are again derived from a smooth solution of the form: Then, we fix the Young's modulus and the Poisson's ratio of the structure to 10 6 [P a] and 0.3 [−], respectively. This example confirms the severity of locking interface phenomena, especially as the shell gets progressively slender. Indeed, we vary the thickness in the range [0.5 , 0.05 , 0.01] [m], where the results are reported in Figure 14. Moreover, we observe that in the trimmed case these detrimental effects are even more pronounced, since more basis functions are involved in the imposition of the constraints compared to the non-trimmed case.     Comparison of a classic penalty method, the scaled version with respect to the problem parameters proposed in [28] (scaled) and our projection approach (proj). (a) Geometry setup and physical parameters.

Remark 9.
In the trimmed case, at any point of a coupling interface γ and for each neighboring patch, we have (p + 1) × (p + 1) shape functions providing a non-zero contribution to the penalty matrices. This is in contrast to the non-trimmed case, where at any point of γ and for each neighboring patch we have at most (p + 1), respectively, 2(p + 1) B-splines involved in the computation of the displacement and rotational coupling terms.

Trimmed astroid
This example is adapted from the shell obstacle course presented in [9]. We consider the computational domain Ω split into three trimmed subdomains Ω i , i = 1, 2, 3 as depicted in Figure 15. In the same figure, the two trimmed interfaces γ , = 1, 2 are defined as quadratic B-spline curves. The domain Ω is characterized by the control points P ij , i, j = 1, . . . , 3 as summarized in Table 1, where the indices i, j are ordered as the parametric coordinates ξ, η represented in Figure 15. The load function and boundary data are computed from the following manufactured solution: where we impose inhomogeneous Dirichlet and Neumann type boundary conditions on the displacements and on the bending moments, respectively, on the entire boundary ∂Ω. Even though this problem is defined on a planar geometry, meaning that bending and membrane responses are decoupled, its investigation is still worthwhile since the solution u is defined as a function of the parametrization. This drastically complicates the derivation of the exact quantities and their stable computation. The convergence results for the error in the L 2 and energy norms for several values of the thickness t = [0.1 , 0.01 , 0.005] [m] are depicted in Figure 16. This example confirms that our       . Comparison of a classic penalty method, the scaled version with respect to the problem parameters proposed in [28] (scaled) and our projection approach (proj). (a) Geometry setup and physical parameters.
x y trimming curves Trimmed multi-patch geometry. Figure 15: Problem setup and initial multi-patch geometry for the trimmed astroid example.
projection method mitigates the detrimental effects linked to interface locking, yielding a significant gain of accuracy per-degrees-of-freedom. This is particularly noticeable as the thickness of the structure becomes smaller, where, for other penalty techniques, locking phenomena hinder the optimal convergence in the pre-asymptotic regime.

Trimmed cylinder
This example is again adapted from the shell obstacle course presented in [9]. We consider the computational domain Ω split into four trimmed subdomains Ω i , i = 1, . . . , 4 as depicted in Figure 17.
where a 3 denotes the covariant vector in the thickness direction. The convergence results for the error measured in the L 2 and energy norms, respectively, are depicted in Figure 18. Similarly to our previous findings, we observe a faster convergence behavior of the projection method in the pre-asymptotic regime, where interface locking is avoided on very coarse meshes. This results in a substantial gain of accuracy per-degree-of-freedom, which is particularly noticeable for quadratic B-splines.

The DTU 10 MW Reference wind turbine blade
In our last example, we perform an isogeometric shell analysis of the DTU 10 MW Reference wind turbine blade [4], whose design was inspired by the NREL 5 MW reference wind turbine [33]. The blade is modeled by 20 non-conforming cubic spline surfaces. As noted in [28], a multi-patch design allows to accurately resolve material discontinuities along the patch interfaces. The outer shell       . Comparison of a classic penalty method, the scaled version with respect to the problem parameters proposed in [28] (scaled) and our projection approach (proj).      Figure 18: Convergence study of the error measured in the L 2 and energy norms for the trimmed cylinder example, B-splines of degree p = 2, 3. Comparison of a classic penalty method, the scaled version with respect to the problem parameters proposed in [28] (scaled) and our projection approach (proj).   Figure 19: DTU 10 MW blade geometry and depiction of the circumferential regions used for the composite materials definition.
of the blade and the internal shear webs are depicted in Figure 19. In the same figure, colored regions are used to define the corresponding composite layup, where each region has a different multidirectional ply stacking sequence and a varying thickness distribution along the spanwise direction. We summarize the most relevant mechanical properties in Table 2. Moreover, in Figure 20, we show the composite layup of the leading panels through the thickness as a function of the spanwise coordinate. For further details on material properties and thickness profiles we refer to [4]. For the analysis, we consider the response of the blade under gravity load, where the blade is modeled as clamped on the rotor side. The L 2 norm of the displacement field and the corresponding deflection of the blade are depicted in Figure 21. The results are obtained by employing a discretization of quadratic B-splines defined on 89528 elements. Note that we directly import the geometry used in [4] for the structural analysis.
Remark 10. The latter is true for every patch except for webs A, B and C, which are obtained by linear extrusion of a generating spline, meaning that one linear element suffices to exactly describe the surface along the corresponding parametric direction. Therefore, h-refinement is performed along the latter direction by introducing 50 equidistributed knots.
We remark that all the results on the blade have been obtained by setting the scaling factor β = p in    the penalty terms to limit the impact of the latter on the condition number of the stiffness matrix.

Simplified topology optimization of webs A and B
This example is meant to show the applicability of the proposed methodology to trimmed geometries obtained by a simplified topology optimization. Note that this numerical test is just a showcase of the flexibility of our computational framework and a realistic topology optimization of the webs is beyond the scope of this work. Furthermore, the geometric operation described in this section are based on an heuristic engineering approach. We trimmed away from the original geometry two holes, one close to the center of the structure and another at the end the web. This design is obtained by adapting the optimized solution presented in [2], where the final geometry is depicted in Figure 22. This design results in a reduction of ≈ 20.6% of the original total mass of web A. Similarly, we perform the same operations on web B. The L 2 norm of the displacement field and the corresponding deflection of the blade with trimmed webs are depicted in Figure 23, where we observe a reduction of ≈ 2.0% in tip displacement related to the loss of total mass of the structure. The results are obtained with quadratic B-splines defined on a total of 84250 active elements.

Conclusions
In this contribution we have extended the methodology presented in [19] to the coupling of trimmed, non-matching surfaces in the context of isogeometric Kirchhoff-Love shells. The strategy is based on the L 2 -projection of suitable penalty terms at the corresponding coupling interface onto a reduced space of degree p − 2 with respect to the chosen approximation space. On one hand, the projection mitigates the detrimental effects related to interface locking starting from very coarse discretization. On the other hand, it gives us insights into the proper scaling of the penalty parameters based on the underlying discretization. Consequently, the proposed coupling method retains the optimal rates of convergence achievable by B-splines, as demonstrated by our findings on an extensive series of benchmark problems. Moreover, building upon the penalty factors studied in [28], our approach is fully parameter-free, since the penalty coefficients are completely determined by the problem setup. Similarly to what was observed in [19], our strategy is particularly suited for spline spaces of moderate degrees p = 2, 3, where the projection turns out to be computationally efficient and the condition number stemming from the super-penalty does not yield a significant deterioration of the solution accuracy. Then, the applicability of our method to tackle complex, industrial-like structure has been studied. In particular, we have performed a static shell analysis of the DTU 10MW Reference wind turbine blade [4], whose design is composed by 20 non-conforming cubic patches. Furthermore, since trimming is naturally incorporated into our methodology, we have carried out a simplified topology optimization of the internal webs. This example demonstrates that the proposed computational framework is readily applicable to tackle industrial optimization loops of engineering relevance.
To conclude, we have numerically validated the wide range of applicability and robustness of the proposed projected super-penalty approach for coupling trimmed Kirchhoff-Love shells, where the method does not suffer from boundary locking and the penalty parameters are automatically derived from the problem setup to attain optimal rates of convergence. Finally, in this contribution only a geometrically linear shell formulation and its static behavior have been considered. Potential future research directions include the extension of the proposed approach to account for geometrical and material non-linearities and the assessment of its performance in dynamics.