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 Δ 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 Δ 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-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.

Fig. 1

Phase-field description of fracture (sketchy): is the crack phase-field

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-22], pressurized fracture in elastic and porous media [23, 24], fracture in films [25] and shells [26-28], and multi-field fracture [29-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-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–5, 7–9, 11–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 , or 3 be an open and bounded domain representing the configuration of a m-dimensional linear elastic body, and let and 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 , such that the displacement and traction loading data are prescribed on the corresponding parts of the boundary, see Fig. 2a.

Fig. 2

a Sketch of geometry and loading setup; b the computed crack phase-field evolution

For the mechanical system at hand, the phase-field formulation of brittle fracture [5] in an incremental variational setting relies on the following energy functionalwithand the related minimization problem at each . In the above, the displacement field and the crack phase-field are the arguments of . As already mentioned, the limiting values of d, namely, and 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 , where, in turn, is the second-order infinitesimal strain tensor, 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 and with and as the principal strains and principal strain directions, respectively. Finally, is the material fracture toughness, and is the regularization parameter that controls the width of the transition zone of d between the two material states.

With defined by (1), the state of the system at a given loading step is then represented by the solution ofwhereis the kinematically admissible displacement space with and denoting the usual Sobolev space, andis the admissible space for d with being known from the previous step. The condition 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 requirement, problem (4) is a constrained minimization problem and its necessary optimality condition which enables computing the solution is a variational inequality. Its partitioned form reads assee e.g. [11, 48], where and are the directional derivatives of the energy functional with respect to 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 obtainThe related counterparts of the standard fourth-order elasticity tensor read in this casewhere is the standard Heaviside function and , is the fourth-order symmetric identity tensor, whereas are the fourth-order tensors obtained by differentiation of with respect to .

Stemming from the irreversibility constraint the variational inequality 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 can be used to prevent a decrease of the phase-field. substitutes the corresponding term in the original , thus yieldingand the system for computing the solution iswhere 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 and d, and solving the corresponding equations until convergence. The algorithm is sketched in Table 1.

Table 1

Staggered iterative solution process for (11) at a fixed loading step l

Note that the equation in Table 1 is strongly non-linear due to the non-linearity of , see equation (8). Therefore, at every staggered iteration with given , a Newton–Raphson procedure is needed to compute , with e.g. being taken as the initial guess, and as a user-defined tolerance. Owing to the ‘nested in’ nature of the Newton–Raphson process, it has to be . In the presented numerical examples we take .

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 defined by (1) is a standard non-overlapping domain decomposition procedure applied to . 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 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 be its open complement (), 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 represents a reasonably small ‘fraction’ of such that . Let also be the interface between and , 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.

Fig. 3

Domain decomposition procedure: a the classic one, when is decomposed into local and complementary sub-domains and , respectively, which do not overlap and are coupled by the interface ; b its extension to the global/local setting, where a fictitious domain is introduced to form the so-called global domain

We now introduce two functions on and , namely, and such thatand assume that the displacement stemming from the solution of problem (4), can be represented asWe furthermore introduce a function such that the phase-field stemming from the solution of (4) has the representationUsing (13) and (14) in the function W in (2), we arrive at the energy functionals in the corresponding sub-domains, namely,andsuch that, also owing to (12), it holdswhere, to recall, is the original reference functional (1). As a result, the domain decomposition variational formulation, which is equivalent to reference formulation (4), readsThe 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, , thus yielding the standard linear stress-strain relation . 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):andwith being the (unknown) Lagrange multipliers, which represent tractions. In this case, however, the ‘’-problem being posed foris under-determined, since no relation is yet specified between and , nor between and .

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 and on , by setting in and Summing the obtained functionals leads toNote 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 which approximates the reference problem (4) then readsOwing to condition (21), formulation (23) is called dual in the literature, see e.g. [54]. The relation between the solution of the reference problem (4) and the solution triple 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 and on . The latter is achieved by introducing the functionalwith representing the (unknown) Lagrange multiplier, which has the dimension of a displacement. Summing and with , and also regrouping the terms, we finally obtainFrom (25), it can be grasped that the introduction of enables also to implicitly provide a weak continuity between and across via an intermediate function (this is in addition to the already incorporated weak continuity between and ). Concluding that , the variational problem for in (25) which approximates the reference problem (4) will be as follows:In this case, the representation of stemming from the solution of problem (4) in terms of the solution triple stemming from (26) reads aswhereas the representation for d in terms of 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 and with the corresponding sub-domains, see e.g. [55-57].

Table 2 briefly summarizes the considered formulations.

Table 2

Domain decomposition formulations of the reference problem (4)

Formulation	Imposed continuity between		Unknowns
	u_C&u_L	λ_C&λ_L
Primal, (18)	Strong	–	(u_C,u_L,d_L)
Dual, (23)	Weak	Strong	(u_C,u_L,d_L,λ)
LLM, (26)	Weak	Weak	(u_C,u_L,d_L,u_Γ,λ_C,λ_L)

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 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 is introduced to ‘fill the gap’ obtained in by removing from it, see Fig. 3b. It is assumed that is constituted by a material with the same linear elastic behaviour as in . It is also assumed that is open (i.e. ). Unification of with and forms the global domain , that is, . The fictitious domain is furthermore assumed free of geometrical ‘imperfections’ which may be present in , see Fig. 3b. Therefore, it is in general , and the constructed global domain should not be confused with the original reference domain .

Summing up the above, the role of the fictitious domain 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 is then straightforwardly suitable for meshing and solving procedures within legacy codes. As it will be also seen below, the use of 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 into . That is, we introduce a function such that and on in the sense of trace. The former also implies that on and on .

Owing to the definitions of and , the first term in (25) is recast as followsand we also substitute for in the third and fourth integrals in (25). This yields the desired global/local representation (approximation) of the reference energy functional in (1), namely,(We used the to distinguish between the previously considered and the constructed .) The resulting global/local variational problem, which approximates the reference formulation (4) readsand the relation between the solution of (4) and the solution triple of (29) is given byIn what follows, for the sake of compactness we set .

Coupled system in weak form

To present the weak formulation of (29), we introduce the directional derivatives of with respect to the various components of . Recalling the function W defined by (2), the ‘main’ three derivatives readwhere , and is the test function;where , and is the test function;where is the test function. The following ‘modified’ version of , 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 are where and are the corresponding test functions.

Using equations (30) and (31), (32), the global and local weak problems are, respectively, formed:andwhereas equations (33), (34) and (35) are used for establishing the (weak) coupling between them: Here, is the vector of the unknowns to be solved for, and 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 ‘straightforwardly’. But the presence of the two domain integrals, namely, over and would imply in this case the need to simultaneously access the corresponding stiffness matrices (in the following, and ), or, in other words, a necessity of modifying —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 in (G). The divergence theorem leads towhere is the unit outward normal vector to . 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 ). The second term can be simplified as follows. In the most general case, 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 illustrating equation (37); b choice of that results in for

with . In this case, due to the following basic propertiesthe corresponding integrals in the right-hand side of (38) are simplified, thus yieldingHere, is the unit normal vector on , outward with respect to . To further simplify the last expression, we note that it is always possible to pick (and, hence, the resulting ) such that , see Fig. 4b, and, as a result, the last surface integral cancels. For (36), this eventually yields:It can now be assumed that, given , there exists such thatholds. In the above, we use the subscript F to indicate, according to (39), the relation of the corresponding quantity to (more precisely, to the restriction of to ).

on ,

on and on ,

on ,

Owing to (39) and (40), we finally arrive at the following partitioned representation of equation (G): with satisfying

Equations (), (), system (L) and coupling equations (), (), () constitute what we term global/local coupled system, which is to be solved for the vector .

Non-intrusive computational scheme

Let be the iteration index. For designing at a fixed loading step l the iterative solution procedure for the global/local system defined by (), (), (L), and (), (), (), the following prerequisites are taken into account:As follows from (c) and (e), elimination of and from the set of unknowns to be originally solved for is achieved. These two quantities, as well as , are the recovered ones.

Since the data are posed on , the process initialization (i.e. iteration ) is started with the solution of global problem (), ().

In order to fit equation () with in the concept of non-intrusiveness, must be treated as a known quantity. This defines the order in which equations () and () are solved at any iteration : the solution of () precedes the solution of (). In this case, as desired, the stiffness matrix remains unaltered; the access to is still required, but only at the stage of solving (), not ().

For solving (), must be also known. At , can simply be taken from the previous loading step. At , we use coupling equation () for the extraction of , assuming is already known. This defines the order in which the global and local problems are solved: at any iteration starting from , the solution of (L) precedes the solution of ().

We also notice that:

Coupling equation () provides the boundary condition for of the local problem (L).

Coupling equation () is used for the recovery of .

The summary of the solution operations to be performed at any iteration n of the procedure, excluding the initialization step (), is as follows:The detailed scheme, including the iteration , is depicted in Table 3. Note that in all equations in the table we omit and .

Table 3

Non-intrusive iterative solution process for (), (), (L), and (), (), () at a fixed loading step l

* See “Staggered process for the local problem” section, ** see “Accuracy/convergence check” section

solution of local problem (L) coupled with (),

recovery phase using () and (),

solution of global problem (),

recovery phase using ().

Staggered process for the local problem

Solution of the local system in Table 3 at the given global/local iteration 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 , 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 ) is chosen as follows. At iteration (and staggered iteration ) the values known from the previous loading step are used as the initial guess. At (and staggered iteration ), we naturally take .

At any fixed iteration , the accuracy check for the solution triple obtained at the staggered iteration is performed as follows:where, to recall, is given by (31). If (41) is fulfilled—note again that in the following numerical test, – the staggered process is stopped, we set and perform .

Accuracy/convergence check

Derivation of the convergence and stopping criteria for the global/local iterative solution process in Table 3 is rather straightforward. Indeed, at any iteration , the solution outcome is denoted as . Plugging this in equations (), (), (L), (), (), () and comparing the obtained outcome with the corresponding equations in Table 3, it is straightforward to locate the imbalanced quantities:andTherefore, the solution accuracy at n is measured by the quantitywhere 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 is recovered (post-processed) from the right-hand side of (43). The stopping criterion for the global/local loop can then be defined aswith to be prescribed. Owing to the ‘nested in’ nature of the staggered process, it has to be . Recalling that , we set .

In our computations (see Table 3), we use a more convenient form of the stopping criterion. Settingwe define , and use this quantity to now checkThis choice of fulfills the requirement , which is stipulated by (44), (45), and the already prescribed above magnitudes of and .

Since the quantity naturally stems from the global/local solution accuracy check , 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 and are only iterative convergence indicators, but none of them provides an adequate check of the solution accuracy. In particular, since measures, though implicitly, the displacement continuity—a match between and across (recall that the traction continuity—a match between and 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 , this reads: given the triple known from the iteration , as well as ‘recovered’ at the iteration n, we solve

for and setWe 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 be a finite element partition of into triangles or quadrilaterals, I be the number of nodes in , and , be the nodal shape function associated with the node i and supported on the collection of elements in that share i. Finally, let a scalar-valued quantity represent the nodal value.

The standard discretization of the solution of the reference problem (using Voigt’s notation) is as follows:where is the interpolation/basis matrix and is the displacement nodal-vector, andwhere is the interpolation/basis matrix and d is the phase-field nodal-vector. constitutes both matrices and . The corresponding representation of the gradients reads:where is a matrix, andwhere is a matrix. The problem test functions and w, as well as their gradients are discretized accordingly.

For the global/local formulation, we assume the existence of the partitions and of and , respectively. Using this and adopting the above notations, the solution discretizations are given bysuch thatNote that the superscripts G and L are introduced in the corresponding definitions of , and , .

To construct the discretization of the Lagrange multipliers , , and the supplementary quantity on we account for the following. In the most general situation, three distinct partitions of may be assumed: the restrictions of and —denoted as and , respectively—which serve to create the corresponding bases for (also ) and , and the ‘independent’ partition to be used for creating the basis for . Introducing the three related basis matrices , and , we writeIn the following, for our numerical example, we assume that: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].

the partitions , and match (this is usually termed a ‘matching case’);

the basis in the global and local domains is identical, that is, ;

the basis on the interface is obtained from by the corresponding restriction, that is, ;

the nodal shape functions composing bases and are piecewise linear.

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.

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

Fig. 7

Comparison of the load-displacement curves

Fig. 12

Comparison of the load-displacement curves

The geometric data are as follows (all given in mm): , , , with the hole diameters and . The material data are: Young’s modulus , Poisson’s ratio and the critical energy release rate . The characteristic length in the phase-field formulation is . We consider the plane-strain situation.

The algorithmic parameters are: the loading with and the increment size , the tolerance magnitudes are , , and .

We recall that we use -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 mm. The former fulfills the heuristic requirement 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.

Fig. 6

Finite element mesh used for the discretization a of the reference domain , b of the global and local domains and , respectively

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.

Fig. 8

Comparison of the computed phase-field profiles

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 , and 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.

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 and

The first important observation is that , 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 (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.) 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.

Fig. 10

Time-displacement curves in terms of ‘accumulated time’

Fig. 11

Time-displacement curves in terms of ‘time per loading step’

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 significantly 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 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 -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 displacementswhere, owing to (), one haswith as a -matrix representation of , andAitken’s method modifies (50) at any iteration introducing the damping factor s.t. as follows:whereas the Quasi-Newton correction modifies (50) at any iteration by replacing the matrix with s.t. , thus resulting inIn (51), we explicitly havewith . In (52), the SR1 update formula implieswhereas the Broyden update reads asIn both cases . 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 ifwith a constant . Otherwise, we simply set . 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.

Fig. 15

Comparison of the time-displacement curves presented in different formats

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

Fig. 13

Convergence behavior of the different versions 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

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 non-uniqueness 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.

Fig. 14

Comparison of the phase-field profiles computed with the various acceleration/relaxation versions of the global/local approach

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 -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 -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 the possibility of the adaptive choice of the local domain . The representation result for d, namely,obtained in Ambati et al. [15], p. 392, equation (34) suggests that , or rather the history field , may serve as an “indicator” for the adaptive choice of . In (56), is an averaging volume, , and is a bell-shaped function, typically a Gaussian, such that holds. These issues are all open for further research.