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

Advanced Modeling and Simulation in Engineering Sciences, 5 (1)

section_content: 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. 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. 

In the past decade, both phase-field and non-intrusive global/local approaches have been extended to deal with a growing number of situations of interest for engineers. The currently available phase-field formulations of brittle fracture encompass static and dynamic models. We mention the papers by Amor et al. [7], Miehe et al. [8, 9], Kuhn and Müller [10], Pham et al. [11], Borden et al. [12], Mesgarnejad et al. [13], Kuhn et al. [14], Ambati et al. [15], Wu et al. [16], where various formulations are developed and validated. Recently, the framework has been also extended to ductile (elasto-plastic) fracture [17] [18] [19] [20] [21] [22], pressurized fracture in elastic and porous media [23, 24], fracture in films [25] and shells [26] [27] [28], and multi-field fracture [29] [30] [31] [32] [33] [34] [35] [36]. Non-intrusive global/local approaches have also been applied to a quite large number of situations: the computation of the propagation of cracks in a sound model using the extended finite element method (XFEM) [37], the computation of assembly of plates introducing realistic non-linear 3D modeling of connectors [38], the extension to non-linear domain decomposition methods [39] and to explicit dynamics [40, 41] with an application to the prediction of delamination under impact using Abaqus [42]. Alternative strategies can be derived from the Partition of Unity Method [43, 44]. 

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.
section_name: Introduction

section_content: 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 traction tl 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

where

is the kinematically admissible displacement space with H 1 ( ) := [H 1 ( )] m and H 1 denoting the usual Sobolev space, and

is the admissible space for d with d l-1 being known from the previous step. The condition d l-1 ≤ d is used to enforce the irreversibility of the crack phase-field evolution. 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

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

and the system for computing the solution (u, d)

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,

, 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 :=
section_name: The phase-field approach to brittle fracture

section_name: Table 1 Staggered iterative solution process for (11) at a fixed loading step l

section_content: 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. We now introduce two functions on L and C, namely,
section_name: Global/local approach in a non-intrusive setting

section_content: and assume that the displacement u ∈ V ūl stemming from the solution of problem (4), can be represented as

We furthermore introduce a function

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, E(u, d) ≈ E(u C, u L, d L, λ), 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 λ

, the variational problem for E in (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).
section_name: Domain decomposition formulation

section_content: 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, 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.
section_name: Global/local formulation

section_content: 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

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

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 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):

with λ F satisfying

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 ).
section_name: Coupled system in weak form

section_content: 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, tl ) 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: 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.
section_name: Non-intrusive computational scheme

section_content: 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, d L, λ L, ) 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

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

Initialization, n = 0: l, u L, l, d L, l, u Γ, l, λ C, l, λ L, l ). * See "Staggered process for the local problem" section, ** see "Accuracy/convergence check" section
section_name: Staggered process for the local problem

section_content: Derivation of the convergence and stopping criteria for the global/local iterative solution process in Table 3 is rather straightforward. Indeed, at any iteration n ≥ 1, the solution outcome is denoted as

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

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 2, and use this quantity to now check

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.
section_name: Accuracy/convergence check

section_content: 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.
section_name: Incremental setting

section_content: 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

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

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,

• 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.
section_name: Finite element discretization

section_content: 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 Fig. 5 a Specimen geometry and loading conditions; sketches of (b) the fracture pattern and c the load-displacement curve with the points of interest 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): The algorithmic parameters are: the loading ūl = l u with l ∈ [1, 110) and the increment size u := 0. 06 • 10 -3 mm, the tolerance magnitudes are TOL NR := 10 -8, TOL Stag := 10 -5, and TOL GL := 10 -6. 

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.
section_name: Results and discussion

section_content: 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. ) 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 displacements 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.
section_name: Reference and global/local results

section_content: 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.
section_name: Conclusions

section_name: Acknowledgements

section_content: Not applicable.
section_name: Acknowledgements

section_name: Funding

section_content: TG was partially funded by the ERC through the Starting Researcher Grant INTERFACES, GA 279439. NN was partially funded by the DFG through the IRTG 1627 "Virtual materials and their validation".
section_name: Funding

section_name: Author details

section_content: 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

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

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 Authors' contributions TG elaborated the concepts proposed and discussed in the manuscript, contributed to the numerical studies and wrote the draft of this paper, NN performed the numerical implementation and experiments, OA contributed to writing the part of the introductory section, and LDL contributed to the numerical experiments and paper revisions. All authors read and approved the final manuscript.

section_content: The authors declare that they have no competing interests.
section_name: Competing interests

section_content: Not applicable.
section_name: Ethics approval and consent to participate

section_content: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
section_name: Publisher's Note