A non-intrusive global/local approach applied to phase-field modeling of brittle fracture

This paper aims at investigating the adoption of non-intrusive global/local approaches while modeling fracture by means of the phase-field framework. A successful extension of the non-intrusive global/local approach to this setting would pave the way for a wide adoption of phase-field modeling of fracture, already well established in the research community, within legacy codes for industrial applications. Due to the extreme difference in stiffness between the global counterpart of the zone to be analized locally and its actual response when undergoing extensive cracking, the main foreseen issues are robustness, accuracy and efficiency of the fixed point iterative algorithm which is at the core of the method. These issues are tackled in this paper. We investigate the convergence performance when using the native global/local algorithm and show that the obtained results are identical to the reference phase-field solution. We also equip the global/local solution update procedure with relaxation/acceleration techniques such as Aitken’s \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta ^2$$\end{document}Δ2-method, the Symmetric Rank One and Broyden’s methods and show that the iterative convergence can be improved significantly. Results indicate that Aitken’s \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta ^2$$\end{document}Δ2-method is probably the most convenient choice for the implementation of the approach within legacy codes, as this method needs only tools already available for the so-called sub-modeling approach, a strategy routinely used in industrial contexts.


Introduction
The variational approach to fracture by Francfort and Marigo [1] and the related regularized formulation of Bourdin et al. [2][3][4][5], commonly referred to as phase-field model of (brittle) fracture, is a widely accepted framework for modeling and computing fracture phenomena in elastic solids. The phase-field framework for modeling systems with sharp interfaces consists in incorporating a continuous field variable-the so-called order parameter-which differentiates between multiple physical phases within a given system through a smooth transition. In the context of fracture, such an order parameter (termed the crack phase-field) describes the smooth transition between the fully broken and intact material phases, thus approximating the sharp crack discontinuity, as sketched in Fig. 1. Gerasimov et al. Adv. Model. and Simul. in Eng. Sci. (2018) 5:14 Page 2 of 30 The evolution of this field as a result of the external loading conditions models the fracture process. The formulation is strongly non-linear and calls for the resolution of small length scales.
The use of phase-field approaches in the case of structures of industrial complexity has been the subject of limited investigations thus far and poses a number of challenges. In this paper, in order to move forward in this direction we advocate the use of non-intrusive global/local strategies initially proposed in [6]. When dealing with large structures, fracture phenomena most often occur in regions of limited extent only. Moreover, in the case of brittle fracture most of the structure behaves elastically. These features are particularly appealing for global/local approaches as they make it possible to first compute the global model elastically, and then determine the critical areas to be re-analyzed, while storing the factorization of the decomposition of the structural stiffness. The local models are then iteratively substituted within the unchanged global one, which has the advantage of avoiding the reconstruction of the mesh of the whole structure. In fact, this is the main motivation of non-intrusive global/local approaches: to avoid the modification of the finite element model used by engineers, the creation of a complex global model being by far the most time-consuming task, a task which is more and more externalized.
The phase-field simulation of fracture processes with legacy codes bears a number of advantages which fit perfectly within the framework of non-intrusive coupling approaches using pre-defined 'fixed' meshes. The most obvious advantage is the ability to track automatically a cracking process by the evolution of the smooth crack field on a 'fixed' mesh which, in the proposed procedure, is the mesh of the local model. This is a significant advantage over the discrete fracture description, whose numerical implementation requires explicit (in the classical finite element method, FEM) or implicit (within XFEM) handling of the discontinuities. The possibility to avoid the tedious task of tracking complicated crack surfaces in 3D significantly simplifies the implementation. The second advantage is the ability to simulate complicated processes, including crack initiation (also in the absence of a crack tip singularity), propagation, coalescence and branching without the need for additional ad-hoc criteria and with very few parameters to be identified. This feature is particularly attractive for industrial applications, as it minimizes the need for time-consuming and expensive calibration tests.
Due to the extreme difference in stiffness between the global counterpart of the zone to be re-analyzed locally and its actual response when undergoing extensive cracking, the foreseen fundamental issues associated with the use of the global/local strategy in combination with phase-field fracture modeling are robustness, accuracy and efficiency of the fixed point iterative algorithm which is at the core of the method. Also, the finite element treatment of the phase-field formulation of brittle fracture is known to be computationally demanding, mainly due to the non-convexity of the energy functional to be minimized with respect to both arguments (the displacement and the phase field) simultaneously [45][46][47]. As a result, the so-called monolithic approach manifests major iterative convergence issues of the Newton-Raphson procedure. A new line-search scheme [46] and modified Newton methods [47] have been recently proposed to tackle this problem. Alternatively, staggered (also termed partitioned, or alternate minimization) solution scheme is widely used. This is based on decoupling of the strongly non-linear weak formulation into a system and then iterating between the equations [2][3][4][5][7][8][9][11][12][13]15,16]. The staggered scheme is proved to be robust, but typically has a very slow convergence behavior of the iterative solution process, see e.g. [15,46,48]. In view of the above, a central question that arises when combining non-intrusive global/local approaches with phase-field modeling of fracture is how additional global/local iterations affect and possibly deteriorate the highly sensitive iterative behavior of the staggered scheme used to solve the phase-field equations. In this paper, we make a first attempt to address these questions.
The paper is organized as follows. In "The phase-field approach to brittle fracture" section, we outline the main concepts of phase-field modeling of brittle fracture and illustrate the specific formulation used in the present paper. "Global/local approach in a non-intrusive setting" section introduces the non-intrusive global/local approach for the solution of the reference phase-field model considered in "The phase-field approach to brittle fracture" section. This is done in several steps. We start by illustrating an intrusive global/local scheme through a domain decomposition formulation in a variational setting well adapted to the phase field formulation. Several options are considered, including the so-called primal, dual and localized Lagrange multipliers based versions. This domain decomposition framework is used afterwards to define some convergence indicators in terms both of incompatibility of the reaction forces and of displacement jumps at the interface between the unchanged global model and the re-analyzed local one. The motivation here is to see which indicator or combination of indicators are the most suited to an appropriate estimation of the quality of the global/local iteration results with respect to the phase-field determination. The third version is then extended to the global/local setting, for which a non-intrusive computational procedure is devised. The numerical results which illustrate the performance of the proposed non-intrusive global/local approach as well as their qualitative and quantitative comparison with the reference solution are reported in "Results and discussion" section. Therein, we also outline and apply three relaxation/acceleration techniques, which are incorporated into the global/local iterative procedure and aim at improving its efficiency. Conclusions and outlook finalize the paper.

The phase-field approach to brittle fracture
In this section, we consider a mechanical system undergoing a brittle fracture process modeled with the phase-field formulation, and term this the reference problem. For this problem, we develop in "Global/local approach in a non-intrusive setting" section a global/local formulation, which is dissected numerically in "Results and discussion" section.
Let ⊂ R m , m = 2 or 3 be an open and bounded domain representing the configuration of a m-dimensional linear elastic body, and let D,0 , D,1 and N,1 be the (non-overlapping) portions of the boundary ∂ of on which homogeneous Dirichlet, non-homogeneous Dirichlet and Neumann boundary conditions are prescribed, respectively. In the following, we consider a quasi-static loading process with the discrete pseudo-time step parameter l = 0, 1, ..., such that the displacementū l and tractiont l loading data are prescribed on the corresponding parts of the boundary, see Fig. 2a.
For the mechanical system at hand, the phase-field formulation of brittle fracture [5] in an incremental variational setting relies on the following energy functional with Fig. 2 a Sketch of geometry and loading setup; b the computed crack phase-field evolution and the related minimization problem at each l ≥ 0. In the above, the displacement field u : → R m and the crack phase-field d : → [0, 1] are the arguments of E. As already mentioned, the limiting values of d, namely, d = 0 and d = 1 represent the undamaged and fully broken material phases. Furthermore, + and − are the so-called 'tensile' and 'compressive' parts of an additive decomposition of the elastic strain energy density function (ε) := 1 2 ε : C : ε = 1 2 λtr 2 (ε) + μtr(ε · ε), where, in turn, ε is the second-order infinitesimal strain tensor, C is the fourth-order elasticity tensor, and λ and μ are the Lamé constants. The decomposition of into + and − is required in order to distinguish between fracture behavior in tension and compression, more precisely, to avoid crack growth and crack faces interpenetration in compression. Here we use the spectral-based split, proposed in [8,9]: where a ± := 1 2 (a ± |a|) and ε ± := 3 I=1 ε I ± n I ⊗ n I with {ε I } 3 I=1 and {n I } 3 I=1 as the principal strains and principal strain directions, respectively. Finally, G c is the material fracture toughness, and 0 < diam( ) is the regularization parameter that controls the width of the transition zone of d between the two material states.
With E defined by (1), the state of the system at a given loading step l ≥ 0 is then represented by the solution of  Figure 2b depicts an example of phase-field pattern resulting from (4). Note that due to the d l−1 ≤ d requirement, problem (4) is a constrained minimization problem and its necessary optimality condition which enables computing the solution (u, d) ∈ Vū l × D d l−1 is a variational inequality. Its partitioned form reads as see e.g. [11,48], where E u and E d are the directional derivatives of the energy functional with respect to u and d, respectively. It is The displacement test space in (5) is defined as V 0 := {v ∈ H 1 ( ) : v = 0 on D,0 ∪ D,1 }. In (6), the components ∂ ± ∂ε are the corresponding 'tensile' and 'compressive' stresses, which are strongly non-linear in ε. In the case of the spectral-based split in (3), we obtain The related counterparts of the standard fourth-order elasticity tensor C read in this case where H + is the standard Heaviside function and H − := 1 − H + , J is the fourth-order symmetric identity tensor, whereas P ± are the fourth-order tensors obtained by differentiation of ε ± with respect to ε. Stemming from the irreversibility constraint d l−1 ≤ d the variational inequality E d ≥ 0 in (5) requires special solution algorithms, see e.g. [49,50]. Here, the irreversibility of d is enforced 'indirectly' via the notion of a history variable, as proposed in Miehe et al. [9]. The idea is that the tensile energy + can be viewed as the driving force of the phase-field evolution. Hence, the maximal + accumulated within the loading history and denoted as H l (x) := max ∀l + (ε(u)) can be used to prevent a decrease of the phase-field. H l substitutes the corresponding + term in the original E d , thus yielding (10) and the system for computing the solution (u, d) ∈ Vū l × H 1 ( ) is where E u is given by (6). We obtain in (11) an equality and unconstrained spaces for d and w.
The staggered solution algorithm for the system in (11) implies alternately fixing u and d, and solving the corresponding equations until convergence. The algorithm is sketched in Table 1.
Note that the equation E u = 0 in Table 1 is strongly non-linear due to the non-linearity of σ(u, d) , see equation (8). Therefore, at every staggered iteration k ≥ 1 with given d k−1 , a Newton-Raphson procedure is needed to compute u k , with e.g. u k−1 being taken as the initial guess, and TOL NR as a user-defined tolerance. Owing to the 'nested in' nature of the Newton-Raphson process, it has to be TOL NR < TOL Stag . In the presented numerical examples we take TOL NR := 10 −8 < TOL Stag := 10 −5 .

Global/local approach in a non-intrusive setting
The starting point towards a non-intrusive global/local approach to the phase-field problem (4) with E defined by (1) is a standard non-overlapping domain decomposition procedure applied to E. The resulting formulation is then extended to a global/local one in the spirit of [39,51], for which the non-intrusive computational scheme is devised.

Domain decomposition formulation
Let L be an open sub-domain of , where cracking (in a general setting: a strong localization effect due to non-linearity) is expected to take place, and let C ⊂ be its open complement ( C := / L ), where the material remains intact and elastic (in a general setting: non-linearity is negligible). In the following, the subscripts L and C always stand for local and complementary, respectively. It is typical to assume that L represents a reasonably small 'fraction' of such that | L | | C |. Let also ⊂ be the interface between L and C , a set with one dimension less than the dimension of , such that With an application to the problem sketched in Fig. 2, this domain decomposition idea is presented in Fig. 3a. We now introduce two functions on L and C , namely, u L ∈ H 1 ( L ) and u C ∈ and assume that the displacement u ∈ Vū l stemming from the solution of problem (4), can be represented as We furthermore introduce a function d L : L → [0, 1] such that the phase-field d ∈ D d l−1 stemming from the solution of (4) has the representation Using (13) and (14) in the function W in (2), we arrive at the energy functionals in the corresponding sub-domains, namely, and such that, also owing to (12), it holds where, to recall, E is the original reference functional (1). As a result, the domain decomposition variational formulation, which is equivalent to reference formulation (4), reads arg min The advantage of 'replacing' (4) with (23) is that one of the two sub-problems stemming from (18), more precisely, the complementary one, will be linear: indeed, W (ε(u), 0) = (ε(u)), thus yielding the standard linear stress-strain relation σ(u) := ∂W ∂ε (ε(u), 0) = C : ε(u). And this, moreover, will take place in a 'large portion' of , since by assumption Due to the strong displacement continuity requirement (12), formulation (18) is called primal in the literature, see e.g. [52]. This requirement may be too restrictive from the computational standpoint [53]. Relaxing, or rather neglecting (12), results in the appearance of the traction-like terms in the corresponding sub-domain energy functionals (15) and (16): and with λ C , λ L ∈ L 2 ( ) being the (unknown) Lagrange multipliers, which represent tractions. In this case, however, the 'argmin max'-problem being posed for is under-determined, since no relation is yet specified between u L and u C , nor between λ L and λ C .
Two standard ways to proceed with (19) and (20), and obtaining a variational formulation equivalent to the original one in (4) are as follows.
Option 1 One imposes a strong continuity between λ C and λ L on , by setting in E 1 and Summing the obtained functionals leads to Note that, in this case, , since the surface integral over the interface provides the weak continuity between the local and complementary displacement fields. λ is the (unknown) Lagrange multiplier. The corresponding variational problem for E which approximates the reference problem (4) then reads Owing to condition (21), formulation (23) is called dual in the literature, see e.g. [54]. The relation between the solution (u, d) of the reference problem (4) and the solution triple (u C , u L , d L ) is given by (13) and (14). Option 2 One preserves the representations (19) and (20), and, in contrast to (21), imposes only a weak continuity between λ C and λ L on . The latter is achieved by introducing the functional with u ∈ H 1 ( ) representing the (unknown) Lagrange multiplier, which has the dimension of a displacement. Summing E 1 and E 2 with E 3 , and also regrouping the terms, we finally obtain From (25), it can be grasped that the introduction of u enables also to implicitly provide a weak continuity between u L and u C across via an intermediate function u (this is in addition to the already incorporated weak continuity between λ L and λ C ). Conclud- (25) which approximates the reference problem (4) will be as follows: In this case, the representation of u stemming from the solution of problem (4) in terms of the solution triple (u C , u L , u ) stemming from (26) reads as whereas the representation for d in terms of d L defined by (14) remains unaltered. In the literature, formulation (26) is sometimes called the localized Lagrange multipliers based formulation (we abbreviate this as LLM), where the term 'localized' is used to associate the multipliers λ C , λ L and u with the corresponding sub-domains, see e.g. [55][56][57]. Table 2 briefly summarizes the considered formulations. Formulation (23) is seemingly less computationally demanding than (26), since there is only one extra field λ to be solved for in the former case, versus the triple (u , λ C , λ L ) of unknown fields in the latter one. The potential advantage of (26) over (23) is a greater flexibility, at the finite element discretization stage, of handling the interface between complementary and local domains.
As follows, we move on with the LLM formulation (26) and extend it to the global/local setting, for which, in turn, a non-intrusive solution procedure is devised. This will lead to a non-intrusive global/local approach to the phase-field formulation (4).

Global/local formulation
As a first step, a so-called fictitious domain F is introduced to 'fill the gap' obtained in by removing L from it, see Fig. 3b. It is assumed that F is constituted by a material with the same linear elastic behaviour as in C . It is also assumed that F is open (i.e. ⊂ F ). Unification of F with and C forms the global domain G , that is, The fictitious domain F is furthermore assumed free of geometrical 'imperfections' which may be present in L , see Fig. 3b. Therefore, it is in general G = , and the constructed global domain G should not be confused with the original reference domain .
Summing up the above, the role of the fictitious domain F is twofold: it replaces the "sub-regions" of a structure (reference domain) containing geometric details (e.g. holes,

Formulation
Imposed continuity between Unknowns inclusions etc.) and/or constitutive non-linearity by there details-free and linearly elastic "counterparts". The obtained global domain G is then straightforwardly suitable for meshing and solving procedures within legacy codes. As it will be also seen below, the use of F is essential to realize the concept of non-intrusiveness of the computational scheme for solving the coupled global/local formulation. Next to this, it is assumed that there exists a continuous prolongation of u C into F . That is, we introduce a function u G ∈ H 1 ( G ) such that u G | C ≡ u C and u G = u C on in the sense of trace. The former also implies that u G = 0 on D,0 and u G =ū l on D,1 .
Owing to the definitions of G and u G , the first term in (25) is recast as follows and we also substitute u G for u C in the third and fourth integrals in (25). This yields the desired global/local representation (approximation) of the reference energy functional E in (1), namely, (We used the a to distinguish between the previously considered E and the constructed E.) The resulting global/local variational problem, which approximates the reference formulation (4) reads and the relation between the solution u of (4) and the solution triple (u G , u L , u ) of (29) is given by In what follows, for the sake of compactness we set (u G , u L , d L , u , λ C , λ L ) =: z.

Coupled system in weak form
To present the weak formulation of (29), we introduce the directional derivatives of E with respect to the various components of z. Recalling the function W defined by (2), the 'main' three derivatives read where where w L ∈ H 1 ( L ) is the test function. The following 'modified' version of E d L , adjusted to account for the irreversibility of the phase-field evolution, will be used in our computations: This is similar to the modification discussed for equation (10). The remaining three variational derivatives of E are where v ∈ H 1 ( ) and β C , β L ∈ L 2 ( ) are the corresponding test functions. Using equations (30) and (31), (32), the global and local weak problems are, respectively, formed: whereas equations (33), (34) and (35) are used for establishing the (weak) coupling between them: Here, (u G , u L , d L , u , λ C , λ L ) is the vector of the unknowns to be solved for, and (v G , v L , w L , v , β C , β L ) is the vector of the corresponding test functions. For the presented system of equations, a computational scheme can already be devised. We should notice, however, that equation (G) in the current form does not fit in the notion of non-intrusiveness yet. Indeed, being a linear one, it can naturally be solved for u G 'straightforwardly'. But the presence of the two domain integrals, namely, over G and F ⊂ G would imply in this case the need to simultaneously access the corresponding stiffness matrices (in the following, K G and K F ), or, in other words, a necessity of modifying K G -a situation that contradicts the concept of non-intrusiveness. Avoiding this can be done in two steps: first, by introducing a partitioning of equation (G), and then, devising the appropriate iterative solution procedure. The former will be presented here, and the latter is addressed in "Non-intrusive computational scheme" section.
We focus on the domain integral over F in (G). The divergence theorem leads to where n ∂ F is the unit outward normal vector to ∂ F . The first term in the right-hand side of (36) can be canceled using the divergence-free assumption for the stress (no body forces in F ). The second term can be simplified as follows. In the most general case, ∂ F is composed of five non-overlapping parts, including . More precisely, as sketched in Fig. 4a, and hence Fig. 4 a The possible complex nature of ∂ F illustrating equation (37); b choice of L that results in In this case, due to the following basic properties • σ(u G ) · n N,0 = 0 on N,0 , • v G = 0 on D,0 and on D,1 , • σ(u G ) · n N,1 =t l on N,1 , the corresponding integrals in the right-hand side of (38) are simplified, thus yielding Here, n is the unit normal vector on , outward with respect to F . To further simplify the last expression, we note that it is always possible to pick L (and, hence, the resulting F ) such that ∂ F ∩ N,1 = ∅, see Fig. 4b, and, as a result, the last surface integral cancels. For (36), this eventually yields: It can now be assumed that, given u G , there exists λ F ∈ L 2 ( ) such that holds. In the above, we use the subscript F to indicate, according to (39), the relation of the corresponding quantity to F (more precisely, to the restriction of u G to F ). Owing to (39) and (40), we finally arrive at the following partitioned representation of equation (G): Equations (G 1 ), (G 2 ), system (L) and coupling equations (C 1 ), (C 2 ), (C 3 ) constitute what we term global/local coupled system, which is to be solved for the vector (u G , u L , d L , u , λ C , λ L ).

Non-intrusive computational scheme
Let n ≥ 0 be the iteration index. For designing at a fixed loading step l the iterative solution procedure for the global/local system defined by (G 1 ), (G 2 ), (L), and (C 1 ), (C 2 ), (C 3 ), the following prerequisites are taken into account: (a) Since the data (ū l ,t l ) are posed on D,1 , N,1 ⊂ ∂ G , the process initialization (i.e. iteration n = 0) is started with the solution of global problem (G 1 ), (G 2 ).
(b) In order to fit equation (G 1 ) with λ F = λ F (u G ) in the concept of non-intrusiveness, λ F must be treated as a known quantity. This defines the order in which equations (G 1 ) and (G 2 ) are solved at any iteration n ≥ 0: the solution of (G 2 ) precedes the solution of (G 1 ). In this case, as desired, the stiffness matrix K G remains unaltered; the access to K F is still required, but only at the stage of solving (G 2 ), not (G 1 ). (c) For solving (G 1 ), λ C must be also known. At n = 0, λ C can simply be taken from the previous loading step. At n ≥ 1, we use coupling equation (C 1 ) for the extraction of λ C , assuming λ L is already known. This defines the order in which the global and local problems are solved: at any iteration starting from n = 1, the solution of (L) precedes the solution of (G 1 ). We also notice that: (d) Coupling equation (C 3 ) provides the boundary condition for u L of the local problem (L). (e) Coupling equation (C 2 ) is used for the recovery of u .
As follows from (c) and (e), elimination of λ C and u from the set of unknowns to be originally solved for is achieved. These two quantities, as well as λ F , are the recovered ones.
The summary of the solution operations to be performed at any iteration n of the procedure, excluding the initialization step (n = 0), is as follows: • solution of local problem (L) coupled with (C 3 ), • recovery phase using (C 1 ) and (G 2 ), • solution of global problem (G 1 ), • recovery phase using (C 2 ).
The detailed scheme, including the iteration n = 0, is depicted in Table 3. Note that in all equations in the table we omit dx and ds.

Staggered process for the local problem
Solution of the local system in Table 3 at the given global/local iteration n ≥ 1 requires an additional nested iterative solution process. In our case, this is the staggered procedure from Table 1, which is adjusted to handle an extra variable λ L , and is also equipped with the appropriate definition of the input (initial guess) data and of the stopping criterion.
The initial guess for the staggered loop (with the iteration index k ≥ 0) is chosen as follows. At iteration n = 1 (and staggered iteration k = 0) the values (u L,l−1 , d L,l−1 , λ L,l−1 ) known from the previous loading step are used as the initial guess. At n ≥ 2 (and staggered iteration k = 0), we naturally take (u n−1 L , d n−1 L , λ n−1 L ). At any fixed iteration n ≥ 1, the accuracy check for the solution triple (u (k) L ) obtained at the staggered iteration k ≥ 0 is performed as follows: where, to recall, E u L is given by (31). If (41) is fulfilled-note again that in the following numerical test, TOL Stag := 10 −5 -the staggered process is stopped, we set (u (k) L ) =: (u n L , d n L , λ n L ) and perform n + 1 → n.

Table 3 Non-intrusive iterative solution process for (G 1 ), (G 2 ), (L), and (C 1 ), (C 2 ), (C 3 ) at a fixed loading step l
Initialization, n = 0:   Table 3 is rather straightforward. Indeed, at any iteration n ≥ 1, the solution outcome is denoted as (u n G , u n L , d n L , u n , λ n C , λ n L ). Plugging this in equations (G 1 ), (G 2 ), (L), (C 1 ), (C 2 ), (C 3 ) and comparing the obtained outcome with the corresponding equations in Table 3, it is straightforward to locate the imbalanced quantities: and Therefore, the solution accuracy at n is measured by the quantity (44) where Res Stag is the staggered residual of the local problem with 'last k' denoting the index of the converged staggered solution (see "Staggered process for the local problem" section for details), and λ n+1 F is recovered (post-processed) from the right-hand side of (43). The stopping criterion for the global/local loop can then be defined as with TOL GL to be prescribed. Owing to the 'nested in' nature of the staggered process, it has to be TOL Stag < TOL GL . Recalling that TOL Stag = 10 −5 , we set TOL GL := 10 −4 . In our computations (see Table 3), we use a more convenient form of the stopping criterion. Setting we define η n := (η n u ) 2 + (η n λ ) 2 , and use this quantity to now check η n ≤ TOL GL := 10 −6 .
This choice of TOL GL fulfills the requirement η n ! < TOL GL − TOL Stag , which is stipulated by (44), (45), and the already prescribed above magnitudes of TOL GL and TOL Stag . Since the quantity η naturally stems from the global/local solution accuracy check E z (z n ; y) = 0, it represents not only the iterative convergence indicator, but also the solution accuracy indicator-a very desired property, since the former is only suitable for tracing the convergence of the corresponding iterative solution process, but, clearly, is not adequate for stopping criterion. The corresponding ingredients η u and η λ are only iterative convergence indicators, but none of them provides an adequate check of the solution accuracy. In particular, since η u measures, though implicitly, the displacement continuity-a match between u G and u L across (recall that the traction continuity-a match between λ C and λ L on -is, in our case, fulfilled automatically), it is also the indicator of a good "gluing" between the two models.

Incremental setting
For later developments ("Results and discussion" section), it proves convenient to reformulate the global equation in incremental form. It is straightforward to see that for a given global/local iteration n ≥ 1, this reads: given the triple (u n−1 G , λ n−1 F , λ n−1 C ) known from the iteration n − 1, as well as (λ n F , λ n C ) 'recovered' at the iteration n, we solve for u G =: ( u G ) n and set We term equation (46) a 'direct update' within the global/local iterative procedure. This is in contrast to the notion of a 'relaxed/accelerated update' to be considered in Section .

Finite element discretization
In the following, for the sake of simplicity, we assume the dimension of the reference problem is 2. Let P be a finite element partition of into triangles or quadrilaterals, I be the number of nodes in P, and N i , i = 1, . . . , I be the nodal shape function associated with the node i and supported on the collection of elements in P that share i. Finally, let a scalar-valued quantity· i represent the nodal value. The standard discretization of the solution of the reference problem (using Voigt's notation) is as follows: where N u is the 2 × 2I interpolation/basis matrix and u is the 2I × 1 displacement nodalvector, and where N d is the 1 × I interpolation/basis matrix and d is the I × 1 phase-field nodalvector. N i constitutes both matrices N u and N d . The corresponding representation of the gradients reads: where B u is a 3 × 2I matrix, and ∇d = ∂d ∂x ∂d ∂y where B d is a 2 × I matrix. The problem test functions v and w, as well as their gradients are discretized accordingly.
For the global/local formulation, we assume the existence of the partitions P G and P L of G and L , respectively. Using this and adopting the above notations, the solution discretizations are given by such that Note that the superscripts G and L are introduced in the corresponding definitions of N u , N d and B u , B d .
To construct the discretization of the Lagrange multipliers λ C , λ L , u and the supplementary quantity λ F on we account for the following. In the most general situation, three distinct partitions of may be assumed: the restrictions of P G and P L -denoted as T G and T L , respectively-which serve to create the corresponding bases for λ C (also λ F ) and λ L , and the 'independent' partition T to be used for creating the basis for u . Introducing the three related basis matrices N G λ , N L λ and N u , we write In the following, for our numerical example, we assume that: • the partitions T G , T L and T match (this is usually termed a 'matching case'); • the basis in the global and local domains is identical, that is, N G u = N L u =: N u ; • the basis on the interface is obtained from N u by the corresponding restriction, that is, N G λ = N L λ = N u = N u | ; • the nodal shape functions N i composing bases N u and N d are piecewise linear.
Because of the matching interface situation, no intricate data transfer methods (the construction of prolongation and restriction operators, generalized inverse matrices etc.) are required. Also, with the above choice of the discretization basis for the Lagrange multipliers, the related inf-sup condition is fulfilled, see e.g. [58,59].
Using expressions (47), (48) and (49) along with the above assumptions, the matrix representation of all equations in Table 3 is straightforward.

Results and discussion
To illustrate the proposed approach, we consider the following benchmark problem. A square specimen with two holes of different diameters is subjected to tension loading, see Fig. 5a. The holes are introduced to weaken the structure and to facilitate the specimen cracking in absence of a stronger singularity such as a pre-existing crack. The holes location is chosen such that prediction of the sub-region where cracking occurs (hence, the local domain for the forthcoming global/local analysis) is feasible. Taking a different size of the holes is intended to obtain a geometrically non-trivial crack pattern, as depicted in Fig. 5b. This, moreover, results in a multi-stage crack propagation process to be manifested by a load-displacement response with two peak points, see Fig. 5c for a sketch, and Figs. 7 and 12 for the actually obtained results. We believe that the present setup, being neither extremely complex, nor trivial, is suitable for the purpose of a qualitative and quantitative comparison between the reference results and results obtained with the proposed global/local approach.
The geometric data are as follows (all given in mm): a = 1, b 1 = 0.197, b 2 = 0.210, b 3 = 0.490 with the hole diameters c 1 = 0.247 and c 2 = 0.0806. The material data are: Young's modulus E = 210 GPa, Poisson's ratio ν = 0.3 and the critical energy release rate G c = 2.7 · 10 −3 kN/mm. The characteristic length in the phase-field formulation is = 1.5 · 10 −2 mm. We consider the plane-strain situation.
We recall that we use P 1 -triangles for approximating all unknown variables both in the reference and global/local formulations. The minimum finite element size in the reference and local domains is 0.004 mm, the maximum element size in the reference and global domains is 0.1 √ 2 mm. The former fulfills the heuristic requirement h < /2 for the element size inside the localization zone (i.e. the support) of d. The reference domain partition contains 18,672 elements. The discretizations of the global and local domains contain 200 and 18,552 elements, respectively. That is, in our case, the reference and global/local problems have a comparable discretization size, as can be grasped from Fig. 6.

Reference and global/local results
We start here with the presentation of the quantitative and qualitative reference and global/local results and their comparison. As desired, the two load-displacement curves in Fig. 7 are identical in the entire range of loading, including the pre-and post-peak behavior. The computed phase-field profiles in Fig. 8 are also in a very good agreement. This is already a good indicator of the potential of the global/local approach with application to systems with strong non-linearity and localization.  For a better insight into the the iterative convergence behavior of the global/local solution process, in Fig. 9 we depict the convergence indicators from "Accuracy/convergence check" section for four given loading steps corresponding to the points 1-4 of our interest sketched in Fig. 5c. Thus, we plot the quantities η u , η λ and η = η 2 u + η 2 λ such that the amount of global/local iterations required for the solution convergence at the step (also in comparison with other steps) can be detected.
The first important observation is that η u , which implicitly measures the displacement discontinuity between the solutions of the global and local problem across the interface, is two orders of magnitude less than η λ . Thus, its contribution to η, which is used not only for tracing the convergence of the iterative solution process, but also for the solution accuracy check, is negligible. This means that a stopping criterion based solely on the use of η u (what seems typical for the global/local approaches in e.g. plasticity) will yield, in our case, erroneous results. Secondly, it can be noted that a quite large amount of global/local iterations is needed, especially at loading steps corresponding to the peak loads of the load-displacement curves in Fig. 7 (the points of interest 2 and 4 from Fig. 5c).
Resulting from the slow convergence of the global/local procedure, the corresponding cumulative computational time turns out to be high, see Fig. 10, where also the time for solving the reference formulation by the staggered scheme is depicted. For the given setup, with a standard machine (Intel(R) Core(TM) i7-3770 OK, CPU 3.5 GHz, RAM 16.0 GB) it takes about one hour of staggered computations vs. approximately four hours required for the global/local approach. (We should note however that our goal was not to gain computational efficiency, but rather to enable computations with legacy codes.) Fig. 9 Convergence behavior of the global/local iterative solution process at four different loading steps (the points 1-4 from Fig. 5c), illustrated in terms of the indicator η, as well as its ingredients η u and η λ High efforts are not surprising, as the global/local problem has a larger discretization size than the reference problem, and three nested iterative processes vs. two for the reference problem. The latter results in a larger time per loading step, as can be seen in Fig. 11.
It can be grasped that the rapid increase of cumulative time in Fig. 10 for both formulation appears at loading steps related to the peak points 2 and 4. Also, regardless of the formulation, the computational time per step in Fig. 11 at these points is sig-nificantly higher (by almost two orders of magnitude, to be more precise) than at the pre-peak loading steps. These observations correlate with the convergence results from Fig. 9.
Non-convexity and non-linearity of the global/local formulation, as well as the complicated multi-level iterative nature of the related iterative solution procedure result in a generically slow convergence of the approach. Another impacting factor that should be noted is that the stiffness matrix of the global problem K G is never updated within the global/local computation process. Incorporation of an incremental update relaxation in this process is thus our next goal, with the objective to obtain an acceleration of the convergence process.

Relaxation/acceleration techniques: Aitken's, SR1, Broyden et al.
Following [39] and [51], we will consider and incorporate two types of relaxation/acceleration techniques into our approach: Aitken's 2 -method (also known as dynamic relaxation, whose efficient implementation in fluid-structure interaction computations has already been reported [60,61]) and Quasi-Newton correction. Within the family of Quasi-Newton correction formulae, we restrict ourselves to the Symmetric Rank One (SR1) and the Broyden update versions.
Technically, both types deal with the global solution update equation (46) and modify it specifically. Let us consider (46) written in terms of the nodal displacementŝ where, owing to (G incr ), one has with C as a 3 × 3-matrix representation of C, and Aitken's method modifies (50) at any iteration n ≥ 2 introducing the damping factor ω n−1 = f (ω n−2 , ( û G ) n−1 , ( û G ) n ) s.t. ω 0 = 1 as follows: whereas the Quasi-Newton correction modifies (50) at any iteration n ≥ 2 by replacing the matrix In (51), we explicitly have with ω 0 = 1. In (52), the SR1 update formula implies whereas the Broyden update reads as In both cases K 1 = K G . Further details about computing the inverse matrices, efficient data storage etc. can be found e.g. in [6,39]. Also, following Conn et al. [62], the SR1 update formula (54) is used only if |(r n ) T ( û G ) n−1 | ≥ c 1 |r n | |( û G ) n−1 |, with a constant c 1 ∈ (0, 1). Otherwise, we simply set K n := K n−1 . This helps preventing the convergence issue of the global/local procedure using the SR1 based relaxation. The results obtained with the relaxation/acceleration techniques are depicted in Figs. 12-15. As can be seen from Fig. 12, all three considered techniques yield identical load-displacement curves, also identical to the curve obtained from the global/local procedure with no relaxation/acceleration.
Similarly to Figs. 9 and 13 presents and compares the convergence of the global/local iterative procedure and its acceleration/relaxation versions at the four loading steps of interest. Here, we only plot the indicator η and not its ingredients. For a given point, the Fig. 12 Comparison of the load-displacement curves amount of iterations required for the convergence of the solution process in all acceleration/relaxation techniques is similar, but is less (in some cases, significantly) than in the original unaccelerated case. Figure 14 compares the phase-field solutions of the global/local formulations computed using the corresponding acceleration/relaxation techniques. It can be observed that even though the load-displacement curves are identical in all cases, the corresponding phase-field profiles are not. This can be explained, first of all, by the solution nonuniqueness of the original reference phase-field formulation, and, secondly, by the fact that the global/local formulation is only the approximation of the reference one.
From the time-displacement curves comparison in terms of both 'cumulative time' and 'time per loading step', Fig. 15, it can be concluded that the desired improvement of efficiency of the original procedure has indeed been achieved. However, in the global time scale, all three techniques have a very similar effect, at least for the considered example.

Conclusions
We combined the adoption of non-intrusive global/local approaches with phase-field modeling of brittle fracture, with the main objective to pave the way for a wide adoption of this framework for industrial applications within legacy codes. We investigate the convergence performance of the fixed-point scheme used for the global/local iterations and showed that the obtained results are identical to the reference phase-field solution. In order to accelerate the observably quite slow convergence behavior, especially close to and beyond the peak point(s) of the load-displacement response, we also equipped the global/local solution update procedure with relaxation/acceleration techniques such as Aitken's 2 -method, the Symmetric Rank One and Broyden's methods. Findings showed that the iterative convergence can be improved significantly, to a similar extent for all investigated methods. Aitken's 2 -method is probably the most convenient choice for the implementation of the approach within legacy codes, as this method needs only tools which are often used for the so-called sub-modeling strategy, which is well known and widely used in industrial contexts. Several extensions and improvements of the proposed framework are foreseen, such as the study of the effect of global modeling choices on the iterative convergence behavior, the investigation of alternative boundary conditions for the coupling, the implementation of mortar-type approaches to enable non-matching meshes at the boundary between local and complementary domains. Moreover, in practical applications when i.e. the evolving localization areas are not known á priori, the global/local approach must be supplied with obtained in Ambati et al. [15], p. 392, equation (34) suggests that + , or rather the history field H l , may serve as an "indicator" for the adaptive choice of L . In (56), V ⊂ is an averaging volume, ξ ∈ V , and g(ξ) is a bell-shaped function, typically a Gaussian, such that 1 |V | V g(ξ) dξ = 1 holds. These issues are all open for further research.