Computational homogenisation of phase-field fracture

In this manuscript, the computational homogenisation of phase-field fractures is addressed. To this end, a variationally consistent two-scale phase-field fracture framework is developed, which formulates the coupled momentum balance and phase-field evolution equations at the macro-scale as well as at the Representative Volume Element (RVE) scale. The phase-field variable represent fractures at the RVE scale, however, at the macro-scale, it is treated as an auxiliary variable. The latter interpretation follows from the homogenisation of the phase-field through volume or a surface-average. For either homogenisation choices, the set of macro-scale and sub-scale equations, and the pertinent macro-homogeneity satisfying boundary conditions are established. As a special case, the concept of selective homogenisation is introduced, where the phase-field is chosen to live only in the RVE domain, thereby eliminating the macro-scale phase-field evolution equation. Numerical experiments demonstrate the local macro-scale material behaviour of the selective homogenisation based two-scale phase-field fracture model, while its non-selective counterpart yields a non-local macro-scale material behaviour.


Introduction
An in-depth understanding of fracture (initiation and propagation) processes in materials is essential for the prediction of fracture-induced failure in engineering structures. To that end, the past century has seen a thrust towards developing theoretical approaches to help gain a deeper understanding of fracture processes. The earliest theoretical approach, developed by Griffith and Taylor (1921) reasoned that fracture propagation occurs if the energy release rate reaches a critical value. Much later, in an alternate approach Irwin (1957) postulated a fracture propagation criterion based on stress-intensity factors. However, both theories were unable to predict the initiation of fracture and explain topologically complex (branching, merging, kinking and curvilinear) fractures. However, these limitations were eliminated with a variational model based on energy minimisation of the fractured continuum (Francfort and Marigo, 1998). The numerical implementation of the same was proposed in Bourdin et al. (2000), motivated by the Ambrosio-Tortorelli regularisation of the Mumford-Shah potential (Mumford and Shah, 1989). An auxiliary variable, the phase-field was introduced that interpolates between the intact and the broken material states. This lends the name phase-field fracture model (PFFM).
In the past decade, there has been an increased interest in PFFM, primarily due to its ability to predict fracture initiation and handle Some other numerical techniques adopted for the PFFM include the use of the dissipation-based arc-length method (May et al., 2015), modified Newton method (Wick, 2017b) and error-oriented Newton method (Wick, 2017a) . While most studies are focused on quasi-static analyses, Borden et al. (2012) adopted a monolithic Newton solver for dynamic (brittle) fracture simulations. Therein, it was reported that the physically limited crack tip velocity prevents full fracture within a single timestep if the timestep sizes are chosen adequately. As an alternative to monolithic solvers, a staggered (alternate minimisation) solver was suggested in Bourdin (2007) in conjunction with 'crackset' irreversibility. Later, in Miehe et al. (2010a), the 'crack-set' based irreversibility was replaced by an implicit 'history variable' based irreversibility. Although the staggered solver is numerically robust owing to the convexity of the energy functional w.r.t displacement and phase-field separately, it is computationally expensive compared to monolithic solvers (Gerasimov and De Lorenzis, 2016). Yet another aspect connected to computational efficiency is the adaptive refinement of the mesh. In particular, the phase-field fracture model requires extremely fine meshes in the phase-field transition zone. In this regard, fixed uniform meshes could be used when the fracture path is not known in advance. However, if it is known, certain sub-domains of the mesh could be pre-refined. More elegant ways in the form of errororiented mesh refinement (Burke et al., 2010;Wick, 2016), refinement based on the phase-field reaching a certain threshold (Heister et al., 2015) and local increase of the tensile energy (Klinsmann et al., 2015), and multi-level hp refinement using the finite cell method (Nagaraja et al., 2019) exists in the phase-field fracture literature. Despite these advancements, the development of robust and computationally efficient solution and meshing techniques are still topics of active research.
So far, the studies pertaining to the PFFM are limited to a single scale. In the context of multi-scale approach, the PFFM has been used in conjunction with the Multi-scale Finite Element Method (MsFEM) to simulate brittle fracture (Patil et al., 2018a), failure in composites (Patil et al., 2018b) and fractures in highly heterogeneous materials (matrix with voids and/or inclusions) (Patil et al., 2019). The MsFEM assumes a fine-scale domain embedded within a coarse macro-element. The fine-scale features (voids, cracks and other heterogeneities) are then captured using multi-scale basis functions, computed numerically onthe-fly. However, these fine-scale features if several magnitudes lower in size than the domain itself renders the fine-scale problem expensive. A cheaper alternative can be formulated on assuming separation of scales which allows a comparatively smaller fine-scale (referred to as sub-scale in this manuscript) domain in a computational homogenisation framework. The separation of scales was assumed in a study involving porous media (He et al., 2020), using the Finite Element-Heterogeneous Multi-scale Method (FE-HMM). However, only the elastic tensor was 'homogenised' owing to the presence of microstructural pores, and the phase-field evolution equation was not solved at the micro-structural level. This indicates that the microstructural fractures/cracks were not accounted for. In yet another study (Fantoni et al., 2019), asymptotic homogenisation of the microstructures were performed offline for varying phase-field values. The homogenised constitutive tensor was then obtained using a closed-form expression based on two-scale asymptotic homogenisation and interpolation of the phase-field variable. Such a method, however, requires that the offline computations include all possible failure topologies of the microstructure. This could be a challenging task in the case of topologically complex microstructural features. An elegant alternative would be to introduce a framework, wherein the coupled momentum balance and phase-field evolution equations are established the macro-scale as well as at the microstructural (RVE) scale, along with adequate computational homogenisation technique. However, to the best of the authors' knowledge, such a framework has not been developed yet.
In this manuscript, a two-scale phase-field fracture framework is developed using the Variationally Consistent homogenisation (VCH) framework (Larsson et al., 2010b) and the relevant computational homogenisation aspects are discussed. The VCH framework provides an elegant procedure to derive pertinent scales for a hierarchical multiscale problem, from its fully resolved fine-scale problem 2 . The critical ingredient of the method lies in the conjunction of the Variational MultiScale method (Hughes et al., 1998) and the separation of scales adopted through classical (first-order) homogenisation. The Hill-Mandel macro-homogeneity conditions (Hill, 1963(Hill, , 1984Nemat-Nasser, 1999) are fulfilled through equivalent Variationally Consistent Macro-homogeneity Conditions. The advantages of the VCH framework lies in its applicability in the homogenisation for a general class of problems, and in establishing scale-bridging strategies. The VCH framework has been used to derive multi-scale models in porous media (Larsson et al., 2010a;Ohman et al., 2013;Jänicke et al., 2020), gradient-enhanced visco-plastic dissipative materials (Runesson et al., 2017), and computational homogenisation of micro-fractured continua using the eXtended Finite Element Method (XFEM) (Svenning et al., 2016b(Svenning et al., , 2017, to cite a few. However, the VCH framework has not been explored yet in the context of smeared-type fracture or damage models. In the view of existing literature on the phase-field fracture model and the VCH framework, discussed in the preceding paragraphs, a twoscale phase-field fracture framework addresses the two-fold research gap, viz., (i.) the lack of a multi-scale framework wherein the coupled momentum balance and phase-field are formulated at the macro-scale and RVE scale, and (ii.) extending the VCH framework to smearedtype (phase-field) fracture model. Moreover, the two-scale phase-field fracture framework is generic in the sense that it allows different choices pertaining to computational homogenisation of the microstructural quantities. This aspect is explored at length in this manuscript, with (i.) volume and surface-average based homogenisation measures, and (ii.) selective homogenisation in the context of the phase-field variable. In particular, the novel contribution of this manuscript are: • the formulation of a variationally consistent two-scale phasefield fracture framework, that allows different models based on computational homogenisation choices; • establishing the space-variational (Euler-Lagrange) equations and pertinent homogenised dual quantities for three different twoscale phase-field fracture models, derived adopting volumeaverage, surface-average and selection homogenisation measures.
The focus of this manuscript lies in the computational homogenisation aspects of the different two-scale phase-field fracture models and not in the representativeness of real random media. Therefore, the RVEs used throughout this manuscript are artificially created and designed to demonstrate the underlying micro-structural features. However, in the case of real random media, the existence and size determination of RVEs (or Statistical Volume Elements) requires careful investigation. For more on this aspect, the reader is referred to Ostoja-Starzewski (2006) and Gitman et al. (2007).
This manuscript is organised as follows: In Section 2, the reader is introduced to the Phase-Field Fracture Model (PFFM), its underlying energy functional and the set of coupled space-variational (Euler-Lagrange) equations. The two-scale phase-field fracture framework is then developed in Section 3. Within this framework, a family of twoscale phase-field fracture models are developed, based on different homogenisation choices. Thereafter, in Section 4, a numerical investigation is carried out on the artificially created RVEs in the context of constraints (Dirichlet, Neumann and Strongly Periodic boundary conditions, and domain or surface constraints) and pertinent upscaled (homogenised) quantities for the different two-scale phase-field fracture models. A model multi-scale FE 2 problem is presented in Section 5 and results from the simulations are discussed. Finally, Section 6 lays down the concluding remarks of this manuscript. R. Bharali et al.

Notation
The following notations are strictly adhered to in this manuscript: • Zero-order tensors (scalars) are represented using italic letters, firstorder and higher order tensors are represented with bold-faced letters. • A function with its arguments , is written in the form ( , ), whereas a variable with operational dependencies , is written as [ , ]. • The volume and surface-average of a quantity, say , are denoted as ⟨ ⟩ □ and ⟨⟨ ⟩⟩ □ . They are defined later in the text, in Section 3.1.

Phase field fracture model
In this section, the reader is introduced to the Phase Field Fracture Model, starting with the Francfort-Marigo energy functional (Francfort and Marigo, 1998), its phase-field regularisation and minimisation. All formulations and derivations are within the small strain continuum framework.

The energy functional
Let ∈ R dim (dim = 2, 3) be the domain occupied by the fracturing solid as shown in Fig. 1a. Its boundary is decomposed into a Dirichlet boundary ( ) D and a Neumann boundary ( ) N , such that = ( ) D ∪ ( ) N and ( ) D ∩ ( ) N = ∅. Furthermore,  denotes the crack set (a single sharp crack in Fig. 1a) in the solid.
The energy of a fracturing elastic solid is described by the Francfort-Marigo functional in Francfort and Marigo (1998) where ( [ ]) is the elastic strain energy density function, p denotes the tractions on ( ) N , and last integral pertains to fracture energy, where is the Griffith fracture toughness. The elastic strain energy density function is defined as where and are the Lame parameters, is a second-order identity tensor, ∶ → R is the displacement, and is the symmetric strain tensor given by, In Fig. 1b, the sharp crack topology is regularised by introducing a diffusive (smeared) fracture zone of width > 0, and an additional scalar auxiliary variable . The fracture surface  is now replaced by the continuous variable ∶ → [0, 1], where 0 corresponds to the intact state and 1 indicates a fully formed crack. Accordingly, the integrand in (1) is replaced by an elliptic Ambrosio-Tortorelli function, ( , cf. Bourdin et al. (2000). The energy functional for the fracturing solid now attains the form In the event of a fracture occurring in a solid, the strain energy of the solid is expected to decrease. Additionally, in this manuscript, it is assumed that fractures occur only under tensile loading. Both these requirements are met upon introducing an additive split of the elastic strain energy density into a tensile part + and a compression part − , such that a monotonically decreasing degradation function ( )+ acts only on + (Miehe et al., 2010a). This results in the modified energy functional where ( ) = (1 − ) 2 and = 1 − 10 (a small term that prevents numerical singularity 3 ). The tensile-compressive split of the strain energy density are given by In the above relation, ± [ ] is defined as where represents the th eigenvalue of the strain, and is its corresponding normalised eigenvector.
The subsequent sections would involve the space-variational (Euler-Lagrange) equations pertaining to the energy functional in (5). In this context, the Cauchy tensile and compressive stresses are defined as (8)

The space-variational formulation
In order to predict the fracture path in a solid occupying the domain , the energy functional in (5) should be minimised w.r.t. the solution variables, vector-valued displacement and scalar-valued phase-field . This has to be further augmented with an additional requirement of fracture irreversibility (no healing of fractures is permitted) and pertinent Dirichlet and/or Neumann boundary conditions. This results in a constrained minimisation problem that reads: Problem Statement 1. Find and for all times ∈ [0, ] such that, ( ), ( ) ( ) and.
Here, [0, ] refers to the time interval of interest. In this manuscript, the time refers to a loading step, instead of the actual time (quasi-static loading). The system in (9)  does not lead to loss of generality of the original problem (4). The space-variational (or the Euler-Lagrange) equations are derived by taking the first variation of the energy functional w.r.t. its solution variables and . This results in the following: The trial and test spaces are defined as The left superscript in (12b) refer to the previous step in (pseudo) time. For brevity, the superscript ( + 1) over the variables and solution fields in the current step in time is avoided. ■ Note that the variational inequality (10b) in Problem Statement 2 stems from the fracture irreversibility requirement ≥ . The treatment of the fracture irreversibility is a widely discussed topic when it comes to developing computationally efficient and robust equality-based solution techniques. In this context, Gerasimov and De Lorenzis (2016) suggested a penalisation approach to (10b). Adopting an alternative approach, Heister et al. (2015) proposed the use of a semi-smooth Newton method developed by Hintermüller et al. (2002). Yet another alternative was suggested in Wick (2017a), where an augmented Lagrangian method was developed using the Moreau-Yoshida regularisation. Note that all of the aforementioned literature advocated the use of a monolithic solver. However, in Miehe et al. (2010a), a staggered (alternate minimisation) solution technique is proposed, where the fracture irreversibility is enforced implicitly using a 'history term' , defined as the maximum accumulated tensile energy over the loading history. Based on the assumption that the fracture is driven by the tensile energy, the authors in Miehe et al. (2010a) postulated that the replacing the tensile energy term + ∶ [ ] in (10b) with  would ensure the fracture irreversibility 4 . Mathematically,  is given by where,  is the history term computed in the previous step in (pseudo) time. Note that substitution of the history term  in place of + ∶ [ ] in (10b) changes the variational inequality formulation in the Problem Statement 2 to a variational equality formulation that reads: The trial and test spaces are defined as The above set of equations are solved using an alternate minimisation algorithm, wherein, (14a) is solved, followed by computation of  using (13) and solving (14b). This sequence is repeated iteratively until the error measure defined as is less than a certain tolerance. In the above relation, represents a degree of freedom, is the number of fields (displacement and phase-field in this manuscript), corresponds to the number of degrees of freedom of type , and the subscript + 1 indicates the current iteration. Moreover, the set of equations are augmented by time-dependent Dirichlet and/or Neumann boundary conditions, stated earlier in Problem Statement 2. Also, note that the trial and test spaces for the phase-field in this equality-based formulation differ from variational inequality-based formulation in Problem Statement 2. ■ In order to have a concise representation of the space-variational Eqs. (14a) and (14b), the quantities dual to the strain , phase-field and its gradient are defined as, respectively. This allows re-stating (14a) and (14b) in the compact form

Variationally consistent two-scale phase-field fracture framework
In this section, a two-scale phase-field fracture framework is developed. The framework is developed using the Variationally Consistent homogenisation (VCH) technique proposed in Larsson et al. (2010b). In brief, the VCH technique replaces a fine-scale problem with a macro-scale problem, such that every macro-scale material point is associated with an RVE. This is made possible upon introducing running average approximations of the integrand in the spacevariational (Euler-Lagrange) equations, and separation of scales using first-order homogenisation. These aspects are treated in detail in the following sub-sections. Later in the text, the computational homogenisation aspects pertaining to volume or surface-average homogenisation measures as well as selective homogenisation of the phase-field variable are discussed at length. These include establishing prolongation/homogenisation rules and deriving the relevant homogenised dual quantities.

Running averages
The VCH technique allows a continuous macro-scale problem in the domain , upon introducing a sub-scale RVE □ | at each macroscale material point ∈ . Any integrand on is approximated as a quantity averaged over □ . For instance, an integrand in is obtained through volume-averaging on □ as Incorporating the volume-averaging definition (19a) and (19b) in (18a) and (18b) yields Note that each term within the angular brackets ⟨⋅⟩ □ are evaluated on the RVEs, located at macro-scale material points (also referred to as Gauss/integration points in a numerical integration scheme). Furthermore, the prescribed tractions p and p are assumed to be appropriately homogenised.

Remark 1.
The VCH framework is generic in the sense that there is no restriction on the definition of the averaging that replaces an integrand. For instance, the integrand, could also be defined through surface-average approximation over the RVE boundary □ as Also, the volume-averaging could be carried out over a part of the RVE domain. An example of such an approach is averaging over a failure-zone (Nguyen et al., 2010) (not pursued in this manuscript).
In the next sub-section, the RVE solution fields , and the corresponding test functions , would be additively decomposed into a macro-scale contribution and an RVE scale fluctuation adopting the first-order homogenisation technique.

Scale transition
Scale transition enables to define the RVE solution fields and their corresponding test functions in terms of their macro-scale counterparts (denoted with an overbar in this manuscript). To this end, first, the solution fields and are additively decomposed into a macro-scale contribution (with a superscript M) and an RVE scale fluctuation (with a superscript s), Thereafter, the macro-scale contributions M and M are assumed to be linearly varying (first-order) Taylor series expansions about the smooth macro-scale solution fields and (an approach, consistent with the first-order homogenisation technique). This results in where = [ ]| , = | and = | . For the sake of brevity, | is dropped in the subsequent text of this manuscript. Note that in (23a), the skew-symmetric part of the displacement gradient is excluded due to rigid body invariance. Consequently, the definition of the symmetric strain in (3) is adopted. Furthermore, the test functions and also follow the same additive decomposition and linearly varying macroscale contributions using first-order Taylor series expansion about their corresponding macro-scale test functions and . This procedure of mapping a macro-scale field to its contribution in the RVE (sub-scale) counterpart is termed as prolongation.

Macro-scale problem
The macro-scale space-variational (Euler-Lagrange) equations for the phase-field fracture problem is obtained upon testing (20a) and the trial and test spaces are defined as Remark 2. In the above formulation, a tacit assumption is made allowing the identification of appropriately homogenised Dirichlet ( p , p ) and Neumann ( p , p ) values, analogous to those used in Problem Statement 3.
Remark 3. Note that the macro-scale phase-field evolution equation (24b) is different from the original formulation (18b), due to the presence of the additional non-local term in the former. This additional term stems from the higher-order term in the prolongation of M (consistent with the first order homogenisation technique).

RVE problem
The RVE space-variational (Euler-Lagrange) equations are obtained upon localising (20a) and (20b) to each RVE domain □ . To this end, (20a) and (20b) are tested with the fluctuating test functions = s and = s .

RVE weak/strong periodicity problem
The canonical form of the RVE problem, according to the weak micro-periodicity format (Larsson et al., 2011) is stated as − 1 with pertinent spaces Note that the RVE phase-field bounds ∈ [0, 1] are self-regulated by the weak form equations and need not be incorporated in the space (28b). Moreover, the Lagrange multipliers ( ) , ( ) and ( ) are related to the macro-scale quantities defined in Problem Statement 4 as = 1 with the jump operator • □ defined as • □ = • + − • − . The superscripts + and − are indicative of the RVE boundaries + □ and − □ respectively, as shown in Fig. 2a. The RVE boundary + □ has a positive outward normal, and − □ has a negative outward normal in a Cartesian coordinate system . ■ Incorporating the constraint Eqs. (27c), (27d) and (27e) in the RVE problem ensures the fulfilment of the Hill-Mandel macro-homogeneity conditions (Hill, 1963(Hill, , 1984Nemat-Nasser, 1999). A formal proof of the same is presented in Appendix A of this manuscript.
Problem Statement 5 allows an independent discretisation of the Lagrange multipliers ( ) and ( ) from that used for the displacement and phase-field. As elucidated in Larsson et al. (2011), theoretically, using the same discretisation for the solution fields and the Lagrange multipliers at the RVE boundary enforces a strongly periodic boundary condition whereas, adopting a single Lagrange multiplier element for the RVE edge as shown in Fig. 2b results in a Neumann boundary condition. However, Svenning et al. (2016a) showed that LBB-stability is ensured only if the solution fields ( and in this manuscript) mesh have at least one node inside each of their corresponding Lagrange multiplier ( ( ) , ( ) ) elements. In this manuscript, strongly periodic boundary conditions are enforced through restrictive enrichment of the displacement and phase-field test and trial spaces. This results in: using the test and trial spaces Remark 4. The RVE Weak/Strong Periodicity problem requires fixing one of the RVE corner nodes (bottom left node is chosen in this manuscript) in order to restrict rigid body translations.

RVE Neumann problem
The RVE Neumann problem arises from choosing trial spaces for the Lagrange multipliers ( ) ∈ T □ and ( ) ∈ Q □ such that Here, and are homogenised quantities dual to and respectively, and is the surface normal. Adopting the aforementioned trial spaces in Problem Statement 5, along with trivial manipulation results the RVE Neumann problem that reads: R. Bharali et al.
using the test and trial spaces Remark 5. The RVE Neumann problem requires fixing one of the RVE corner nodes (bottom left node is chosen in this manuscript) in order to restrict rigid body translations. Furthermore, the RVE Neumann problem allows a small discrepancy in the context of consistency with the initial values. However, this discrepancy exists only in the first step of a fully coupled two-scale analysis.

RVE Dirichlet problem
The RVE Dirichlet Problem results from choosing to enforce displacement and phase-field values on the RVE boundary. This results in: using the test and trial spaces Remark 6. Enforcing a Dirichlet boundary condition as stated in Eq. (36b) on the phase-field would lead to a 'undesirable' conflict in the presence of initial fractures on the RVE boundary. It is presented in this manuscript solely as a proof of concept and for the sake of completeness.
So far, the macro-scale kinematic quantities , and have been defined as the volume-average of their RVE counterparts. Therefore, the macro-scale problem in Section 3.3 and the RVE problems (Problem Statements 5-8) derived in this section constitute a volume-average based two-scale phase-field fracture model. In the next sub-section, a surface-average based two-scale phase-field fracture model is introduced that defines the macro-scale phase-field as the surface-average of its RVE counterpart.

Surface-average based two-scale phase-field fracture model
The surface-average based two-scale phase-field fracture model defines the macro-scale phase-field as the surface-average 5 of the RVE phase-field , keeping the other kinematic quantities [ ] and same as in the volume-average based two-scale phase-field model. This results in the constraints (27e), (30c), (33e) and (35c) being replaced by with evaluated at a material point ∈ . Consequently, the term ( ) ⟨ ⟩ □ is replaced by ( ) ⟨⟨ ⟩⟩ □ in the RVE phase-field evolution equation. The macro-scale equations in Problem Statement 4, however, remain unchanged. The macro-scale and RVE problem statements for the surface-average based two-scale phase-field fracture model are not explicitly stated in this manuscript for brevity.

Remark 7.
Note that the macro-scale phase-field being defined as the volume or surface average of its RVE counterpart, is not indicative of fracture on the macro-scale. Rather, it must be treated as an auxiliary macro-scale variable.

Selective homogenisation based two-scale phase-field fracture model
Yet another variant of the two-scale phase-field fracture model is proposed in this sub-section, based on 'selective homogenisation' of the solution fields. 'Selective homogenisation' refers to the selective upscaling of the solution variables from the sub-scale to the macro-scale. In this regard, a simple choice would be to discard any notion of the phasefield variable at the macro-scale scale, i.e., the phase-field is assumed to live only on the RVE domain. This 'special case' is not new in the computational homogenisation literature. For instance, the pressure field was assumed to live only on the RVE domain in liquid-phase sintering (Ohman et al., 2013), Stokes' flow  and fluid transport in fractured media (Pollmann et al., 2020) problems, to cite a few.
For the phase-field fracture problem, assuming the phase-field to live only on the RVE domain leads to the non-existence of the macroscale phase-field evolution equation (24b), thereby circumventing the need to extract the homogenised quantities dual to the macro-scale phase-field and its gradient. The absence of the macro-scale phase-field evolution equation is expected to reduce computational cost compared to the volume-average and surface-average based two-scale phase-field fracture models. However, assuming the phase-field only as an RVE quantity would result in an RVE-based local material model at the macro-scale, similar to the local damage model in continuum damage mechanics. The RVE-based local dissipative material model would render the macro-scale problem mesh sensitive (refer to de Borst et al. (1993) for more on this aspect).
As far as the RVE problems (Problem Statements 5-8) are concerned, considering the phase-field only as an RVE quantity would eliminate the need for constraints (27e), (30c), (33e) and (35c). However, constraints on the RVE phase-field must be enforced such that the Hill-Mandel macro-homogeneity conditions are satisfied. This is achieved through Neumann and Periodic boundary conditions with a zero macro-scale phase-field gradient.

Single-scale RVE numerical study
The single-scale numerical study extracts the homogenised dual quantities for the different two-scale phase-field fracture models, discussed earlier in Section 3. To this end, a set of numerical experiments are carried out on artificially created RVEs. The RVEs differ in material constituents and/or initial fracture topology. The initial fractures are modelled by defining interfaces within the RVE domain and prescribing = 1 on these surfaces. All material and geometric parameters pertaining to the RVEs are addressed in the next sub-section. The subsequent sub-sections conduct a three-fold numerical investigation, where • Study I computes the homogenised dual quantities pertaining to the selective homogenisation based two-scale phase-field fracture model (refer to Section 3.6), • Study II compares the volume-average and surface-average based two-scale phase-field fracture models, based on their homogenised stress-strain response, and • Study III involves a parametric study in order to ascertain the influence of the macro-scale phase-field gradient on the homogenised dual quantities in the volume-average and surfaceaverage based two-scale phase-field fracture models.

Artificially created RVEs
Three different artificially created RVEs are considered in this manuscript with varying initial fracture topology and/or material constituents, as shown in Fig. 3. All of them are two-dimensional unit squares (in mm). Fig. 3a shows an RVE with an initial vertical fracture. This RVE is symmetric w.r.t the fracture topology. The second RVE in Fig. 3b is devoid of initial fractures, instead, the matrix is embedded with randomly placed inclusions of varying size (shown in dark blue colour). These inclusions fulfil wall-periodicity as they are allowed to penetrate through the RVE boundary and appear on the opposite edge. As such, material periodicity is invoked. Finally, Fig. 3c shows an RVE with random initial fractures, that fulfils wall-periodicity. Note that the latter two RVEs are not symmetric as far as the material and fracture topology are concerned.
The material and geometric properties for the different RVEs are presented in Table 1. Note that the matrix material remains the same in all the RVEs, and the inclusion properties apply only to the RVE in Fig. 3b. Throughout the entire numerical investigation, the RVEs are subjected to a strain-loading in the -direction. The loading is quasi-static, and the solution-based error measure (16) is adopted to terminate the iterations with a tolerance 1 − 3.

Study I
Study I pertains to the selective homogenisation based two-scale phase-field fracture model that considers the phase-field only as an RVE (sub-scale) solution field. (refer to Section 3.6 for details). As such, the problem is driven only through a quasi-static strain-loading in the -direction. Table 2 presents the strain increments adopted for the different RVEs. When the strain-loading is enforced through DBC and NBC, the phase-field evolution equation is augmented with NBC. However, when the strain loading is applied through SPBC, the SPBC is also enforced on the phase-field.
In the selective homogenisation based two-scale phase-field fracture model, the macro-scale phase-field evolution equation ceases to exist. As such, the homogenised stress (dual to the homogenised strain ) is the only macro-scale quantity that requires upscaling. Fig. 4 presents the homogenised stress-strain curves for the three RVEs with different displacement boundary conditions (DBC, NBC, and SPBC). Each subfigure corresponds to a single RVE, while the curves of different colour represent the different boundary conditions. It is observed from all the sub-figures that the phase-field variable implicitly contained in the definition of the homogenised stress (see Eq. (25a)) manifests in the form of a dissipative-type behaviour. Furthermore, the DBC is found to yield a stiffer stress-strain response in comparison to the NBC, while the SPBC stress-strain curve lies in between the DBC and NBC response. The stiff behaviour of the DBC, owing to the rather restrictive enforcement of linearly varying displacements, is established in the computational homogenisation literature. Fig. 5 shows the phase-field fracture topology at failure for the RVE with a single initial fracture. Irrespective of the applied boundary conditions, this fracture topology remains the same, i.e., through elongation of the initial vertical crack. This explains the closeness homogenised stress-strain curves in Fig. 4 with different boundary conditions. However, in the case of DBC, the fracture is not allowed to reach the RVE boundary, rather it spreads horizontally as seen in the red curve in Fig. 5a. This prevents the total loss of material integrity and results in an artificial stiffening (evident from the horizontal plateau). The artificial stiffening is, however, not observed with the NBC and SPBC as observed from the green and blue curves in Fig. 4.
Figs. 6 and 7 show the phase-field at the fracture initiation stage and at the final step of the analysis respectively, for the RVE with stiff inclusions. It is observed that NBC and SPBC results in fracture initiation on the RVE boundary (see Figs.6b and 6c) which propagate into R. Bharali et al. the RVE with increase in loading until total loss of integrity of the RVE (see Figs.7b and 7c). However, for the DBC, fracture initiation occurs inside the RVE domain and not on the RVE boundary as observed from Fig. 6a. Furthermore, similar to the RVE with single initial fracture, total loss of integrity is not achieved as the fracture is not allowed to develop at the RVE boundary. This manifests in the form of an artificial stiffening in the stress-strain curve shown in Fig. 4b. Moreover, the restrictive nature of linearly varying displacements enforced by the DBC in conjunction with stiff inclusions on the RVE boundary yields a stiffer response compared to NBC and SPBC in the pre-peak regime of the stress-strain curve.
The phase-field fracture topologies at failure for the different RVEs with varying boundary conditions poses a question as to which of them are reasonable. In this context, the DBC that results in an unphysical artificially stiffened response is ruled out. Next, the NBC circumvents the issue with the artificially stiffened response resulting in a realistic fracture pattern for the RVEs with no initial fractures on the boundaries as observed from Figs. 5b and 7b. However, when the RVE has initial fractures at the boundary, the NBC leads to widening of these existing fractures as seen in Fig. 8b, resulting in an unrealistic response. The SPBC, however, circumvents both the artificial stiffening and widening of existing boundary fractures, at the cost of wall-periodicity (see Figs. 5c, 7c and 8c). Therefore, subsequent studies in this manuscript (i.e., Study II and III ) involve only the SPBC.
Next, in Fig. 9, the macro-scale phase-field (obtained as a postprocessing step) is plotted against the homogenised strain xx (indirection) for the SPBC. The blue and the red curve correspond to the volume and the surface-averaged definition of respectively. In either case, is far below one, even after the total loss of material integrity. Thus, the macro-scale phase-field is not an indicator of a fully developed fracture. Rather, upon reaching the total loss of material integrity, the curve flattens to form a horizontal plateau. The formation of the  plateau signifies a halt in the formation of new fracture or propagation of existing fractures. Later, in Study II and III, the blue curve is used to enforce the constraint (30c) for the volume-average based two-scale phase-field fracture model. Likewise, the constraint (38), pertaining to the surface-average based two-scale phase-field fracture model is enforced using the red curve.

Study II
This sub-section concerns the volume-average and surface-average based two-scale phase-field fracture models. The numerical aspects of both models (space-variational equations and constraints) are discussed in Section 3. Similar to Study I, a strain-loading is applied to the three RVEs (cf . Tables 2 and 3), albeit using the only SPBC. Additionally, the constraint (30c) in the volume-average based two-scale phasefield model is enforced using parametrised by the blue curves in Fig. 9. Likewise, for the surface constraint (38) in the surface-average based two-scale model, is parametrised by the red curves in Fig. 9. Apart from the aforementioned constraints, the macro-scale phase-field gradient is set to zero. Fig. 10 shows the homogenised stress-strain response obtained for the different RVEs. In all the sub-figures, the blue curves correspond to the volume-average based two-scale phase-field fracture model, while the red curves belong to the surface-average based two-scale phasefield fracture model. It is observed that blue and the red curves are comparable (maximum relative difference in stresses ≈ 6% in Fig.  10b) when the surface and volume-average phase-field is imposed in a consistent manner using the curves in Fig. 9. Moreover, the fracture at the final time-step also remains similar to those presented in Figs. 5c, 7c and 8c.

Study III
This sub-section extends Study II in order to assess the influence of zero/non-zero macro-scale phase-field gradient on the RVE homogenised dual quantities. To this end, numerical experiments are carried out on the RVE with inclusions (see Fig. 3b). The RVE loading conditions remain the same as presented in Table 3, the only change being, the SPBC on the RVE phase-field is enforced with a , which is not explicitly set to zero. Instead, is parametrised as where x [mm −1 ] and y [mm −1 ] are constants, and ( xx ) is chosen as a linear function of the homogenised strain. 6 Based on the choice of these quantities, different parametrisations of is achieved. For instance, choosing x = y = 0 or ( xx ) = 0 results in = 0. In this study, y is set to zero and x is chosen randomly, such that the macroscale phase-field does not result in ∉ [0, 1]. Appendix B explains this aspect in detail. Fig. 11 presents the homogenised dual quantities for the volumeaverage based two-scale phase-field fracture model. The homogenised dual quantities are defined in the Problem Statement 4 (see Eqs. (25a)-(25d)). The homogenised stress shown in Fig. 11a is dual to the homogenised strain. It is observed that the stress-strain response is objective w.r.t. the chosen x values. As the macro-scale phase-field gradient is parametrised using x , the aforementioned observation indicates that the homogenised stress-strain is not influenced by the macro-scale phase-field gradient. The dual quantity , defined in (25d) represents the volume-average of the imbalance between the fracture driving and resisting forces, excluding the gradient term. It is dual to the macro-scale phase-field . Fig. 11b shows that is objective w.r.t the chosen macro-scale phase-field gradient parametrisation. Finally, the homogenised quantity x + x 7 dual to the macro-scale phasefield gradient x is presented in Fig. 11c. This dual quantity does exhibit a dependence on the chosen macro-scale phase-field gradient parametrisation. This behaviour is attributed to varying local phasefield gradients within the RVE in the vicinity of the fracture zone, with different values of x . Moreover, on comparing Fig. 11c with Fig. 11d, it is observed that x is the dominant term in the overall homogenised quantity x + x , dual to x . Fig. 12 presents the homogenised dual quantities for the surfaceaverage based two-scale phase-field model. The homogenised dual quantities are defined in (25a)-(25d). Fig. 12a presents the homogenised stress-strain response, which is found to be objective w.r.t the chosen parametrisation of the macro-scale phase-field gradient. This observation is similar to one with the volume-average based twoscale phase-field fracture model in Fig. 11a. However, the homogenised quantity , dual to the surface-averaged macro-scale phase-field does exhibit a dependency on the chosen macro-scale phase-field gradient parametrisation. It is important to note that evolution in the surfaceaverage based model differs from the volume-average based model (cf. Figs. 12b and 11b) since they are dual to different quantities, volumeaveraged phase-field and surface-averaged phase-field. Furthermore, the homogenised quantity dual to the x-component of the macro-scale phase-field gradient x also exhibits a dependency on the chosen x values, as seen from Fig. 12c. This behaviour is similar to that observed in the case of the volume-average based two-scale phase-field model (cf. Figs. 11c and 12c). The reason for this behaviour is mentioned in R. Bharali et al.   the previous paragraph. Finally, for the surface-average based two-scale phase-field model too, x remains the dominant term in the overall quantity x + x (cf. Figs. 12c and 12d). Both, volume and surface-average based two-scale phase-field fracture model yield dual quantities to the macro-scale phase-field and its gradient, in addition to stress. In this context, it is imperative to carry out a fully coupled two-scale simulation in order to ascertain the effect of these model choices on the macro-scale structural behaviour (for instance, the load-displacement relation). The next section deals precisely with this aspect.

Multi-scale FE numerical study
In this section, the two-scale phase-field fracture models based on selective homogenisation and volume-average homogenisation of the phase-field (presented in Section 3) are investigated in the context of a fully coupled two-scale application. To this end, a one-dimensional uniaxial strain macro-scale problem is set up as shown in Fig. 13a. The one-dimensional bar is discretised with four linear elements, 1 metre each in length. The bar is fixed at the left end and loading is applied at the right end in the form of prescribed displacement. Moreover, the cross-sectional area is set to unity apart from the element adjacent to the fixed boundary, where the area has been reduced by 10%. This has been done in order to induce a localisation in that element. Finally, note that all lateral strains are set to zero, in order to ensure a one-dimensional continuum behaviour.
As shown schematically in the two-scale problem in Fig. 13a, each macro-scale Gauss point is associated with a two-dimensional RVE. In this regard, the RVEs with stiff inclusions and random fractures are not chosen for this study as they would require pre-refinement of the mesh along rather complex fracture path to reduce computational expense. Instead, the RVE with a single vertical fracture is chosen for this study 13a and 3a). The material properties remain the same as in Table 1.
The RVE mesh is pre-refined in the expected fracture propagation subdomain as shown in Fig. 13b in order to reduce the computational expenses. The element-size in the sub-domain containing the fracture is set to half of the length-scale parameter in accordance with the recommendations put forward in Miehe et al. (2010b). Moreover, for the volume-average based two-scale phase-field fracture model, a stationary analysis is carried out solely using the phase-field evolution equation to ascertain the initial value of the macro-scale phase-field. Finally, the solution-based error measure (16) is adopted to terminate the iterations with a tolerance 1 − 3.
Fig. 14 presents the macro-scale load-displacement curves for the two-scale phase-field fracture models with selective and volumeaverage based homogenisation of the phase-field. Furthermore, for the markers in these curves, the corresponding macro-scale phase-field at the macro-scale Gauss Points (GPs) are presented in Fig. 15. Note that for the selective homogenisation based two-scale phase-field fracture model, is computed as a post-processing quantity as the macroscale phase-field evolution equation does not exist. Such a modelling choice also renders a local macro-scale behaviour as observed from Fig. 15a, where grows only in one element beyond the peak load (indicated with a cyan and blue markers in Fig. 15a). However, in the case of the volume-average based two-scale phase-field fracture model, is distributed across all the elements from the peak load until failure (indicated with a green, orange and purple markers in Fig.  15b), thereby exhibiting a non-local material behaviour at the macroscale. This non-local macro-scale material behaviour manifests in the form of higher peak load and prescribed displacement at failure for the volume-average based two-scale phase-field fracture model compared to selective homogenisation based model, as far as the macro-scale load-displacement curves are concerned.
The numerical investigation in this section formally establishes the proof of concept pertaining to solvability of the selective homogenisation and volume-average based two-scale phase-field fracture models in the context of solvability of the fully coupled macro-scale and RVE problems. Furthermore, for the volume-average based two-scale phasefield fracture model, the macro-scale phase-field and its gradients is self-regulated and there is no need for artificial bounds on the macroscale phase-field gradient while solving the RVE problems. In this regard, it is important to note that is an auxiliary variable regularising the macro-scale problem, and is not indicative of failure. For instance, the RVE attached to the element close to the fixed boundary incurs a total loss of integrity when ≈ 0.1415% (shown in Fig. 15b).
Remark 8. The macro-scale length-scale for the non-selective volumeaverage based two-scale phase-field fracture model is a priori unknown. Numerical methods for estimation of this length-scale and choosing appropriate discretisation thereafter would be a part of future work.

Concluding remarks
A novel two-scale phase-field fracture framework is proposed for computational homogenisation of fractures in complex microstructures (RVEs). The framework has been developed using the Variationally Consistent Homogenisation technique (Larsson et al., 2010b), and it allows the use of several homogenisation measures (volume-averaging, surface-averaging, or selective homogenisation). Within this framework, a family of two-scale phase-field fracture models are developed using the different homogenisation measures w.r.t the phase-field variable. In this context, the macro-scale phase-field is defined as the volume-average and surface-average of its RVE counterpart, resulting in a 'volume-average based two-scale phase-field fracture model' and 'surface-average based two-scale phase-field fracture model' respectively. In both models, the phase-field represent fractures in the RVE (sub-scale), while at the macro-scale, it is treated as an auxiliary variable. The macro-scale phase-field is not indicative of material point failure (does not reach a value ≈ 1 on the total loss of integrity), however, its evolution reaches a horizontal plateau, indicating a halt in the initiation of new fracture(s) or propagation of existing fracture(s). For both, volume and surface-average based two-scale phase-field fracture models, the pertinent coupled momentum balance and phase-field evolution equation are formulated at the macro-scale and sub-scale, along with macro-homogeneity conforming prolongation/homogenisation rules. Furthermore, numerical studies on artificially created RVEs indicate that the homogenised stress-strain response is similar for both models, even though the homogenised dual quantities in the macro-scale phase-field evolution equation differ. In this regard, it is observed that prolongation of the phase-field through first order homogenisation results in a higher order term, which has a dominant contribution compared to the conventional boundary flux term. Furthermore, for a single-scale parametric RVE study, the macroscale phase-field gradient is required to be bounded in order to obtain physically acceptable meaningful results, i.e., ∈ [0, 1] everywhere within the RVE. This manuscript provides an initial estimate of the upper and the lower bound of the macro-scale phase-field gradient. The authors would like to stress that the bounds remain relevant only for parametric studies on RVEs and not in an FE 2 (Feyel, 1999) analysis.
Yet another two-scale phase-field fracture model is developed based on selective homogenisation of the phase-field variable. By construction, this model yields a local material behaviour at the macro-scale, similar to local damage model in continuum mechanics. This phenomenon has been demonstrated in this manuscript using a fully coupled two-scale application. On the contrary, in the same application, the volume-average based two-scale phase-field fracture model yielded a non-local macro-scale material. This behaviour is attributed to the presence of a macro-scale phase-field evolution equation which regularises the macro-scale phase-field. Nonetheless, the fully coupled two-scale application provides a numerical proof of concept that the macro-scale and RVE equations are solvable without the need for bounds on the macro-scale phase-field gradient.
Future studies may involve the determination of the macro-scale length-scale for the volume and the surface-average based two-scale phase-field fracture model. Also, another homogenisation measure could be incorporated in the current framework which results in a macro-scale phase-field that is indicative of material point failure (reaches a value ≈ 1 on the total loss of integrity). In this regard, the failure-zone averaging scheme proposed in Nguyen et al. (2010) offers a good starting point. Another extension could be the incorporation of weak micro-periodicity constraints (or weakly periodic boundary conditions) proposed in Larsson et al. (2011), in order to circumvent the enforcement of periodic fractures. The RVE problems in this manuscript were of saddle point nature, owing to the use of Lagrange multipliers, and were solved using a direct solver. The use of iterative solvers with an exploration into the preconditioning techniques offer yet another research dimension. Finally, the two-scale phase-field fracture framework may be extended to complex multiphysics problems (e.g., fluid flow, cement hydration) and validation studies may be carried out.

Software implementation and data
The RVE studies in the Section 4 were carried out in the software package COMSOL Multiphysics 5.5. The multi-scale FE 2 studies in Section 5 were carried out in the open-source software package openFE2 (https://github.com/rbharali/openFE2). Additional data will be made available upon request.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  The macroscopic power for the two-scale phase-field fracture given by Utilising the constraint equations in rate form (A.5a)-(A.5c) together with (A.6a) and (A.6b) results in This concludes the Hill-Mandel macro-homogeneity proof.

Appendix B. Bounds on the macro-scale phase-field gradient
In this section, the bounds for the macro-scale phase-field gradient is established that ensures that the RVE response remains realistic. To this end, Fig. B.1 presents the homogenised dual quantities pertaining to the volume-averaged two-scale phase-field fracture model with arbitrarily chosen x , while y is set to zero. The homogenised dual quantities are defined in (25a)-(25d).
It is observed that for x = 1 + 2 (green curve), the post-peak branch develops sooner compared to x = 0 and 1 + 0. This behaviour is attributed to unrealistic phase-field values at the RVE boundaries, evident from Fig. B.2a. The phase-field ∉ [0, 1], and as such the simulation results are bogus. This observation asserts the fact that in a 'single-scale' RVE analysis, the macro-scale phase-field gradient cannot be chosen arbitrarily for a parametric study. Rather, the macro-scale phase-field gradient must be chosen such that, ∈ [0, 1], everywhere in the RVE domain.  The parametrisation of the macro-scale phase-field gradient may be carried out adopting a trial and error method to arrive at a set of admissible values of x and y . However, such a procedure could be tedious in the absence of a good initial guess. This problem is circumvented using the DBC for the (36b), and requiring ∈ [0, 1] everywhere on the RVE boundary. This results in