A framework for efficient isogeometric computations of phase-field brittle fracture in multipatch shell structures

We present a computational framework for applying the phase-field approach to brittle fracture efficiently to complex shell structures. The momentum and phase-field equations are solved in a staggered scheme using isogeometric Kirchhoff–Love shell analysis for the structural part and isogeometric secondand fourth-order phase-field formulations for the brittle fracture part. For the application to complex multipatch structures, we propose penalty formulations for imposing all the required interface constraints, i.e., displacement (C0) and rotational (C1) continuity for the structure as well as C0 and C1 continuity for the phase field, where the latter is required only in the case of the fourth-order phase-field model. All involved penalty terms are scaled with the corresponding problem parameters to ensure a consistent scaling of the penalty contributions to the global system of equations. As a consequence, all coupling terms are controlled by one global penalty parameter, which can be set to 103 independent of the problem parameters. Furthermore, we present a multistep predictor–corrector algorithm for adaptive local refinement with LR NURBS, which can accurately predict and refine the region around the crack even in cases where fracture fully develops in a single load step, such that rather coarse initial meshes can be used, which is essential especially for the application to large structures. Finally, we investigate and compare the numerical efficiency of loosely vs. strongly staggered solution schemes and of the secondvs. fourth-order phase-field models. c ⃝ 2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The modeling and simulation of crack initiation and propagation for structural members is a challenging research topic of high industrial relevance, as it has applications both during the design process and for the inspection and maintenance of in-service structures. The phase-field approach has emerged, in the last two decades, as a promising approach for the modeling and computation of fracture. The approach consists in the approximation of sharp cracks by a continuous field called phase-field, which allows for the description of arbitrarily complex crack patterns. Initially proposed for brittle fracture, the method is based on the variational formulation by Francfort and Marigo [1] of Griffith's theory [2], later regularized by Bourdin et al. [3]. Phase-field fracture models may also be categorized as gradient damage models [4,5], which can recover Griffith's solution through Γ -convergence [6]. A variety of contributions to the original phase-field fracture formulation have been proposed. A fundamental aspect is here the split of the elastic strain energy into "active" and "inactive" (often called "positive" and "negative" or "tension" and "compression") parts to avoid fracture developing under compression, for which different approaches have been proposed. Miehe et al. [7] proposed a split based on the spectral decomposition of the strain tensor, Amor et al. [8] a split based on volumetric and deviatoric parts of the strain tensor. An alternative split based on the symmetric and antisymmetric parts of the strain tensor was proposed by Freddi et al. [9], while Steinke and Kaliske [10] introduced a decomposition of the stress tensor with respect to the crack orientation. Besides the mostly used quadratic degradation function, Borden [11] and Kuhn et al. [12] studied the effects of cubic and higherorder degradation functions. Another fundamental aspect is fracture irreversibility, which was introduced by Miehe et al. [7] through a history field, which corresponds to the maximum of the positive part of the strain energy density, while Gerasimov and De Lorenzis [13] ensure irreversibility through a penalty method. The phase-field approach has been successfully applied to the description of brittle fracture [7,8,14,15], ductile fracture [16][17][18][19], dynamic fracture [20][21][22], fracture in anisotropic media [23,24], and to the description of fatigue failure [25][26][27][28] among many others. Various works have been devoted to the validation of the numerical simulations through experiments [29][30][31][32][33]. Phase-field formulations for fracture have been discretized with standard finite elements as well as isogeometric analysis (IGA). Borden et al. [34] were the first to exploit the continuity properties of isogeometric discretizations in this context and proposed a higher-order phase-field formulation for fracture. Furthermore, phase-field formulations have also been successfully implemented in the commercial finite element software Abaqus [35][36][37][38].
For modeling fracture in thin-walled structures, various approaches have been presented for coupling phasefield models with plate and shell formulations. Solid-shell elements were adopted by Ambati et al. [39], who used through-thickness isogeometric discretization of phase-field and displacement field, and Reinoso et al. [40], who used, for the phase-field variable, a linear interpolation between the top and bottom surfaces. When adopting reduced-order models in which the kinematics is based on midsurface variables, special attention is required for the varying stress distribution through the shell thickness related to bending. Areias et al. [41] used two different phase-fields for describing fracture at the top and bottom of shell surfaces. Amiri et al. [42] and Ulmer et al. [43] used a single phase-field variable defined over the shell midsurface for describing the fracture behavior. However, in the first contribution, the tension-compression split of the elastic strain energy was not considered, while in the second one the split was applied only to the membrane part of the strain tensor. In the work of Kiendl at al. [44] the tension-compression split, based on the spectral decomposition of the strain tensor, was performed at various points through the thickness and numerical integration of the positive and negative strain energy contributions was performed. Paul et al. [45] used a volumetric/deviatoric strain decomposition and employed thickness integration only for the bending part of the strain energy.
In order to correctly resolve the high gradient transition of the phase field around the cracked area, phase-field models require fine meshes at least in the fracture zone. The use of adaptive mesh refinement has been explored with different approaches. Burke et al. [46], Artina et al. [47] and Del Piero et al. [48] showed that local mesh refinement does not influence the propagation of the crack and has little influence both on the energy curves and on the fracture field, justifying the use of adaptive mesh refinement algorithms for saving computational resources when the fracture path is not known in advance. In the first two works, the refinement was defined by an a posteriori error estimate, while the third one used the phase-field value as an indicator for the refinement. The same choice was made by Wick [49] and Heister et al. [50], who focused their attention on a predictor-corrector scheme for the adaptive mesh refinement algorithm, allowing for recomputation of load steps after the refinement. Klinsmann et al. [51] employed an approach for refinement or coarsening of the mesh whose indicator depends on the value of the tension part of the elastic energy, the norm of the phase-field gradient and the element size. Nagaraja et al. [52] used an adaptive multi-level hp-refinement approach.
In the context of IGA, due to the tensor-product structure of NURBS (Non-Uniform Rational B-splines), local refinement is not possible for standard NURBS-based discretizations. Different techniques have been introduced for solving the problem of local refinement, including Hierarchical B-splines [53,54] later developed into Truncated Hierarchical B-splines [55], T-Splines [56,57] and Hierarchical T-splines [58]. Dokken et al. [59] introduced Locally Refined (LR) B-splines, which have been used in the IGA context by Johannessen et al. [60]. LR B-splines have also been shown to produce better matrices, in terms of sparsity and condition number, with respect to other technologies developed for the local refinement of splines [61]. LR B-splines have been extended to LR NURBS by Zimmerman et al. [62]. Adaptive local refinement with T-splines for phase-field fracture analysis was presented by Borden et al. [21] in a simplified version that involved successive runs of the same analysis with meshes refined according to the result of the previous run, thus avoiding the issue of solution transfer between the meshes. Recently, Paul et al. [45] used LR NURBS adaptive mesh refinement for simulating dynamic phase-field fracture in thin shells. Adaptive refinement using Truncated Hierarchical B-splines and Bézier extraction was presented by Kästner and Hennig et al. [63][64][65], in addition to the possibility of mesh coarsening. These works include detailed discussions about the transfer of the solution between the meshes, distinguishing between quantities related to control and integration points. If only refinement is employed, IGA leads to an error-free projection regarding control point variables, while for the quantities related to the Gauss points, different strategies can be used.
The aforementioned contributions focused their attention on rather simple geometries that can be simulated using single patch isogeometric models. However, most of real-world industrial structures present such a complexity that multipatch models with non-matching discretizations at patch interfaces are needed. If higher-order formulations are employed, like Kirchhoff-Love shells [44] for the structural part or the fourth-order phase-field formulation [34] for the fracture part, continuity over patch boundaries needs to be ensured for the primal unknowns and their first derivatives, which poses additional challenges.
In this paper, we present an isogeometric approach for the efficient simulation of brittle fracture in complex multipatch shell structures, using Kirchhoff-Love shells coupled with second-and fourth-order phase-field formulations. In particular, we present an algorithm for adaptive local refinement, which is crucial for efficient phase-field fracture simulations, especially when applied to large structures. We propose a multistep predictor-corrector algorithm for adaptive refinement with LR NURBS, which can accurately predict and refine the region around the crack even in cases where fracture fully develops in a single load step. Furthermore, we present patch coupling formulations for multipatch structures with non-matching meshes. We propose a penalty approach for both structural and phasefield coupling, considering both C 0 and C 1 continuity. The different penalty contributions involved are scaled with the corresponding problem parameters, such that all coupling terms can be controlled by one global penalty parameter, which can be chosen in a problem-independent fashion. This leads to a robust and accurate framework for fracture simulations on complex multipatch structures. Finally, for optimizing the efficiency of the simulations, we investigate and compare the computational cost of loosely vs. strongly coupled staggering strategies and of second-vs. fourth-order phase-field formulations.
The paper is structured as follows. Section 2 reviews the phase-field formulation for brittle fracture and how it is coupled with an isogeometric rotation-free Kirchhoff-Love shell formulation, including the discretization with LR NURBS and the penalty coupling for the extension to multipatch structures. Section 3 contains the details regarding the solution algorithm, with a focus on the adaptive mesh refinement procedure and on the used staggering scheme. The applicability and efficiency of the proposed approach are tested through several numerical experiments in Section 4, ranging from standard fracture benchmarks to complex multipatch problems. Finally, conclusions are drawn in Section 5.

Phase-field model of brittle fracture for shell structures
This section briefly presents the main features of the phase-field approach for modeling brittle fracture and how this formulation can be coupled with plate and shell models, with a focus on Kirchhoff-Love shells which are considered in this contribution.
According to Griffith's theory [2], the equilibrium of a crack is controlled by a potential energy term and a term related to the work required for creating new surfaces. The same two terms are present in the variational formulation of brittle fracture by Franfort and Marigo [1], in which the entire cracking process is described by the minimization of the energy functional The elastic strain energy is computed by integrating its density ψ e , which directly depends on the strain tensor ε, over the domain Ω \ Γ , while the crack energy surface term includes the material fracture toughness G c . The drawbacks of this formulation, specifically the fact that the discrete crack domain Γ ⊂ Ω is unknown and evolves during the analysis, are overcome by the regularized formulation by Bourdin et al. [3]. In this approach, the crack surface energy term is approximated as by introducing a fracture energy density ψ s , defined over the whole domain Ω , which depends on the phase-field variable s. The continuous variation of s, ranging from 1, corresponding to intact material, to 0, corresponding to fully cracked material, approximates the crack topology in the domain. The expression of the fracture energy density originally proposed in [3] ψ s,2 nd (s, ∇s) = G c leads to a strong form of the phase-field evolution equation in which second-order derivatives of s are present, so that the formulation with this choice is termed "second-order phase-field theory". Alternatively, Borden et al. [34] developed a different expression of the fracture energy density which leads to a "higher-order phase-field model" (also known as "fourth-order" formulation, again referring to the order of the derivatives present in the strong form of the phase-field evolution equation). The higher derivatives present in the latter formulation require at least C 1 continuity between elements. The scalar term ℓ in Eqs. (3) and (4) is a length scale parameter which controls the width of the smeared crack in the phase-field approximated model. For ℓ → 0, the approximation converges to the sharp crack solution by Franfort and Marigo (in turn related to Griffith's solution) in the sense of Γ convergence [6]. In order to reproduce the physical asymmetry of the material behavior in tension and compression, the strain tensor ε is additively decomposed into its positive (tensile) and negative (compressive) components: Among the possible approaches proposed in the literature, we adopt a tension-compression split based on the spectral decomposition of the strain tensor as proposed by Miehe et al. [7] in which ε i and n i represent the eigenvalues and eigenvectors of the strain tensor, respectively. ε + and ε − are then obtained from the positive and negative principal strains as having ⟨x⟩ ± = (x ± |x|)/2. According to the split of the strain tensor, the strain energy density and the stress tensor are decomposed into tensile and compressive parts as follows, with λ and µ as the Lamé constants and I as the identity tensor. The positive terms ψ + e and σ + are then degraded by a degradation function g(s): The standard quadratic degradation function, including a positive small factor η ≈ 0 to avoid zero stiffness of the material in a fully cracked state, is adopted: According to the approach exposed above, the energy functional from Eq. (1) can be rewritten in its regularized version as: Fracture irreversibility, meaning that the crack does not heal if external loads are removed, is enforced according to Miehe at al. [7] by replacing ψ + e with the so-called history variable H, defined as the maximum of the positive part of the strain energy density over the pseudo-time of the analysis: A critical feature of the phase-field model in combination with plate and shell formulations is the tensioncompression split introduced in Eq. (5). At each point of the continuum, the strain tensor is defined as the sum of the membrane and bending parts as follows, the latter term depending on the curvature κ and linearly varying along with the thickness coordinate θ 3 . The spectral split requires additional attention because, due to bending, the strains can vary between tension and compression through the shell thickness t, as depicted in Fig. 1. We adopt the approach from Kiendl et al. [44] where, as usually done for Kirchhoff-Love shells, the model is reduced to the behavior of the midsurface variables. So, we define the strain energy surface density Ψ e , which expresses the strain energy per unit area of the midsurface. Including the tension-compression split of the strain energy density, the positive and negative parts of Ψ e are computed as: The dependency of the strain tensor on the thickness coordinate θ 3 , as shown in Eq. (15), leads to a nonlinear distribution of the in-plane stress σ , whose tensile component is degraded according to Eq. (11). For this reason, the integral in Eq. (16) needs to be computed by adopting numerical integration. In each thickness integration point, the total strain is computed as the sum of membrane and bending parts, as in Eq. (15), and the spectral split is then performed. According to this approach, which can be adopted independently of the specific shell formulation, it is possible to describe a nonlinear degradation of stresses and strain energy through the shell thickness (see Fig. 1), assuming only one value of the phase-field variable s and of the degradation function g(s) at the midsurface. The final expression of the energy functional for brittle fracture problems in thin shells, including degradation only of the positive part of the strain energy surface density, becomes where Ψ s is defined similarly to Ψ ± e by through-thickness integration of the fracture energy density ψ s . The stationarity condition of (17) with respect to s leads to the phase-field evolution equation for shells. The weak form of the equation, which will be used for the solution of the system, read, for the second-and fourth-order phase-field formulation respectively, ∫ with δs representing a test function for s.

Isogeometric formulation for Kirchhoff-Love shell
The focus of this paper is on thin shell structures and an isogeometric rotation-free Kirchhoff-Love shell formulation, as described by Kiendl et al. [66], is adopted and extended to local refinement through LR NURBS.

Shell kinematics
In Kirchhoff-Love shell theory, which includes thin plates as a special case, segments initially perpendicular to the shell midsurface remain straight and perpendicular after deformation. Therefore, transverse shear strains can be neglected and the kinematics of the shell can be fully described by the displacement field of the midsurface.
For describing the shell kinematics, a curvilinear coordinate system is considered, with θ 1 and θ 2 as the parametric coordinates used for defining the midsurface, and θ 3 ∈ [−t/2, +t/2] as the through-thickness coordinate. Greek indices α = 1, 2 and β = 1, 2 are adopted for denoting the in-plane components and (·), α = ∂(·)/∂θ α indicates the partial derivatives with respect to θ α . Considering a point r(θ 1 , θ 2 ) on the shell midsurface, a covariant coordinate system can be defined by the tangent base vectors a α = r ,α . The dual contravariant vectors are defined by a α · a β = δ α β , being δ α β the Kronecker delta. The total strain at each point of the shell, see Eq. (15), depends on the membrane strain ε m = ε m αβ a α ⊗ a β and curvature κ = κ αβ a α ⊗ a β , whose covariant components can be expressed as follows, where a 3 represents the unit vector normal to the surface: The expressions ε m and κ are the linearized version of the more general and nonlinear strain measures defined for the considered shell formulation, about which details can be found in [66]. This is appropriate in the context of brittle fracture, where usually failure occurs without the development of large displacements. Isogeometric analysis is a favorable choice for the discretization of such a model because it provides the continuity required by the presence of second-order derivatives in the curvature expression (C 1 continuity).
For implementation and solution of the system of coupled equations, the variational formulation based on the virtual work principle, corresponding to the weak form of the momentum equation, is chosen in which the internal virtual work is defined as follows, with δε and δκ computed according to Eqs. (20) and (21) from δu, which can be interpreted as virtual displacement field. n and m represent the effective stress resultants, for tension and bending respectively, computed by numerical integration of the stresses through the thickness.

Discretization of the geometry using LR NURBS
In traditional NURBS-based isogeometric analysis, surfaces can be parametrized using two parametric coordinates (ξ, η), two sets of knot vectors , a net of n × m control points P i, j , and a tensor product of univariate B-spline basis functions (N i, p and M j,q ) of degree p, q as follows, where w i, j represent the weights associated with each control point. The univariate B-spline basis functions are defined recursively, with respect to the degree p, by the Cox-de Boor formula, starting from piecewise constants ( p = 0): The global knot vector net of lines defines the so-called Bézier mesh of the geometry.
In phase-field fracture analyses, a discretization employing a fine mesh around the crack region is needed to correctly resolve the smeared crack profile, which is often very steep and whose width is controlled by the length scale parameter ℓ (see Fig. 19). The tensor product property of standard NURBS surfaces allows only a global refinement of the geometry, thus resulting in a high number of elements needed and therefore in computationally expensive analyses. For this reason, a local refinement technology needs to be adopted to be able to use small elements only in the crack area. We choose to use LR NURBS [62], an extension of LR B-splines [60]. The idea behind LR splines is the fact that each univariate B-Spline basis function has support in the parametric space only over a limited number of knot spans, i.e. [ξ i , ξ i+ p+1 ], whose knots constitute a "local knot vector" Ξ i . Bivariate LR B-splines are defined as the product of univariate B-splines over local knot spans as follows, and have support over a portion of the domain corresponding to a bivariate local knot vector Ξ k in parametric space. The LR NURBS shape functions are defined as and a surface is parametrized as in which n C P is the number of control points. The tensor product properties are maintained only at the level of the single function, allowing for an unstructured configuration of the control points. The term γ k in Eq. (29) is introduced in order to maintain the partition of unity property after splitting of the shape functions. Without entering into the details of the algorithm for the local refinement of the mesh (which can be found in [60]), a LR B-spline or LR NURBS can be split into two new basis functions by inserting a mesh line over a local knot vector, which generates a new couple of control points and weights (in order to maintain the partition of unity) that will replace the preceding structure. The application and the repetition of this procedure in the two parametric directions of the surface allow the local refinement of the mesh. Different refinement strategies have been proposed in [60]. "Full span" (Fig. 2a) and "minimum span" (Fig. 2b) approaches are based on the choice of an element to be refined. In the first case, all the basis functions having support on the element are split by inserting a couple of mesh lines, one for each parametric direction, that span from the minimum to the maximum knot of the support of all the functions supported by the element to be split. A smaller footprint of the refinement is obtained by adopting the second strategy, in which a couple of as short as possible mesh lines, just long enough to split at least one basis function, are inserted through the marked element center. The "structured mesh" strategy ( Fig. 2c) consists in the split of all the knot spans of the support of a chosen basis function. Among the three methods, the last one is the one used in this work since it provides a more regular mesh when multiple refinement steps are performed and, moreover, keeps the aspect ratio of the elements in the parametric space constant.

Penalty formulations for patch coupling
In order to model the fracture of complex multipatch shell structures, a penalty-based methodology for coupling NURBS patches is hereby presented. The approach couples both the structural and the phase-field behavior across the patch interfaces and uses a single dimensionless penalty parameter α = 10 3 which is scaled with the structural and phase-field problem parameters.

Structural coupling
According to the approach proposed by Herrema et al. [67], we impose displacement and rotational continuity by augmenting the virtual work formulation, see Eq. (23), with new contributions: The term δW pd corresponds to the virtual work from the penalization of the relative displacement between the two patches A and B ( u A − u B ) along the interface curve L: The second added term δW pr , instead, is devoted to preserving rotational continuity by penalizing the relative rotations between the coupled edges of the patches where a n is the unit vector lying in the plane of the patch and orthogonal to L and a 3 is the unit vector perpendicular to the surface. The notation( ·) indicates the undeformed configuration of these geometric variables. For more details regarding the methodology, see Herrema et al. [67]. The formulation works for smooth and non-smooth patch connections as well as for matching and non-matching meshes at the interface. This is an important feature when we employ adaptive local refinement, which might cause different mesh densities at the two sides of a coupling interface even if the meshes were initially matching. Moreover, the penalty formulations can be used for weakly imposing clamping or symmetry boundary conditions at patch edges by considering only the components relative to the first patch in the penalty virtual work contributions. According to [67], the penalty parameters α d and α r are computed from a global penalty parameter α = 10 3 by scaling it with the membrane and bending shell stiffness values, respectively. For the extension to phase-field fracture analysis, we further multiply the penalty parameters α d and α r with the degradation function g(s) to ensure a consistent scaling of structural stiffness and penalty stiffness in the fractured zones. For uniform isotropic material configurations, the proposed scaling is as follows where E is Young's modulus, ν is Poisson's ratio and t is the shell thickness, while h is the average element length along the coupled edge having the finer discretization.

Phase-field coupling between patches
In order to enforce C 0 continuity of the phase-field between the patch interfaces, we add the following term to the left-hand side of the weak form of the phase-field equation (18) and (19). In analogy with the approach adopted for the structural coupling, the term penalizes the difference of s at the two sides of L: If the higher-order phase-field formulation is employed, we propose an approach for imposing the required C 1 continuity of the phase-field across smooth patch interfaces, in addition to the aforementioned C 0 continuity, by penalizing the relative changes in the directional derivative of the phase-field along a n between the two patches: This term needs to be added to the left-hand side of Eq. (19). In (37), we assume that the normal vectors a A n and a B n point in the same direction. The approach is applicable to smooth patch interfaces, while it cannot be extended straightforwardly to patch connections forming a kink. In the latter case, the fact that the vectors a A n and a B n belong to two different planes makes the choice of the directional derivative of s to be penalized ambiguous. In such situations, where the C 1 continuity of the phase-field cannot be imposed, we adopt the second-order phase-field formulation, which requires only C 0 continuity of s, see the example in Section 4.6. The term in (37), considering only the components relative to s A , can be used also for weakly imposing symmetry conditions along one edge.
The penalty parameters are chosen in the same fashion as done for the structural coupling (Section 2.3.1). The global penalty coefficient α = 10 3 is scaled by terms that maintain dimensional consistency with the phase-field equation and ensure that α C 0 P F and α C 1 P F are large enough to guarantee the satisfaction of the imposed continuity constraint without creating ill-conditioning in the phase-field stiffness matrix. The weak forms of the phase-field equations (18) and (19) suggest also the importance of including the term related to the history field H, which is the driving force of the phase-field equation and which becomes numerically predominant when fracture develops. So, for the imposition of phase-field C 0 continuity, the penalty scaling term is defined as follows, while for C 1 continuity we define where H max corresponds to the current maximum value of the history field over all the integration points of the structure.

Adaptive local refinement and staggering schemes
For performing phase-field fracture simulations on complex structures with a reasonable numerical effort, the use of an adaptive local refinement scheme is crucial. In the following, we present a strategy featuring a multistep predictor-corrector algorithm in order to refine the mesh only where needed, i.e. around the crack area, without a priori knowledge of the evolution of the crack. Brittle fracture simulations often involve fast-growing cracks that (fully) develop in a single load step. In this situation, the crack may grow even outside of the mesh region just refined. If this happens, the current load steps or n load steps, including the current one, need to be recomputed with a new further refined mesh, until "convergence of the refinement" is achieved (a typical example of unstable crack growth can be found in Section 4.1).

Adaptive mesh refinement algorithm employing LR NURBS
As indicator for the refinement, we use the value of the phase field s, similarly to Borden et al. [21]. In this work, it was shown that using the quadratic degradation function (12), crack nucleation starts at s = 0.75, and, therefore, s t = 0.75 was also used as threshold value for refinement. In this paper, we use a slightly higher threshold value, namely s t = 0.80. In order to check whether each element needs to be refined, the value of the phase-field s measured at the center of the element is compared with s t . In general, one can state that the higher s t leads to a slightly larger area of refinement, but, on the other hand, it leads to fewer recomputation steps due to crack growth out of the refined zone.
Regarding the refinement typology, we choose to employ the "structured mesh" refinement strategy (Section 2.2.2), for which all the knot spans of the support of certain NURBS basis functions are split. These basis functions are selected among all the NURBS having support on each element marked for refinement as the ones which do not include, in their support, any element having already the minimum mesh dimension (see Section 3.2) for the refinement round. This approach guarantees a regular mesh and a smooth transition between zones with different refinement levels.
The optimal number of steps to be recomputed each time the mesh is refined (n) depends on the type of the problem and the type of solver. When convergence of the staggered iterations is achieved in each step (see Section 3.3), relatively large load steps can be used since their size does not affect the accuracy of the results, but only the frequency in capturing the response of the system. Moreover, the smaller the steps, the higher becomes the computational cost of the analyses, as observed also by Gerasimov et al. [13]. For this reason, we choose to use "large" load increments and to set n = 1. If smaller steps are used, we suggest values of n between 3 and 5 as good alternatives. In case n is chosen larger than 1, the check on the need for mesh refinement is performed every n steps.
The adaptive mesh refinement algorithm is summarized in Fig. 3 in the form used in this work, i.e. with n = 1.

Transfer of field and history variables from coarse to refined mesh
Considering an initial non-refined mesh M 0 and m levels of refinement required, there will be m intermediate refinement rounds producing meshes M 1 , M 2 , . . . , M m−1 , M m , which can be discarded, except for the last one, at the end of the refinement process (see for example Fig. 4). During each refinement round, some basis functions will be refined, i.e. some elements will be split. In order to avoid excessive mesh refinement by splitting elements that have already been refined during the previous steps of the analysis, we set a minimum element dimension for the m th refinement round equal to h 0 /2 m , where h 0 is the characteristic element dimension in the initial non-refined mesh M 0 . The evaluation of the element size is done in the parametric space so that the mesh distortion does not influence the refinement strategy.
For each refinement round, all the state variables defined over the coarse mesh M m need to be transferred to the refined mesh M m+1 . For field quantities defined at the control points (displacement u and phase-field s), the projection occurs according to the same algorithm used for determining the coordinates of the control points in the refined mesh, as outlined in [60]. Regarding the variables stored at the integration points, nominally the history field H, the transfer is based on interpolation of the variable between coarse and refined meshes in a fashion similar to the approached used by Caseiro et al. [68] for transferring strain quantities from integration points to an alternative set of points, see Fig. 5. For each element to be refined, we define locally a set of bivariate Bernstein polynomials B m constituting the basis functions for a Bézier element with the same polynomial order of the adopted NURBS parametrization. So, the value of the history field in a (ξ, η) point of the coarse mesh element can be computed as whereĤ m indicates the (unknown) values of the history field at the control points of the Bézier element. If (ξ m ,η m ) indicates the set of the local coordinates of all the integration points within the Bézier element, the previous expression can be rewritten as: In the latter equation, B m collects the value of the Bernstein polynomials in all the integration points. For each of the sub-elements in which the elements is split, the coordinates of the integration points in the parametric space of the fine element are denoted by (ξ m+1 ,η m+1 ). Analogously to Eq. (41), the projected value of the history field onto the integration points of the refined mesh elements can be found as: Since the refinement occurs by splitting in half each element along both the parametric directions, the coordinates of the integration points are the same in the Bézier parametric space of each element to be split, and so the term B m (ξ m+1 ,η m+1 )B −1 m (ξ m ,η m ) of Eq. (42) is the same for all elements (at each refinement level) and can be precomputed at the beginning of the analysis. For the solution of the non-linear system of coupled equations, the staggered solution idea presented by Bourdin et al. [3] and Miehe et al. [7] consists in solving, for every pseudo-time step of the analysis, first the weak form of the momentum equation (23) for the displacement field u, and then the weak form of the phase-field equation ((18) or (19)) for s. During each solution, the field which is not solved for is kept frozen. We follow the same idea and, after each staggered iteration, we check the convergence of the solution so that the two equations are recomputed, if necessary, until convergence of the staggered iterations is reached. This guarantees that the structure is in equilibrium at the end of each load step. As the momentum equation is nonlinear due to the strain energy split, an iterative solver which finds the equilibrium solution through successive inversion of the structural system is required. In the previously described staggering, which will be called "loosely coupled" scheme, each staggered iteration involves multiple iterations of the momentum equation, followed by a solution of the linear phase-field equation. We adopt instead a "strongly coupled" staggering scheme, in which, during each staggered iteration, only a single iteration of the equilibrium equation is performed. Then, the solution of the phase-field equation is computed, before checking the convergence of the staggered iteration process with an adequately low tolerance for the residual (see Box 1). By reducing the number of iterations of the momentum equation, which are computationally more expensive than the linear solution of the phase-field equation, the strongly coupled staggered approach is able to decrease the computational effort of the simulations, while obviously providing the same result of the loosely coupled staggering scheme, as shown in Section 4.2.

Numerical examples and benchmark tests
The features of the formulation presented in the previous sections will be tested on selected numerical examples. Subsections will be dedicated to the adaptive local refinement algorithm using LR NURBS (Section 4.1), to the loose and strong coupling staggering schemes (Section 4.2), to the standard and higher-order phase-field formulations (Section 4.3) and to the penalty approach for patch coupling (Section 4.4). In these studies, a focus is set on comparing the computational cost of the different approaches. Classical in-plane and plate bending benchmarks will be investigated. The final Sections 4.5 and 4.6 are dedicated to two examples regarding complex shell structures.
For all the simulations, we use quadratic LR NURBS with 3 through-thickness Gauss points for each of the 3 × 3 Gauss points per element. The minimum element size adopted is always equal to h = ℓ/2.

Adaptive local refinement examples
In this section, it will be shown how the adoption of the adaptive local refinement algorithm (Section 3.1) dramatically reduces the computational cost of the analyses and guarantees flexibility in the description of the crack path. Two different examples will be investigated: a plate bending case, in which the phase-field crack nucleates in an unnotched geometry and then fully develops in a single load step, and a 2D plane stress case where the crack evolves and grows gradually during the analysis.
We consider a square plate with geometry outlined in Fig. 6 and material parameters E = 190 × 10 3 N/mm 2 , ν = 0.29, G c = 0.295 N/mm and ℓ = 0.02 mm as in [44]. The plate is simply supported at all the four sides and it is subjected to constant pressure in the out-of-plane direction, resulting in a state of biaxial bending. We compare the results obtained using a uniformly pre-refined NURBS mesh (Fig. 7c) with the ones achieved by employing adaptive local refinement starting from a coarser mesh and prescribing two levels of refinement, in order to obtain the same minimum element size as in the first case. The simulations provide the same result in terms of crack path, as it can be noted by comparing Figs. 7a and 7b, and in terms of load-displacement curve, see Fig. 8. Due to the prescribed convergence of the staggered iterations during the solver, the crack develops in a single load step after the nucleation, as typically happens in the context of brittle fracture. This behavior can be noted from the shape of the load-displacement curves, where a deviation from the linear elastic behavior, corresponding to nucleation, is followed by an abrupt pressure drop, corresponding to instantaneous crack growth. The obtained results are in agreement with previous investigations (see for example [44,69,70]), which show the nucleation of the crack at the center of the plate and its growth towards the corners of the plate. In Fig. 7d the final mesh resulting from the adaptive local refinement algorithm is depicted, with the mesh refined only where needed, i.e. in the crack area. The considerably lower number of elements (Fig. 9a) leads to a significant reduction of computational cost, even if the algorithm requires the step in which the crack develops to be recomputed multiple times (Fig. 9b).  The second example is represented by a square plate containing a notch at one of its sides, where an imposed displacement at the upper edge induces a global tensile or shear state in the specimen. This case, which is a classical benchmark for fracture, has already been investigated in multiple publications (see for example [3,7,13,15,21,[36][37][38][39]64]). We start by considering the shear situation, since in this case the crack grows gradually, while the tensile case will be investigated in the following parts of the paper. Geometry and boundary conditions are shown in   Fig. 10b (the thickness is t = 0.1 mm), while material parameters are E = 1 × 10 9 N/mm 2 , ν = 0.3, G c = 2 N/mm and ℓ = 0.05 mm (as in [44]). Because of the NURBS discretization, which does not allow for discontinuities in the patches, the initial notch needs to modeled through the phase field, as in [21]. For this reason, it is necessary to pre-refine the zone around the initial crack. Three levels of refinements are applied to reach an adequate minimum element size. We want to compare the results of a simulation performed using adaptive local refinement with the ones from a pre-refined mesh. In order to limit the computational time of the analysis with a pre-refined mesh, we use fine elements only in the region where the crack is expected to grow, i.e. the lower-right part of the specimen, making sure to do it in an area large enough not to influence the fracture path (Fig. 11b). The simulations provide the same results in terms of crack path, which is depicted in Fig. 11a. Fig. 11c shows the adaptively refined mesh, having elements of small dimension only around the crack and a smooth transition between fine and coarse mesh. The results of the simulations are the same also in terms of load-displacement curve, see Fig. 12. Regarding the efficiency of the analyses, Fig. 13b and 13c show how the use of the adaptive mesh refinement dramatically reduces the computational time. The relative gain in CPU time for the analysis with adaptive mesh refinement is higher during the first load steps, which do not show any crack development and so do not need to be recalculated, but remains remarkable also in the final steps. The computational cost reduction is obviously related to the lower number of elements required (Fig. 13a) and would be even more evident if all the geometry was to be uniformly pre-refined.

Comparison of loosely and strongly coupled staggering schemes
In the following section, the results and efficiency of the "loose" and "strong" coupling staggering schemes outlined in Section 3.3 are compared. For the sake of brevity, we consider only the square plate example (having  the same input parameters as in 4.1), but the same results have been observed in all the cases where the two approaches have been compared. In order to study only the effect of the staggering strategy, we adopt a uniformly pre-refined mesh as in Fig. 7c and, for both the staggering schemes, the final crack path is the same as the one depicted in the same figure. The load-displacement curves in Fig. 14 confirm that both coupling strategies provide identical results, given that at each load step iterations are performed until convergence. Concerning the efficiency of the simulations, the strong coupling strategy requires more solutions of the phase-field equation but fewer solutions of the equilibrium equation (Fig. 15a). Even if this approach needs more staggered iterations, in each of them only one solving of the structural equation is required. Recalling that the solution of the momentum equation is more computationally expensive than the one of the phase-field equation, it becomes clear why the strong staggering approach reduces the computational time of the analysis, as shown in Fig. 15b. The time saving becomes more evident in the final steps of the simulations, where more staggered iterations are required for convergence. For the sake of efficiency, the strongly coupled staggered scheme is adopted in all the following simulations.

Comparison of second-and fourth-order phase-field formulations
We present a comparison between the results of simulations employing second-and fourth-order phase-field formulations, see Eqs. (3) and (4). The choice of a higher-order formulation is natural in the context of isogeometric analysis, but we want to investigate and compare the efficiency of the two models, especially in the context of adaptive local refinement.
We consider the single-edge notched specimen presented in Section 4.1, this time under tension load, and we repeat the simulation for the two phase-field formulations. The geometry of the specimen and the setup of the analysis are shown in Fig. 10a. The example is solved employing adaptive local refinement of the mesh, which was pre-refined only around the initial crack (Fig. 16a). For both the formulations, the final result of the simulations consists in the abrupt fracture of the specimen due to the straight growth of the crack from the initial notch, as shown in Fig. 16b. The difference in the load-displacement curves, showing a higher critical load for the secondorder phase-field formulation (Fig. 17), is in agreement with what reported by Borden et al. [34] and by Weinberg et al. [71], the first one highlighting that the residual "cohesive tractions" are more accurately described, i.e. they are smaller, in the fourth-order phase-field model.
In the analysis performed adopting the second-order phase-field model, the crack fully develops in two steps, while, in the case of fourth-order formulation, the crack needs only one step for completely propagating.  shows the number of recomputations of the analysis steps corresponding to crack growth (all the steps that are not shown do not need to be recomputed by the algorithm). It is possible to notice that, with the second-order phasefield model, the crack propagation requires a total number of step recomputations, 27, divided into two load steps, much higher than the number of step recomputations required for the higher-order formulation, 16. The reason for this can be found in the higher accuracy in the description of the crack residual stresses of the fourth-order phase-field model, which allows the crack to grow longer for a fixed increase of the imposed displacement. The faster convergence of the higher-order model is confirmed by the lower number of staggered iterations required in each load step, especially in the ones correspondent to the abrupt development of the crack and which therefore need a higher number of iterations (Fig. 18c).
The unidimensional phase-field profiles extracted from the final state of the previous analyses are shown in Fig. 19 and compared with the theoretical ones, which are obtained by analytically solving the unloaded phasefield equations in one dimension after imposing the initial condition s(x = 0) = 0. A closed form exponential solution of these profiles exists for both the phase-field formulations [34]. For the second-order model we have:  while for fourth-order theory: It can be noted that the higher-order formulation shows a slightly narrower profile, meaning that the area that needs to be refined can be reduced. Therefore, the number of elements used for the solution of the fourth-order analysis is decreased, as can be seen in Fig. 18a. The combination of these factors results in a lower total computational time for the higher-order formulation (Fig. 18d), especially in the context of adaptive refinement.

Examples involving penalty patch coupling and its use for modeling pre-cracks
We focus our attention on the penalty patch coupling formulation described in 2.3 and we apply it to the examples presented in the previous section. Adaptive local refinement and both second-and fourth-order phase-field models are employed. For the fourth-order phase-field simulations, also, C 1 continuity of the phase-field across patch boundaries is imposed, according to (37).
We start by considering again the square plate example described in 4.1, but this time the geometry is divided into four patches, as shown in Fig. 20a, which are coupled in their structure and the phase-field by the penalty approach presented. We also consider a setup including only 1 patch, corresponding to one-fourth of the plate, involving two symmetric boundary conditions again imposed by penalty (see Fig. 20b).
The results shown in Fig. 21a for the setup with four patches and considering the second-order phase-field model correspond to those of Fig. 7, which were obtained with a single patch arrangement. The results are identical for the three different setups (one patch corresponding to the full model, four patches, one patch with double symmetry corresponding to one quarter of the model), as it can be observed from the load-displacement curves in Fig. 22, which also confirm that the critical load computed using the higher-order phase-field model is slightly lower than the one predicted by the second-order formulation.  We show that the patch coupling approach and the adaptive local refinement algorithm are also able to replicate the results obtained in single-patch structures not only when the crack nucleates, as in the square plate example, but also in cases where a crack represented by the phase field develops. For this purpose, the single-edge notched example under tensile load described in Section 4.3 is considered. The geometry is subdivided into three patches as shown in Fig. 23a, so that the crack is expected to grow and propagate across the boundary connecting two patches. The outcome of the simulation is displayed in Fig. 23b, which corresponds to the results of the single-patch model (Fig. 16b). Good agreement of the load-displacement curves for the two cases is also shown in Fig. 24.
The penalty-based patch coupling can be used also for modeling pre-existing cracks as discontinuities in the structures. This idea will be shown by considering the example of the single-edge notched specimen, both in tension and shear, introduced in Section 4.1, whose behavior is determined by the presence of an initial crack. Since we want to compare our results with the ones obtained with classical finite element theory and due to the lower continuity of the Lagrange polynomials, only the second-order phase-field formulation can be considered. The division of the geometry in patches in order to allow the discrete modeling of the pre-crack (along which no patch coupling occurs) is shown in Fig. 25. The mesh is pre-refined at the discrete crack tip, in order to better describe the phase-field crack which is expected to grow from this location. Fig. 26 shows the final stages of the analyses with the pre-crack modeled by using patch discontinuity. The load-displacement curves in Fig. 27 highlight how the different treatment of the pre-crack affects the prediction of the critical fracture load. Regarding the tensile case, it appears that a discrete pre-crack is able to better model the brittle behavior of the structure, with the crack developing and growing in only one step. A difference in the critical load is present also for the shear case, possibly due to residual stresses in the phase-field pre-crack. For comparison, we computed both cases also with C 0 finite elements, in whose framework pre-cracks are usually modeled as discontinuities in the mesh that can be easily handled without the use of any patch coupling technology. Linear shape functions were employed with minimum element dimension h = ℓ/4. The results show that the IGA discrete crack approach featuring penalty coupling of patches can reproduce the load-displacement curves of the classical finite element framework. Another advantage of the discrete crack modeling is the fact that it does not need pre-refinement of the area in which the initial crack is defined (see for example Fig. 16a for comparison with phase-field induced pre-crack), requiring fewer elements and so allowing for faster analyses. However, the modeling of pre-cracks with patch discontinuities needs to be treated carefully since it might involve contact between the uncoupled patches edge, i.e. the crack surfaces. The current approach does not prevent penetration and a contact formulation should be added to the model for avoiding the problem and correctly simulate critical loads. For the examples presented in this work no penetration of crack surfaces happens, but this should be checked when the methodology is applied. In case the pre-crack is modeled using phase-field, the strain energy split and the degradation only of the "tensile" component avoid the penetration problem.

Multipatch example with smooth patch interfaces: pressurized cylinder with hemispherical caps
We apply the presented formulation and solution algorithm to a shell example featuring a pressurized cylinder with hemispherical end caps and an initial crack. The example resembles the one considered by Borden et al. [21], who simulated the dynamic brittle fracture of the structure using solid elements. In our case, the shell formulation previously introduced is employed, since it can be favorably used for the discretization of the geometry, as shown in Fig. 28. Two alternative geometries are considered, one with uniform shell thickness equal to 10 mm and one with thicker hemispherical caps having t = 18 mm. We take advantage of the symmetry of the problem for modeling only half of the structure (the symmetry condition is imposed by penalty). We use the multipatch approach for modeling the structure and the initial notch with a discrete pre-crack represented by a discontinuity between the patches (see Fig. 29). The material parameters adopted for the analysis are E = 190 × 10 3 N/mm 2 , ν = 0.3, G c = 22.13 N/mm and ℓ = 2.5 mm. Since all the patch interfaces are smooth, the fourth-order phase-field formulation, with C 0 and C 1 continuity of s imposed by penalty between patch boundaries, is adopted. Figs. 30 and 31 show the final crack state for the two situations. In the uniform thickness case, the crack grows straight towards the center of the hemispherical caps, while in the setup with thicker caps the crack branches and grows in circumferential direction along the juncture of cylinder and end-caps. Fig. 32 displays the load-displacement curves for the two cases. 4.6. Multipatch example with non-smooth patches interfaces: two-field beam with I-profile As a final example, we consider a beam having I-shaped cross section, with three supports and subjected to two point loads, each of them acting at the middle of the two beam sections. Geometry and setup of the problem are shown in Fig. 33, while material parameters are E = 190 × 10 3 N/mm 2 , ν = 0.30, G c = 0.295 N/mm and ℓ = 7.5 mm. Each section of the beam is modeled using three patches, one for the web and two for the flanges. In this   situation, the imposition of C 1 continuity for the phase-field at the interface between web and flange patches is not possible using the proposed approach and the term in (37), because the direction of the a n vector cannot be chosen univocally for the patches corresponding to the flanges. Accordingly, we adopt for this example the second-order phase-field formulation, which requires only C 0 continuity of s. The simulation is run in arc-length control.   shows how the fracture develops at different stages of the analysis. The first crack develops in the longer beam span, from the lower flange in tension towards the web across the patch interface. This failure of the structure corresponds to the first drop in the load-displacement curve presented in Fig. 35. After this point, the beam is still able to carry some load until a second crack develops above the middle support, this time nucleating at the top flange, which is in tension at this location. Fig. 36 shows a close-up view of the adaptively refined mesh around the first crack, which emphasizes the importance of adaptive refinement when applying the phase-field approach to larger structures, where the fractured area may cover only a very small portion of the total domain. In this simulation, the initial mesh includes 5112 elements and the final one after local refinement around the cracks counts 2.92 × 10 4 elements. For comparison, a uniformly pre-refined mesh with the same minimum element size would instead require 1.31 × 10 6 elements.

Conclusion
We presented an efficient simulation framework to apply the phase-field approach to brittle fracture to complex shell structures. Our framework is based on isogeometric analysis, using a rotation-free Kirchhoff-Love shell formulation for structural analysis and both second-and fourth-order phase-field formulations for fracture.
For the application to multipatch structures, we proposed penalty formulations which can handle arbitrary patch connections with matching or non-matching interfaces. The various penalty terms involved are scaled with the corresponding problem parameters such that they all can be controlled by one global penalty parameter, which is set to α = 10 3 independently of the problem parameters. Moreover, the structural penalty terms are degraded with the phase-field degradation function in order to ensure a consistent scaling of structural stiffness and penalty stiffness in the fractured zones.
Furthermore, we presented a multistep predictor-corrector algorithm for adaptive refinement based on LR NURBS, using the phase-field value as refinement indicator. We proposed a local Bézier projection for mapping  the history field, which is defined at integration points, from coarser to finer meshes. The corresponding projection matrix depends only on the polynomial degree and is identical for all elements of a mesh throughout all load steps, hence it can be precomputed at the beginning of the analysis, making this mapping approach computationally inexpensive. We also demonstrated that the proposed approach can accurately predict and refine the crack path even in cases where fracture fully develops in a single load step, yielding identical results to simulations with static pre-refined meshes. This is crucial for the application to large structures without a priori knowledge of the crack path, since uniform refinement with an element size resolving the internal length scale would be prohibitive due to the excessively high number of elements required.
Furthermore, we investigated and compared the computational efficiency of loosely vs. strongly coupled staggering strategies and of second-vs. fourth-order phase-field formulations within our framework. From our simulations, we found that a strongly staggered scheme is generally more efficient since it requires in total fewer solutions of the computationally more expensive momentum equations (the structural problem is intrinsically nonlinear due to the tension-compression split of the strain tensor). Moreover, we found that the fourth-order phasefield model is in general more efficient than the second-order model, although the formulation is more involved, due to different reasons. Firstly, the faster convergence of the method leads to fewer staggered iterations of the coupled problem, i.e., fewer solutions of the momentum equations. Secondly, it has a narrower crack profile which decreases the necessary number of elements to resolve the crack. Finally, the more accurate crack prediction within the proposed predictor-corrector scheme for adaptive refinement leads to fewer recomputation steps compared to the second-order formulation. Therefore, we chose to use the fourth-order model whenever possible, i.e., when the domain is smooth. For non-smooth multipatch structures, however, we resort to the second-order formulation, since the orientation of the directional derivatives at the patch interfaces is not always unique in such cases.
For future research, we plan to extend the presented approach to ductile fracture, envisaging the application to steel structures, as well as the extension to dynamic fracture and fatigue.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.