A numerical-homogenization based phase-field fracture modeling of linear elastic heterogeneous porous media

Computational Materials Science 176 (2020) 109519

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A numerical-homogenization based phase-field fracture modeling of linear elastic heterogeneous porous media Bang Hea, Louis Schulera,b, Pania Newella, a b

T

⁎

Department of Mechanical Engineering, University of Utah, Salt Lake City, UT 84112, United States Department of Mechanical Engineering, Ecole Normale Supérieure Paris-Saclay, 94230 Cachan, France

ARTICLE INFO

ABSTRACT

Keywords: Fracture modeling Porous media Computational homogenization Finite element method Phase-field

Most porous media, such as geomaterials and biomaterials are highly heterogeneous in nature, and they contain large variations of microscopic pore structures, such as pore sizes, pore distribution, and pore shapes. The oscillation of microscopic structures is a substantial challenge in theoretical characterization and is usually ignored in continuous modeling. However, mechanical behavior of porous media such as deformation and failure, are essentially impacted by the microscopic heterogeneity which needs to be considered in modeling a porous media. This research proposes a numerical modeling framework with a capability to investigate the effect of microscopic heterogeneity on the macroscopic fracture behavior in porous media by using a numerical homogenization technique, combined with the phase-field fracture modeling method. This numerical modeling strategy computes a homogenized elasticity tensor based on microscopic heterogeneous pore structures heterogenous porous domain by solving boundary value problems at microscopic domain. The strain energy and subsequent propagation of macroscopic fractures will be updated using homogenized stiffness information. Using this numerical scheme, the microscopic pore structure’s impact on the fracture behavior through the homogenized elastic tensor will be taken into account. This multiscale technique is benchmarked against classical problems. The results highlight the importance of the underlying pore structure and reveal that both fracture strength and propagation path can be influenced by the microscopic heterogeneity.

1. Introduction Many natural and man-made substances are considered as porous media, simply because of the presence of voids in their structure. Porous media are usually characterized by the porosity that is defined as the fraction of volume of voids over total volume [49]. For some porous media, mechanical properties can be characterized based on the porosity, pore structure and skeleton properties, but such characterization is usually challenging [20,23,78,69]. This is mainly due to the fact that many porous media are highly heterogeneous and macroscopic properties1 oscillate with variation of microscopic pore structure2. For example, geomaterials and biomaterials contain large variation of micro-pores (e.g,. in term of pore diameter, 0.01–2 mm in sandstone, 0.005–0.1 mm in shales, 10–50 µ m in bone.) Some studies have illustrated the effect of micro pore structure (size, pore distribution,

orientation, etc) on the overall stiffness of porous media [89,46,83,47,90,50,48]. Fracture is the central focus of many material failure studies [21]. The study on fracture in porous media was rarely presented until most recently [21]. One possible reason is that the fracture in porous media is a multiscale and multiphysic phenomenon with uncertain micro heterogeneity which leads to difficulty in finding solution [21]. Simoni et al., have elaborately discussed multifield fracturing (thermo-mechanical fracture, fluids pressure-induced fracture, chemical reactionsinduced fracture etc.) in porous media [76]. The investigations on the influence of saturated fluid on the fracture porous media can also be found in [59,74,72]. In Cao et al.’s work [15], the stepwise fracture propagation concomitant with flow pressure fluctuation was explained. Reversely, the fracture geometries influence on the saturated flow in porous media was also explored in some studies [42,17]. Furthermore,

Corresponding author. E-mail address: [email protected] (P. Newell). For convenience, from now, we refer to macroscopic variables or quantities as macro. For instance, macro pore, macro element, macro scale, macro heterogeneity macro stiffness matrix 2 For convenience, from now, we refer to microscopic variables or quantities as micro. For instance, micro pore, micro element, micro scale, micro heterogeneity micro stiffness matrix ⁎

1

https://doi.org/10.1016/j.commatsci.2020.109519 Received 13 October 2019; Received in revised form 4 January 2020; Accepted 6 January 2020 0927-0256/ © 2020 Elsevier B.V. All rights reserved.

Computational Materials Science 176 (2020) 109519

B. He, et al.

the heterogeneity generated from mixed component and interface friction coefficients were also explored [87,40]. It shows that micro heterogeneity is critical factor affecting fracture behavior in porous media. However, in the past decade, there has been a large amount of interest arising in this subject, mainly driven by the thrives of numerical methods for fracture modeling such as extended finite element method (XFEM) [5,61,10,63,32], isogeometric analysis [13,19], damage gradient method [22,71], Discrete Element Method [66,65], Molecular Dynamics [6,82,86,14] and phase-field fracture modeling [56,58,55,33,88,85]. Many of these numerical methods have been adopted to study the fracture in porous media. For example, de Borst [21] and Simoni et al. [76] have discussed the state-of-the-art numerical methods (XFEM, phase-field, damage gradient, and discrete crack approach) in their books on how to address fracture in porous media. Many other recent numerical works aimed to explore this subject can also be found in [68,75,70,57,27,62]. Among these candidates, the phase-field method shows appealing features and has been widely used in the investigation of fracture in porous media. Phase-field fracture modeling is developed based on the pioneering work completed by Francfort and Marigo [12], Bourdin et al. [29], where they introduce a scalar variable to diffuse the sharp crack zone. This variable regularizes the zero width crack zone to a diffusive zone with finite width. Subsequently following the Griffith’s law [31], the linear elastic fracture (brittle fracture) can be modeled as a coupled mechanical and phase-field problem by minimizing the total energy [56]. Sutula et al.’s recent work [79,81] also shows that the linear elastic fracture is a minimum energy problem. The appealing feature of this so-called variational method is that it overcomes the discontinuity problem in fracture modeling. With phasefield method, the fracture path can continuously propagate in the entire domain without any special treatment to local discretization. Furthermore, phase-field method can be implemented into the conventional finite element codes such as Abaqus [64], COMOSOL [93], and Fenics [35]. Also, phase-field framework can be successfully integrated with the adaptive refinement techniques to reduce the required elements and improve the computational efficiency [34,36]. In phase-field framework, one common assumption is that the fracture starts to propagate as the surface energy is over the threshold value, usually energy release rate Gc [91,58,55,16]. Therefore, the fracture theory basis of phase-field fracture modeling of porous media is comparable to Griffith’s or Irwin’s theory for brittle or linear elastic fracture [31,12]. The literature review of phase-field fracture modeling in porous media [56,58,55,33,88,84,91,16,93] shows that the homogenized macro material property, particularly the elasticity tensor is usually used. In the linear elastic fracture mechanics (brittle fracture), the elasticity tensor plays a critical role in the evaluation of strain energy, which might impact the fracture behaviors in porous media as phasefiled fracture modeling is based on the energy minimization [73,56]. Therefore, we hypothesize that the micro pore structures impact the fracture behaviors in porous media, because the pore structures can change the stiffness of porous media. From the discussion above, macro homogenized elasticity tensor is usually difficult to obtain in the heterogeneous porous media because of the oscillatory micro heterogeneity. Many homogenization technologies have been developed to calibrate material constants, through either porosity or pore network [39,45]. Some recent research work used probability distribution of Young’s modulus to account for micro heterogeneity in their fracture modeling framework for porous media. [93]. Unfortunately, to the best of our knowledge, no published research work incorporated numerical-homogenized elasticity tensor in the fracture modeling. Numerical homogenization is an efficient tool to gather information from micro scale in heterogeneous materials. In past several decades, a set of numerical homogenization strategies have been

proposed [60,25]. Even though, it is still unclear that which one is the most suitable for given problems, the Finite Element Heterogeneous Multiscale Method (FE-HMM) shows appealing features for the homogenization of porous media. In FE-HMM, the discrete analysis is usually used to solve micro partial differential equations (PDEs), and macro heterogeneous properties are obtained using micro solution based on the homogenization theories [2,7,18,41]. Within these frameworks, macro behavior of heterogeneous material can be described without knowing their macro homogenized material property [26]. In this paper, we propose a numerical homogenization-based fracture modeling strategy based on phase-field fracture technique and FEHMM, to clarify the effect of micro structure on fracture behavior in linear elastic heterogeneous porous media. In this scheme, FE-HMM is used in the evaluation of elasticity tensor in porous media to account for the variation of micro heterogeneity impact. A robust approach is provided to extend this homogenization method to general quadrilateral elements. We also present a set of classical fracture benchmark problems to illustrate the capability of proposed approach to capture the micro heterogeneity effect. The rest of the paper is organized as follows. The phase-field fracture modeling method for linear solids is briefly introduced in the second section. The FE-HMM and associated homogenization of elasticity tensor is discussed in the third section. The concept of homogenization-based fracture modeling and verification are presented in the fourth section. Subsequently, a set of numerical examples are discussed in the fifth section, followed by the conclusion section. 2. Phase-field fracture modeling for linear solids Phase-field fracture modeling is a powerful method for numerical modeling of linear elastic fracture, which uses the diffusive variable to regularize sharp crack to a diffusive zone. In this section, the basic phase-field crack diffusion concept and governing equation for phasefield fracture modeling are introduced. 2.1. Phase-field diffusion of crack shape If we consider a one dimensional domain with a crack on surface at x = 0 , as shown in Fig. 1, the damage state of the domain can be represented by a diffusive variable , equal to zero if the material is unbroken and equal to one if it is totally broken (physically, the damage phenomenon is not a discrete phenomenon and crack width is approximately zero) [22]. Diffusive variable, can be expressed as an exponential function with respect to the position variable x , regularized using a length scale parameter l 0 :

= exp(

(1)

x /l 0)

In Eq. (1), is equal to one at the crack location and decays toward zero as the distance to the crack location increases. However, is never null. The length parameter controls the size of the diffusive area. The larger value of l 0 indicates the larger damaged area around the crack. If l 0 approaches zero, then the diffusive damage will approach its discontinuous definition. In this case, the dissipated energy in phase-field model converges to the fracture energy of sharp fracture so the convergence was reserved. As shown by Miehe et al. [56], the exponential function in Eq. (1) is the solution of the following differential equation:

(x )

l02

(2)

(x ) = 0.

Phase-field variable

takes the following extreme values:

Fig. 1. Infinite bar cut by a crack of cross section 2

at x = 0 .

Computational Materials Science 176 (2020) 109519

B. He, et al.

{10

:

as x = 0 as x ±

the crack surface l0 , and the bulk energy is the integral of the strain energy 0 over the unbroken domain. The total potential energy is equal to:

(3)

Eq. (2) is the Euler-Lagrange equation of the variational principle,

= Arg{ inf L ( )}, where W = { functional:

L( ) =

1 2

2

2

+ l02

)d .

=e

x /l 0

=

1 2

+

(5)

2 x /l0

2e

dx = l 0 .

1 1 L( )= l0 2l 0

+ l 02

2

2

d .

)=

where solid:

l0

l 0 is

0(

d ,

2

l + 0 2

(7)

(8)

.

1 2

T C0

g (1) = 0 . (12) 0(

) is calculated as:

,

(13)

=

ij

ij

Wext =

+

.

(14)

b ud

+

h udS,

(15)

where b and h are the volume and surface forces applied to the solid. As shown in Fig. 2, the surface of the body contains a part h , where the surface force is prescribed by Neumann-type boundary condition, and a part u where the displacement is imposed by Dirichlet-type boundary condition u = u 0 . From the equilibrium of the internal and external works Wint Wext = 0 , and applying the divergence theorem, we obtain the weak form of the coupled problem, for all admissible and u of the phase-field and displacement field that satisfy the homogeneous form of the Dirichlet-type boundary conditions [51]:

In the presence of a crack, the total potential energy of a linear elastic solid body is equal to the sum of the bulk energy u and the dissipated energy :

+

) Gc d ,

On the other hand, the variation of the external work is:

(9)

2.2. Governing equations of phase-field fracture model

u

)=

Wint =

From Eq. (9), the crack surface density function depends on the phase-field damage and its gradient.

=

,

(11)

2

.

l0 (

where C0 is the elastic tensor. In this paper, we did not consider splitting energy release into tension and compression portions, because most of the numerical examples presented in this paper are primarily under tension stress. Also, the decomposition of strain energy leads to the unsymmetrical stiffness matrix [54,64] in Finite Element (FE) implementation. The computation would be highly expensive. Therefore, in this study, we use full strain tensor to evaluate the strain energy. However, the energy release splitting can be added to handle compression cases without significantly impacting our framework [54,64]. From the potential energy in Eq. (11), one can obtain the variation of the internal energy:

the crack surface density function per unit volume of the

1 = 2l 0

,

l0

,

+

In Eq. (11), the strain energy

(6)

The above concept of diffusing the crack can be easily extended to the higher dimensions, such as a 2D problem shown in Fig. 2, by introducing the 2D diffusive variable, = (x, t ) , where x denotes the spacial coordinates and t is the time. From Eq. (7), the 2D regularized crack surface functional is: l0 (

)d

g ( ) = (1 )2 + k g (1) = 0, g (1) = 0, g (0) < 0, g is monotonically decreasing

Therefore, the regularized crack surface associated with the lengthscale parameter l 0 can be defined as: l0 ( ) =

0(

where the integral of l 0 represents the crack surface, g ( ) is the degradation function, representing the reduction of stiffness due to the crack in material [9] and is defined as:

This functional can be obtained by integrating a Galerkin weak formulation of the differential Eq. (2). With d = dx, we obtain

L

g( ) u

( ± ) = 0} and L is defined as the following

(0) = 1,

(

=

(4)

W

(10)

The critical energy release rate Gc is the energy dissipated during fracture per unit of new fracture surface area. For a crack to propagate, the strain energy provided to the crack must be greater or equal to this critical value. The crack surface energy dissipated by the formation of the crack is equal to the integral of the critical energy release rate over

div

+b

2(1 +

) [(1

h

0(

u )

Gc

1 l0

)2 + k ]

div (l0 ·n

h

)

udS +

d Gc l 0

·n

dS = 0. (16)

From this weak form, one reaches the strong form of the phase-field problem:

div ( ) + b = 0 Gc

(

1 l0

l0

)

2(1

)

0(

) = 0,

(17)

where the Cauchy stress tensor is defined as the derivative of strain energy with respect to strain and degraded due to the crack propagation. Similar to the strain energy, stresses are not divided into tension and compression parts.

Fig. 2. Schematic of a 2D domain. Dirichlet boundary conditions are imposed on the displacement and Neumann boundary conditions are imposed on the surface force. The phase-field is subject to Dirichlet and Neumann conditions = 1 on and ·n = 0 on . 3

Computational Materials Science 176 (2020) 109519

B. He, et al.

=

)2 + k ] C0: ,

= [(1

finite element method [53,60,38], multigrid with homogenization [67], and many more. The detailed description and comparison of each method are available through literature [60,25]. In this paper, we primarily use homogenization method for problems modeled as multiscale partial differential equations (PDEs). The homogenization theory is an analytical treatment for the multicalse PDEs. It has been studied for more than 30 years and still remains as an active research area. The homogenization theory deals with the macro description of micro heterogeneous systems and concentrates on the averaging process [7,18,41]. However, the analytical solution of multiscale solution is usually difficult to obtain in the case of an arbitrary domain and complex boundary conditions, so the numerical homogenization methods have been widely developed and proposed [1,2,3,26]. In this paper, we adopted a numerical approach within the heterogeneous multiscale method (HMM) framework, FE-HMM. HMM is originally proposed by E and Engquist [25,24] for general multiscale problems. As an efficient tool of gathering information from micro scale, HMM is featured with following advantages [2].

(18)

with boundary conditions:

· n = h on u = u 0 on ·n = 0 on

h u

(19)

.

The governing equation of phase-field in Eq. (17) demonstrates that the strain energy 0 determines the value of damage ϕ (the fracture propagates if the strain energy is larger than Gc ). However, because the strain energy can decrease, this approach does not guarantee the irreversibility of fracture, i.e. (20)

t.

t+ t

In this paper, the staggered approach is adopted to approximate the solution as suggested by Miehe et al. [54]. In this method, a history variable H (x, t ) is introduced to prevent reversibility issue of fracture described above. It is simply defined as the maximum of strain energy over time:

Ht +

=

t

Ht 0 ( ) if 0 ( ) Ht otherwise.

• works for different type of problems and operators, • allows for algorithm without specific assumption on the small scale, • allows for flexible discretization.

(21)

Thus the Karush-Kuhn-Tuker condition for loading and unloading process is satisfied [54]: 0

H

0,

Replacing (17) leads to:

Gc

1 l0

l

H 0

0,

H(

0

H ) = 0.

The FE-HMM is fully discrete analysis using finite element method (FEM) for multiscale PDEs with oscillating stiffness coefficients, established based on the framework of HMM, where the separation of both scales is honored [24]. It has been used to solve elliptic [3] and parabolic problems [4]. The comprehensive introduction of FE-HMM can be found in Abdulle’s paper [2]. The FE-HMM essentially roots in the mathematical homogenization theories based on the asymptotic expansion, so it has a sound theoretical basis. From the published documents, the FE-HMM framework has been used to solve problems including heat diffusion problem on rough surfaces [4], heat conduction [3], Darcy’s flow [2], linear elasticity [1], etc. In this paper, we primarily focus on the linear elastic solids. The modeling procedure of linear elasticity problems within the FE-HMM framework is summarized in the next section. In practice, this multi-scale method consists of solving two boundary value problems (BVP) with finite element method (FEM), the macro BVP and the micro BVP at each macro integration point. These BVPs are prescribed by the balance of linear momentum equations, as well as suitable boundary conditions. Instead of knowing the macro material properties, the stiffness of macro BVP is gathered from the solution of micro BVP. The concept of the two-scale homogenization is shown in Fig. 3.

(22)

by H in the governing equation of phase-field in Eq.

2(1

) H = 0.

(23)

H is now the driving term of the phase-field value. From Eq. (23), when H is sufficiently large so that the value of > 0 , the fracture starts initiating. The fracture propagates as the driving force H increases due to the external load defined in Eq. (17). Also, taking the maximum value of the strain energy guarantees the increase of H value over time and consequently guarantees the fracture irreversibility expressed in Eq. (20). For the discretized formulas in FEM implementation of the above phase-field fracture model, we heavily refer to the paper by MartnezPaneda et al. [51]. His supplementary document provides procedures using staggerd approach to decouple the governing Eq. (17). 3. Two-scale numerical homogenization In this section, a brief review of numerical homogenization will be presented, followed by a concise introduction of the Finite Element Heterogeneous Multiscale Method (FE-HMM), specific for the linear elastic solids. The computation of the homogenized elasticity tensor within the FE-HMM framework is discussed. Finally, a robust approach to obtain the unite strain setting in general quadrilateral elements is introduced.

3.2. Heterogeneous multi-scale finite element method for linear elastic solids The FE-HMM for linear elasticity problem has been discussed in details in the studies by Abdulle [1], Eidel and Fischer [26,28]. Abdulle presented an analysis of using FE-HMM to solve the linear elliptic

3.1. Numerical homogenization Multiscale modeling is increasingly important in many areas dealing with micro heterogeneous materials, including geosciences, biosciences, medicine, and material sciences, as the full system modeling is almost impossible [2]. The arising challenges include the effect of micro structure on the behavior of macro structures [26], distribution of flow pressure as it penetrates into porous media [2], prediction of bone fracture [8], etc. In past decades, a set of multiscale modeling methods (analytical and numerical) have been proposed and a few of them have been widely explored, including multiscale finite element method [37], FE2 method [44,43], heterogeneous multiscale method [24,25], general

Fig. 3. Schematic of numerical homogenization scheme solving two-scale BVPs. 4

Computational Materials Science 176 (2020) 109519

B. He, et al.

equation of elasticity problems. Eidel and Fischer described detailed numerical procedure of FE-HMM for homogenization of linear elastic solids with multiple numerical examples, and completed a comprehensive investigation on convergence and error analysis for different micro coupling conditions. The modeling procedure presented here is a concise summary of Eidel’s work. In the context of HMM, homogenization of linear elastic solids assembles the macro input data from a micro problem by averaging micro solutions [1]. The micro problem is defined as: xj

Cijlm

ul

Cijlm

ul

= fi in

xm

ui = u¯ i on

, P is the polynomial function greater than micro domain size, H space on macro elements, p is the polynomial function degree, K is the macro element. Obviously, R D ( , TH ) is subspace of V . The standard FEM form in Eq. (30) can be solved using numerical quadrature integration if elasticity tensor C0 (x) is given. For example, the bilinear form in Eq. (30) can be evaluated as shown in Eq. (32). However, for a heterogeneous material, C0 (x) is usually unknown. In the context of FE-HMM, the virtual work in bilinear form of Eq. (32) is computed from the known elasticity tensor C (x) at micro scale as shown in Eq. (33): Nqp

BH uH , v H =

D j

= t¯i on

N

(24)

Nqp

where indicates the characterized small length scale of heterogeneitydependent quantity, ijlm is the elasticity tensor for the micro domain, u¯ i is the prescribed displacement on Dirichlet boundary D , t¯i is the traction on the Neumann boundary N , fi is the body force, nj is the outward normal vector of the domain, . The corresponding homogenized macro problem is defines as:

K TH l = 1

xj

xm

0 Cijlm

ul0 xm

ui0 = < u¯ i > 0 Cijlm

ul0

j

xm

on

N

V

f· vdV +

VN

t¯ ·vdA

v

Nqp

=

V, (26)

V = v;v

v

BD

uH =

C0 (u0)· (v) dV =

V

f·v +

VN

t¯ · vdA

v

V. (29)

f· v H dV +

N

t¯ · v H dA

where the macro element space R

R

D(

, TH ) = {uH = u¯ ;

H1 ( );uH uH

K

D(

u

R

D

, TH ,

Nqp

= l=1

(30)

, TH ) is defined as: K

TH },

Kl

Kl

0 BTT I C (x K l ) BJ d

(Luh (I ) Kl )T C (x) LuhK(l J ) dV ,

(35)

Mmac

NIH dH I ,

uh =

Mmic i=1

Nih dih,

(36)

w Kl h (I ) T mic h (J ) (d ) K K d , K

(37)

(38)

where Nqp is the number of integration points at macro element K , T is the micro element, mic is the number of nodes in micro domain K , dh (I , xi) is the displacement vector after solving the micro minimization problem defined in Appendix A (Eq. (A.2)). The solution of the micro problem in Eq. (A.2) is discussed in Appendix A. e, mac After assembly keI,,mac can be exJ , the macro stiffness matrix k K pressed as:

D

(P q (K ))ndim ,

(34)

keI,,mac = BNe [NI , NJ ] J

If using the piecewise linear continuous FEM for both micro and macro problems, the homogenized problem in Eq. (29) can be written in following form: BH uH , v H =

,

where NIH and Nih are the macro and micro shape functions, Mmac and Mmic are the number of macro nodes in macro domain and micro nodes h in a micro domain, respectively, dH I and di are displacement vector for macro node I and micro node i, respectively. Replacing uh (I ) in Eq. (35) with expanded expression in Eq. (36), the stiffness matrix at macro node I , J can be calculated as:

(28)

V

w Kl

l=1

For macro BVP, the bilinear form is:

B0 u0, v =

(33)

M

B is the bilinear form on micro domain, x is the position vector, v is the test function, is the elastic strain tensor and can be expressed by using linear differential operator L with respect to displacement u as: (u ) = Lu .

,

H {Nih}i =mic 1 , the solution for macro problem, u , and solution for micro problem, uh can be expanded as:

(27)

= 0,

C (x) (u hK ): (v hK ) d

(32)

where keI,,mac is the elementary stiffness matrix at macro node I , J , LI is J the linear differential matrix of shape function, B is the strain and displacement matrix, BI = LNI . M Using the macro space basis {NIH } I =mac and the micro space basis 1

where V is the function space satisfied the Dirichlet boundary conditions:

)d ,

d

keI,,mac = BHe (NhI , NH J ) = J

l=1

H1 (

v H (x Kl)

·

where is the size of micro domain, . It is also noteworthy that FEHMM honors scale separation, so length scale ratio of micro domain 1. Eqs. (32) and and macro domain should be much less than one, / l (33) are the energy equivalence condition between macro and micro domains, which is equivalent to the Hill’s condition in FE2 method. The bilinear form in Eq. (32) in terms of shape functions lead to the stiffness matrix keI,,mac at macro node I , J : J

(25)

C (x) (u )· (v) dV =

K

1 1 , 2 2

K = x Kl +

where superscript “0” indicates the homogenized/macro quantity, the symbols <> and <> indicate the averaging operator for f, u¯ and t¯ . As discussed above, FE-HMM is a fully discrete analysis method and two BVPs at micro and macro scales will be solved using FEM. To solve these two BVPs, it is convenient to reformulate the macro and micro BVPs in Eqs. (24) and (25) in variational form: V

uH (x Kl)

where h indicates the micro element size, x Kl refers to the integration point in a macro element K , l indicates lth integration point, w Kl is the weight factor of quadrature integration, K is the micro sampling domain at integration point x l and is defined as:

D

on

1 K

w K l·

= < fi > in

= < t¯i >

B u,v =

w K l · C0 x k l K TH l = 1

(31)

superscript H indicates the macro element size, which usually is much 5

Computational Materials Science 176 (2020) 109519

B. He, et al.

Fig. 4. Schematic of superimposition of nodal unit displacement by Eidel [26] (Permission of using this figure in Eide's paper [26] was granted via Elsevier). Nqp

keK, mac =

w Kl T mic TK K k TK , K

l=1

wisely by evaluating the corresponding virtual work over the micro domain. In terms of the FE-HMM for linear elastic solids described in Section 2.2, the solution of micro problem for calculating the macro stiffness matrix is also available for the evaluation of numerical homogenized elasticity tensor. For the simplex, the micro displacement can be directly used to evaluate strain and elastic tensor coefficients, as the imposed unit displacement of shape function at macro nodes already build the unit strain state in the macro element. However, for quadrilateral elements, the combination of shape function is needed to create the unit state as discussed in following section.

(39)

where is the stiffness matrix of macro element K , is assembled stiffness matrix for micro domain K , TTK is the assembled nodal displacement of micro domain K and can be expressed as:

keK, mac

K mic k

(40)

TK = [[[dh (I , xi) ]i = 1, … , d ]I = 1, … , Nnode ]

Assembling elementwise, we can get the stiffness matrix for the macro problem defined in Eq. (30) and compute the macro solution.

keK, mac

3.3. Homogenized elasticity tensor

3.4. A robust method to generate unit strain states

In the FE-HMM discussed above, the stiffness of the macro problem is calculated from solution of the micro problem, so we can circumvent to compute stiffnes by using the constitutive model and homogenized elasticity tensor, which is usually unknown for heterogeneous materials. However, it is noteworthy that, instead of extracting macro stiffness from the micro solution, the homogenized elasticity tensor C0 can be approximated by the numerical approximate C0, h from the micro domain in FE-HMM framework as well. The detailed derivation of approximate tensor can be found in the literature [1]. The error analysis introduced by C is discussed in [1]. It is noteworthy that the derivation in the reference is based on the simplex elements. We recall the concise homogenization formulas as follows. For a given reference simplex, its vectorial basis can be defines as H H h N 0, j = (1 1 … d ) j , N i, j = ij, i , j = 1, …, d and ui, j is the solution of the micro problem in [1]. The numerical homogenized elasticity tensor can be computed as:

1 V

V

C (x ) (uih, j ): (ulh, m) d = h

1 Y

h

Y

It is noteworthy that computation method in above is limited to simplex element. For the quadrilateral element, the imposed unit displacement at macro nodes in FE-HMM cannot create a required unit strain state for numerical homogenization of elasticity tensor. Eidel in his paper [26] already pointed this problem out and gave the procedure to combine the displacements at macro nodes in order to make the unit strain state in a rectangular macro elements, simply superimposing the unit displacement at two nodes and multiplying them with the element length as illustrated in Fig. 4. However, the limitation in Eidle’s method is that the superimposition of unit displacement works only for rectangular elements. The proposed method cannot create the unit strain state in the general quadraliteral element as shown in Fig. 13(b). From Eq. (42), the unit strain state in the macro element is necessary for the numerical homogenization of elasticity tensor. To extend the homogenization method to general quadrilateral elements, a new method to generate unit strain state is proposed here. (See Fig. 5). In continuum mechanics [52,77], homogenized deformation of a continuous matter can be expressed as an affine map from reference configuration to current configuration:

h

C0, h (Ni, j ): (Nl, m ) d y

h

= C0, h (Ni, j ): (Nl, m)

(41) (42)

where V is the area or volume of micro domain, Y is the area or volume of the cell in the cell problem of asymptotic expansion homogenization theory [7]. In Eq. (42), a unit strain state is created in terms of the given shape functions in the reference simplex, as H

(Ni, j ) =

1 ei 2

ej + ej

ei , i, j = 1, …, d

where X is the material point position in reference configuration, x is the material position in current configuration, F is the deformation gradient, c is the translation vector. From the relation between elastic strain and displacement, the unit strain can be expressed as

(43)

=

where ei and ej are the space unit vectors, is the elastic strain tensor. Therefore, from Eq. (42), the coefficients of C0, h can be computed as followings: 0 ijlm

0, h ijlm

=

1 V

K

C (uih, j ): (ulh, m) d x ,

i, j , l, m = 1, …, d.

(45)

x = FX + c,

=

1 2

1 ei 2

ej + ej

u+(

ei

u)T

(44)

where main, is elasticity tensor in micro element, is the strain tensor, V is the volume of micro domain, and i , j, l, m refer to the index. In Eq. (42), by using the selected macro shape functions at macro nodes, a unit strain state is created in the macro element. This unit state allows us to compute the numerical homogenized tensor coefficient0 ijlm is the homogenized elasticity tensor, K is the micro do0, h ijlm is the homogenized elasticity tensor from micro domains, C

Fig. 5. (a) Rectangular element, (b) general quadrilateral element. 6

(46)

(47)

Computational Materials Science 176 (2020) 109519

B. He, et al.

Fig. 6. Four different micro samples: micro1 is used to verify the comparability between the pure macro model and homogenized model; micro 2–4 are used to explore the micro heterogeneity effect on macro fracture behavior.

=

1 F + FT 2

2I ,

micro pore network and fracture toughness of porous media is not considered in this paper because it is still under study. To avoid the pore effect, the fracture is assumed to be much larger than the pore size. For example, the fracture is in order of macro element size while pore size is usually 10 orders of magnitude smaller than macro element size.

(48)

where I is the identity matrix. From Eq. (48), one simple way to create the unit strain state is to impose the quadrilateral element with the following deformation. Note that the symmetry of F is not necessary.

F= 2 1 . 1 2

4.2. Verification between pure macro model and homogenization-based model

(49)

To illustrate the impacts of micro pore structures on fracture propagation of heterogeneous porous media, four different micro pore structures as illustrate in Fig. 6 have been selected. In this section, Micro 1 will be used to demonstrate the comparability of the homogenization-based fracture modeling framework with pure macro fracture modeling. The other three micro pore structures will be used in the next section to illustrate the impact of micro pore structures on the fracture behavior of porous media. To satisfy the scale separation requirement in FE-HMM, the size of each domain is designed sufficiently small so that ratio of micro and macro domain sizes is very small, /l 1. The dimension size of each micro domain is shown in Table 1. For the demonstration of comparability, the classical single-edge crack model is used to compare the simulation result of the pure macro model and homogenization-based fracture model. The geometry and loading conditions of the single-edge crack model is illustrated in Fig. 7. The bottom edge is fixed and the displacement is imposed on the top edge in vertical direction. For both the pure macro and the homogenized models, we used following material properties: Young’s Modulus E = 210 GPa, Poisson ratio = 0.3, critical energy release rate Gc = 2.7 N/mm, and characteristic length l 0 = 0.015 mm. The tensile load of u = 10−4 mm is applied at the top edge of the plate for first 300 steps. Then u = 10−5 mm is imposed for the rest of the computational steps. In our literature review, we noticed that pre-existing crack in the single-edge crack problem can be modeled either by detached nodes [64,51] or geometrical slot [54,85]. In our numerical studies, we created a pre-existing crack by defining detached nodes as well as slots with width w = 0.01 mm and w = 0.005 mm as shown in Fig. 7. Fig. 8 shows the reaction force and displacements. As one can see, cases with slot are in good agreement with Miehe’s result [54]. It is also shown that the geometrical slots led to lower fracture strength compared with the case with detached nodes (e.g., without slot). Also, our results indicates that wider slot leads to slightly lower fracture strength. Although the case with detached nodes overestimates the load-displacement curve comparing with Miehe’s work, we decided to move forward with this modeling method, so we do not have to deal with the uncertainty associated with the variation in the slot width.

With above deformation gradient, from Eq. (45), the current configuration can be determined if translation c is given. For the quadrilateral element, one node can be assumed to be fixed in the unit strain state, so that c can be computed using Eq. (45). After F and c are determined in Eq. (45), the displacement at each node of the quadrilateral element can be calculated as:

uI = x

X

(50)

One example of generating unit strain state for general quadrilateral element is provided in Appendix B. 4. Homogenization-based phase-field fracture modeling scheme and verification In this section, the framework of homogenization-base fracture modeling is discussed. Subsequently, a verification experiment is conducted to illustrated the alignment of solution between two-scale homogenization model and pure macro model. 4.1. Fracture modeling based on the two-scale homogenization As shown in Eqs. (13) and (23), in the phase-field fracture modeling framework, elasticity tensor directly relates to strain energy, and fracture propagation is driven by the evolution of energy term, H. As discussed in the introduction, for heterogeneous porous media, micro pore structure significantly impacts the elastic stiffness. However, homogenized elasticity tensor is usually unknown due to the oscillatory micro heterogeneity. The numerically-homogenized elasticity tensor within FE-HMM provide an opportunity to approximate the macro elasticity tensor without knowing the macro material constants. This advantage of FE-HMM has inspired us to bring homogenized elasticity tensor into phase-field fracture model and create a homogenizationbased fracture simulation to include the micro pore structure impacts on elastic stiffness as well as its further influence on fracture. In the modeling procedure, the homogenized tensor C0, h will be used to replace C0 in Eq. (17). The stiffness reduction due to fracture propagation is represented by a degradation function (Eq. (12)). Fully fractured elements ( = 1) will not contribute to the updated stiffness matrix. Therefore, fractured elements will be removed from computing stiffness. In this paper, our focus is on linear elastic porous media, so we assume that homogenized elasticity tensor will not change along the deformation and is approximated only at the beginning of the simulation. The study of nonlinear elasticity fracture will be presented in our future publication. It is worth mentioning that fracture modeling of porous media in this paper does not account for the pore skeleton effect on fracture, as that is out of the scope of this study [58]. Also, the relation between

Table 1 The dimensions of micro domains and pore sizes. Micro domain

Micro 1

Micro 2

Micro 3

Micro 4

Domain size (mm)

Length: 2E−4, Width: 2E−4 0 1E−4

Length: 2E−4, Width: 2E−4 19.7% 2E−5

Length: 2E−4, Width: 2E−4 20% 2E−5

Length: 2E−4, Width: 2E−4 32.2% 2E−5

Porosity (%) Element size (mm)

7

Computational Materials Science 176 (2020) 109519

B. He, et al.

Fig. 7. Schematic of single-edge crack model with three different pre-existing crack models.

Furthermore, a few research has shown the mesh-dependency of phase-field technique for fracture modeling [64,93]. For instance, Molnar et al. [64] showed that the reaction force is a function of mesh size h. To illustrate this dependency, we chose two different macro mesh sizes (e.g., h = 0.003 mm and h = 0.005 mm). These reaction force curves are plotted in Fig. 9a and the fracture propagation contour is present in Fig. 9b. From Fig. 9a, one can see that the reaction force of pure macro model and homogenization-based model are almost identical for both models regardless of mesh sizes. Additionally, Fig. 9b shows that the fracture pattern in both pure macro and homogenization-based models are also very similar. This example demonstrates that the homogenized simulation can produce sufficiently exact result as the pure macro simulation if no voids exist in the micro domain. In other words, the homogenized model is comparable to pure macro phase-field fracture model. A closer look at Fig. 9a reveals that the fracture strength is slightly higher for coarser mesh, indicating that different discritization influences the reaction force. This trend is consistent with what has been reported in the literature [93,64].

5. Case studies In this section, a set of numerical examples are tested using the proposed homogenization-based approach to illustrate the impact of micro heterogeneity on fracture behaviors in elastic porous media. 5.1. Single-edge crack plate We start with the classical single-edge crack model to investigate the effect of different micro pore structure on fracture behavior. Material properties, element size, and boundary condition remain same as the model with h = 0.003 mm used in the verification Section 4.2. The fracture propagation patterns in all two-scale homogenization cases (micro 2, micro 3, and micro 4) are almost same to the pure macro simulation as shown in Fig. 9b. However, the reaction force reduces with the increase of porosity, as shown in Fig. 10. Through comparison with the pure macro model, one can conclude that the forces reduce 15%, 13%, 30% with respect to the pure macro case, for single circular pore, single square pore, and five circular pores homogenized cases, Fig. 8. Plot of reaction force and displacement for the single-edge crack cases using detached nodes and slots with width w = 0.01 mm and w = 0.005 mm.

8

Computational Materials Science 176 (2020) 109519

B. He, et al.

Fig. 9. Comparison of (a) reaction force and (b) fracture propagation, for fully macro model and homogenized model. Solid lines are pure macro results while dashed line represents the two-scale results.

respectively. In this classical example, the fracture propagation path is not influenced by the micro pore structures as the energy release pattern is almost identical in this simple geometry even varying the underneath micro pore structures. However, there is a pronounced variation in the fracture strength as shown in Fig. 10. This is mainly due to the fact that the voids in the micro domains reduce the elastic stiffness. From Fig. 10a, one can see the peak reaction force in the homogenized model with five circular holes is about 500 N, which is much lower than 740 N for the base case (pure macro case). This is mainly because of an increase in the porosity. The peak force for both single square and single circular micro domain, are about 645 N and 624 N, respectively. They are lower than the pure macro model but higher than the model of five circular pores. However, it is interesting that the peak force for single square is slightly higher than the peak force of single circular micro domain, even though the porosity of single square is slightly higher as shown in Table 1. This might be caused by the shape effect in micro domain as the deformation of micro elements varies with different void shapes. Fig. 10 also shows that the softer model need more imposed loads to propagate the crack. For example, the pure macro model takes about 0.0065 mm to be fully broken while the single square and circular take about 0.008 mm to be fully broken and five circular model takes about 0.0095 mm to be fully broken. However, it is also interesting that the single square model take less tensile load than single circular to be fully broken although it shows softer behavior than single circular model. This difference may also be caused by the micro pore shape as the square shape leads to more deformation.

From this simple and classical fracture model, one can see that both the porosity and pore shape affect the macro fracture behavior in porous media. While the fracture strength is mainly related to the porosity, the pore shape influence the fracture pattern more. It will be demonstrated through following two complicated fracture models. 5.2. Asymmetrical double notches with holes The model of asymmetrical double notched plate with two holes has been widely investigated with different numerical methods, such as adaptive re-meshing FEM [11], numerical manifold method [92], and XFEM [80,81]. Although the experimental data is still not available, these different simulation results reached a good agreement that is summarized in [85]. Unfortunately, the corresponding phase-field modeling discussed in [85] cannot produce these consistent fracture patterns. In our phase-field simulation, the pure macro model shows a similar crack pattern as reported in the literature [11] (see Fig. 11b). In this section, we use this model to investigate the role of micro heterogeneity on the macro fracture behavior. The geometry and loading condition for this problem is shown in Fig. 11a. The bottom edge is fixed and the tensile load is applied on the top edge in Y direction. The material properties are: Young’s modulus E = 210 GPa, Poisson’s ratio = 0.3, energy release rate Gc = 1 N/mm. The mesh size is h = 0.08 mm, and the characteristic length is l 0 = 0.08 mm. The tensile load is u = 5 × 10−4 mm for first 100 steps and it is reduced to u = 2 × 105 mm for the remaining steps. The result of pure macro model is shown in Fig. 11b, while the result of homogenized-based cases are shown in Fig. 11c–e.

Fig. 10. Reaction force for pure macro model, homogenized model with micro 2–4. 9

Computational Materials Science 176 (2020) 109519

B. He, et al.

Fig. 11. (a) Schematic of geometry and loading condition for the model of plate with asymmetrical double notches with holes, (b) fracture patterns in pure macro model, (c)–(e) fracture pattern in the homogenized models based on micro domain 2–4 (as shown in Fig. 6).

From Fig. 10, one can see that the pattern in the pure macro model is smoother than two-scale cases as the elasticity tensor is assumed to be homogeneous. While the homogenized tensor based on the micro domains with pores may vary with the deformation of macro elements, these variations lead to the change of path for the energy release. As shown in Fig. 11c–e, the fracture pattern shows zigzagged patterns along the propagation path. Furthermore, we also notice that the fracture propagation of two crack paths are not fully symmetrical. Particularly for the homogenized model using micro case with square pore (e.g., Fig. 6), lower segment of crack shows the branching of cracks, which is not observed along the upper segment of crack. At this moment, the reason for this phenomena is not clear. However, the mesh-dependency of the phase-field fracture modeling and sensitivity of homogenized stiffness on the micro pore shapes might be the two contributing factors.

inhomogenized energy evolution, which may cause the variation of crack coalescence. In this section, we adopted a fracture model with three pre-existing flaws to test the micro heterogeneity effect on the fracture coalesce [51]. The Fig. 12a shows the geometry and load conditions. The imposed tensile load is applied to left and right edges. Unlike the model described in the literature [51], the pre-imposed damage around each flaw was not considered in our model since our main purpose is to investigate the effect of micro heterogeneity [51]. The material properties are: Young’s modulus E = 200 GPa, Poisson’s ratio = 0.3, energy release rate Gc = 0.09 N/mm. The mesh size is h = 0.1 mm, and characteristic length is l 0 = 0.6 mm. The imposed displacement is 0.12 mm. The imposed tensile load for the first 200 steps is u = 7.5 × 10 4 mm and u = 1.5 × 10 4 mm for the remaining steps. Fracture pattern of the pure macro model in Fig. 12b matches the result reported in the literature [51]. However, fracture patterns of homogenized models show deviations from the pure macro model even though the over all patterns are similar and the same macro meshes are used. Similar to the previous case discussed in Section 5.2, all fracture paths in homogenized models show zigzagged patterns as illustrated in Fig. 12. In particular, fracture pattern in the case with square pore does not flow smoothly as the pure macro model. These variations might be caused by the strain energy changes due to homogenized stiffness tensor. The fracture patterns in homogenization-based models with

5.3. Crack coalescence Crack coalescence is an important issue when multiple damages exist, because it often leads to the failure of structure. As proposed by Sutula, et al. from Griffith’s law, the crack growth in linear elastic solids is a energy minimization problem [79,80]. They also pointed out that minimum energy is the fracture growth criteria [81]. Therefore, the crack coalescence usually grows along the path that minimizes the total energy [30]. However, the variation of stiffness in porous media leads to 10

Computational Materials Science 176 (2020) 109519

B. He, et al.

Fig. 12. (a) Schematic of geometry and loading conditions for the model with multiple pre-existing damages, (b) fracture pattern of pure macro model, (c) - (e) fracture patterns of homogenized models based on micro domain 2, 3, and 4 (as shown in Fig. 6).

micro 2 and 4 are quite similar as they have similar pore shapes. Also, the fracture pattern is closer to the pure macro model than the case with square pore. It should be noted that the fracture patterns still show zigzagged shape, as highlighted in the sub-figures. As discussed above, this is caused by the variation of homogenized elasticity tensor based on micro domains. We also observe that the width of fracture diffusive zones reduce in the homogenized models as illustrated in the highlighted sub-figures (Fig. 12c-e). The measured width are: ld = 6.6 mm , for pure macro model, ld = 5.7 mm for single circular pore, ld = 4.5 mm for square pore, ld = 6 mm for multiple circular pores. These numerical comparison of fracture behavior during the crack coalesce shows that both the propagation path and fracture diffusive zones will be affected by the micro heterogeneity. However, these observed findings still need to be further validated upon availability of experimental data.

6. Conclusion In this paper, we introduce a computational homogenization-based phase-field fracture modeling scheme to include the micro heterogeneity effect on the fracture behaviors in the linear elastic porous media. Through this approach, the oscillation of micro heterogeneity in porous media can be considered in fracture modeling using the homogenized stiffness information without need to know the macro stiffness. In the proposed framework, the homogenized elasticity tensor is estimated with the FE-HMM and used to calculate the strain energy within phase-field fracture modeling framework where the energy drives the evolution of fracture. To allow the implementation of the above modeling scheme for general quadrilateral elements, a robust approach to generate the unit strain state is developed. From the benchmark tests of a set of fracture models, our fully macro models 11

Computational Materials Science 176 (2020) 109519

B. He, et al.

achieved a good agreement with the reference literature. We used there different micro domains to homogenize the elasticity tensor and update strain energy in phase-field fracture modeling. Comparison with the pure macro result revealed that the micro heterogeneity, not only affects the fracture strength but also the fracture propagation patterns. Even though these numerical result still need experimental validation, these finding are still providing valuable insights for the study of fracture behavior in porous media. In our future work, a more complex framework will be developed for nonlinear porous media to investigate the fracture behaviors under compound loads.

Validation, Visualization, Writing - original draft, Writing - review & editing. Louis Schuler: Software, Writing - original draft, Writing review & editing. Pania Newell: Conceptualization, Funding acquisition, Supervision, Project administration, Writing - review & editing, Visualization.

Data availability

Acknowledgements

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

This work was partially supported by Dr. Newell’s start-up fund provided by the University of Utah and partially supported by the Multi-Scale Fluid-Solid Interactions in Architected and Natural Materials (MUSE) Project, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award #DE-SC0019285.

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.

CRediT authorship contribution statement Bang He: Methodology, Investigation, Software, Formal analysis, Appendix A. Solution of micro problem in FE-HMM

The discrete version of the micro problem, as described by Eq. (26) in each micro domain K at each integration point l = 1, …, Nqp can be reformulated as: Find u hK such that the bilinear form and the preset micro-macro coupling condition are satisfied.

(uhK

)

H ulin ,K

R q (K , Th)

B K (u hK , w hK )·= u hK

K

(A.1)

(u hK )· (w hK ) dV = 0

S q (K , Th )

(A.2)

where S q (K , Th) = {uh W (K );uh (P q (T ))ndim , T defined at the quadrature point x K as: H H ulin , K = u (x Kl) + (x

For the coupling of

u KH (I , xi) =

Mmic

x Kl)· H ulin , Kl

Th}. W (K ) is the boundary condition for micro domain K , The linearization of uH is

uH (x Kl). and

u hK(I , xi)

(A.3) in one micro domain, they are expanded in the same space

NihMmici = 1

H (I , xi) Nmh , K dm

u hK(I , xi) =

, Th) : (A.4)

m=1 Mmic

of

S q (K

h (I , xi) Nmh , K dm .

(A.5)

m=1

The solution of micro problem described in Eq. (A.2) is essential to approximate the homogenized macro stiffness matrix in Eq. (39). To obtain the micro displacement dI , xi , the micro problem can be reformulated as the a energy minimization problem: K

(Lu hK(I , xi) )T C (Lu hK(I , xi) ) dV

min (A.6)

BK

H R q (Kl , Th ) , where BK is the virtual work of internal force on micro over all function over functions u hKl S1 (K , Th ) and satisfies (u hKl ulin , Kl ) domain K . For the purpose of implementing different boundary conditions on the micro domain, the Lagrange multipliers method is employed to solve the minimization problem in Equation A.6 [3]. Subsequently, we yield the total energy fictional:

L dh (I , xi) ,

(I , xi)

=

1 h (I , xi) T mic h (I , xi) + (d ) K K d 2

(I , xi) T G

dh (I , xi)

dH (I , xi)

(A.7)

where G is kinematics coupling conditions for micro and macro domains [26]. For 3 D problems, it reads

G=

b1 0 0 … b Mmic 0 0 0 b1 0 … 0 b Mmic 0 0 b1 … 0 0 B Mmic ¯ G

(A.8)

12

Computational Materials Science 176 (2020) 109519

B. He, et al.

¯ is the L non-redundant periodic conditions conditions for micro nodes on the opposite edge/face of micro domain. The detail discussion where G about L non-redundant boundary condition refers to [3]. The coefficient bm for first d rows in G is the normalization conditions for the macro/micro coupling [26]. They are associated with the d Lagrange multipliers. coefficients bm can be calculated using the kinematics coupling condition defined in Eq. (A.2). The detailed derivation refers to [3,26]. Here, we just give the final calculation formulas.

bm = = =

T 1 4 1 4

h Nm dV =

1 wT

T

T for quadrilateral elements and q = 1, T for hexahedral elements and q = 1.

(A.9)

is the stiffness matrix of micro domain K . It is assembled from the micro elemental stiffness matrix standard Gauss quadrature integration method:

K mic K

nqp

e, mic Kmn

eT wk Bm

k,

k, k

e C Bm

k,

k, k

det J

k,

k, k

e, mic Kmn .

e, mic Kmn

is evaluated using the

,

(A.10)

k=1 e = LNm, J ( k, k , k ) is the Jacobian of isoparmateric mapping, k, k , k are the coordinates variable in reference cube element. where Bm The first variation of L with respect to variables dh (I , xi) and h (I , xi) yields the matrix equation for micro problems:

K mic GT K

0 dh (I , xi ) = I , xi GdH (I , xi ) G 0 for I = 1, …, Nnode ; i = 1, …, d,

(A.11)

where dh (I , xi) is the displacements in the micro domain for the imposed unit macro displacement at macro node I, I , xi is the associated Lagrange multiplier. Solving Eq. (A.11), we can obtain micro displacements dh (I , xi) . As the stiffness matrix K mic K and coupling boundary matrix G are constant for each macro unit displacement, the micro and macro displacements , dH (I , xi) dH , and dh (I , xi) T. can be easily augmented to full matrix, I , xi Appendix B. Example of unit state strain for a general quadrilateral element To illustrate the calculation process to generate the unit strain state in one general quadrilateral element, one example is presented here. The geometry and coordinates of the element is shown in Fig. 13. From the discussion in Section 3, the deformation gradient is assumed to be:

F= 2 1 . 1 2

(B.1)

In order to get translation vector c , let node 1 to be fixed during the deformation, so c can be calculate as:

C = X1

(B.2)

FX1

= 2 , 2

(B.3)

where X1 is the coordinates of node 1. Using Eq. (45), the coordinates of each node in unit strain state will be:

x1 =

1 , 1

x2 =

2

1

,

x3 =

5 , 6.5

x4 =

2.5 . 4.5

(B.4)

Therefore, displacement for each node can be calculated following Eq. (50) as:

u1 =

0 , 0

u2 = 1 , 1

u3 =

4.5 , 4.5

u4 =

3 , 3

(B.5)

where uI , I = 1, 2, 3, 4 is the displacement of each node.

Fig. 13. Example of one general quadrilateral element.

13

Computational Materials Science 176 (2020) 109519

B. He, et al.

References

[37] T. Hou, X.-H. Wu, Z. Cai, Convergence of multi-scale finite element method for elliptic problems with rapidly oscillating coefficients, Math. Comput. 68 (1999) 913–943. [38] T.J. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Eng. 127 (1) (1995) 387–401. [39] H. Jelitto, G. Schneider, A geometric model for the fracture toughness of porous materials, Acta Mater. 151 (2018) 443–453. [40] Y. Jiang, H. Lian, V.P. Nguyen, W. Liang, Propagation behavior of hydraulic fracture across the coal-rock interface under different interfacial friction coefficients and a new prediction model, J. Natural Gas Sci. Eng. 68 (2019) 102894. [41] V. Jikov, S. Kozlov, O. Oleinik, Homogenization of Differential Operatorsand Integral Functionals, Springer-Verlag, Heidelberg, Berlin, 1994. [42] N. Khalili, Two-phase fluid flow through fractured porous media with deformable matrix, Water Resour. Res. 44 (5) (2008). [43] V. Kouznetsova, M. Geers, W. Brekelmans, Multi-scale first-order and second-order computational homogenization of microstructures towards continua 1, 371–386. [44] V. Kouznetsova, M. Geers, W. Brekelmans, Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy, Comput. Methods Appl. Mech. Eng. 193 (2004) 5525–5550. [45] J. Kováčik, Correlation between Young’s modulus and porosity in porous materials, J. Mater. Sci. Lett. 18 (1999) 1007–1010. [46] J. Kubik, Macrodescription OP micropore structure in regard to fluid flow through porous media, in: K. Unger, J. Rouquerol, K. Sing, H. Kral (Eds.), Characterization of Porous Solids. Vol. 39 of Studies in Surface Science and Catalysis, Elsevier 1988, pp. 345–353. [47] H. Laubie, S. Monfared, F. Radjaï, R. Pellenq, F.-J. Ulm, Disorder-induced stiffness degradation of highly disordered porous materials, J. Mech. Phys. Solids 106 (2017) 207–228. [48] D. Leguillon, R. Piat, Fracture of porous materials – influence of the pore size, Eng. Fract. Mech. 75 (7) (2008) 1840–1853. [49] X. Lu, M. Viljanen, 10 – fibrous insulation materials in building engineering applications, in: R. Fangueiro (Ed.), Fibrous and Composite Materials for Civil Engineering Applications, Woodhead Publishing Series in Textiles. Woodhead Publishing, 2011, pp. 271–305. [50] B. Markicevic, N. Djilali, Two-scale modeling in porous media: relative permeability predictions, Phys. Fluids 18 (3) (2006) 33–101. [51] E. Martínez-Pañeda, A. Golahmar, C.F. Niordson, A phase field formulation for hydrogen assisted cracking, Comput. Methods Appl. Mech. Eng. 342 (2018) 742–761. [52] G. Mase, R. Smelser, G. Mase, Continuum Mechanics for Engineers, CRC Press, Boca Raton, 2010. [53] A. Matache, I. Babuka, C. Schwab, Generalized p-fem in homogenization, Numer. Math. 86 (2000) 319–375. [54] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits, Comput. Methods Appl. Mech. Eng. 199 (45) (2010) 2765–2778. [55] C. Miehe, S. Mauthe, Phase field modeling of fracture in multi-physics problems. Part III. Crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media, Comput. Methods Appl. Mech. Eng. 304 (2016) 619–655. [56] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations 83 (2010) 1273–1311. [57] E. Mikaeili, B. Schrefler, XFEM, strong discontinuities and second-order work in shear band modeling of saturated porous media, Acta Geotech. 13 (6) (2018) 1249–1264. [58] A. Mikelic, F. Wheeler, T. Wick, A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium, SIAM Multiscale Model. Simul. 13 (2015) 367–398. [59] E. Milanese, T.D. Cao, L. Simoni, B.A. Schrefler, Fracturing in dry and saturated porous media, in: E. Oñate, D. Peric, E. de Souza Neto, M. Chiumenti (Eds.), Advances in Computational Plasticity: A Book in Honour of d. Roger j, Owen, Springer International Publishing, Cham, 2018, pp. 265–288. [60] P. Ming, X. Yue, Numerical methods for multiscale elliptic problems, J. Comput. Phys. 214 (1) (2006) 421–445. [61] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Numer. Meth. Eng. 46 (1) (1999) 131–150. [62] T. Mohammadnejad, A. Khoei, An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model, Finite Elem. Anal. Des. 73 (2013) 77–95. [63] T. Mohammadnejad, A.R. Khoei, Hydro-mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method, Int. J. Numer. Anal. Meth. Geomech. 37 (10) (2013) 1247–1279. [64] G. Molnár, A. Gravouil, 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture, Finite Elem. Anal. Des. 130 (2017) 27–38. [65] A. Munjiza, K. Andrews, J. White, Combined single and smeared crack model in combined finite-discrete element analysis, Int. J. Numer. Meth. Eng. 44 (1) (1999) 41–57. [66] A. Munjiza, D. Owen, N. Bicanic, A combined finite-discrete element method in transient dynamics of fracturing solids, Eng. Comput. 12 (2) (1995) 145–174. [67] N. Neuss, W. Jäger, G. Wittum, Homogenization and multigrid, Computing 66 (2001) 1–26. [68] V.P. Nguyen, H. Lian, T. Rabczuk, S. Bordas, Modelling hydraulic fractures in

[1] A. Abdulle, Analysis of a heterogeneous multiscale FEM for problems in elasticity, Math. Models Methods Appl. Sci. 16 (2006) 615–635. [2] A. Abdulle, The finite element heterogeneous multiscale method: a computational strategy for multiscale PDEs, GAKUTO Int. Ser. Math. Sci. Appl. 31 (2009) 135–184. [3] A. Abdulle, A. Nonnenmacher, A short and versatile finite element multiscale code for homogenization problems, Comput. Methods Appl. Mech. Eng. 198 (37) (2009) 2839–2859. [4] A. Abdulle, C. Schwab, Heterogeneous multiscale fem for diffusion problems on rough surfaces. SIAM Journal on, Multiscale Model. Simul. (2005) 3. [5] T. Belytschko, T. Black, Elastic crack growth in finite elements with minimal remeshing, Int. J. Numer. Meth. Eng. 45 (5) (1999) 601–620. [6] T. Belytschko, S.P. Xiao, Coupling Methods for Continuum Model with Molecular Model 1 (1) (2003) 12. [7] A. Bensoussan, J.-L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structure, North Holland, Amsterdam, 1978. [8] P. Bhattacharya, M. Viceconti, Multiscale modeling methods in biomechanics, Wiley Interdiscip. Rev.: Syst. Biol. Med. 9 (3) (2017) 1375. [9] M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J. Hughes, C.M. Landis, A phase-field description of dynamic brittle fracture, Comput. Methods Appl. Mech. Eng. 217–220 (2012) 77–95. [10] R. Borst, J. Réthoré, M.-A. Abellan, A numerical approach for arbitrary cracks in a fluid-saturated medium, Arch. Appl. Mech. 75 (2006) 595–606. [11] P. Bouchard, F. Bay, Y. Chastel, Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria, Comput. Methods Appl. Mech. Eng. 192 (35) (2003) 3887–3908. [12] B. Bourdin, G. Francfort, J. Marigo, The variational approach to fracture, J. Elast. 91 (1) (2008) 5–148. [13] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, Berlin, Heidelberg, 1991. [14] P. Budarapu, R. Gracie, S.-W. Yang, X. Zhuang, T. Rabczuk, Efficient coarse graining in multiscale modeling of fracture, Theoret. Appl. Fract. Mech. 69 (2014) 126–143. [15] T.D. Cao, F. Hussain, B.A. Schrefler, Porous media fracturing dynamics: stepwise crack advancement and fluid pressure oscillations, J. Mech. Phys. Solids 111 (2018) 113–133. [16] P. Chakraborty, Y. Zhang, M.R. Tonks, Multi-scale modeling of microstructure dependent intergranular brittle fracture using a quantitative phase-field based method, Comput. Mater. Sci. 113 (2016) 38–52. [17] Y. Chen, H. Lian, W. Liang, J. Yang, V.P. Nguyen, S.P. Bordas, The influence of fracture geometry variation on non-Darcy flow in fractures under confining stresses, Int. J. Rock Mech. Min. Sci. 113 (2019) 59–71. [18] D. Cioranescu, P. Donato, An Introduction to Homogenization, Oxford (1999) UniversityPress. [19] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, first ed., Wiley Publishing, 2009. [20] J. Crawford, The relationship between structure and the hydraulic conductivity of soil, Eur. J. Soil Sci. 45 (4) (1994) 493–502. [21] R. de Borst, Chapter 1 – introduction, in: R. de Borst (Ed.), Computational Methods for Fracture in Porous Media, Elsevier, 2018, pp. 1–12. [22] R. de Borst, V. Verhoosel, Gradient damage vs phase-field approaches for fracture: similarities and differences, Comput. Methods Appl. Mech. Eng. 312 (2016) 78–94. [23] T. Dutta, S. Tarafdar, Fractal pore structure of sedimentary rocks: simulation by ballistic deposition, J. Geophys. Res.: Solid Earth 108 (B2) (2003). [24] E. Weinan, B. Engquist, Z. Huang, Heterogeneous multiscale method: a general methodology for multiscale modeling, Phys. Rev. B 67 (9) (2003) 092101. [25] E.W. Engquist, B. Li, X. Ren, W. Vanden-Eijnden E., Heterogeneous multiscale methods: a review, Commun. Comput. Phys. 2 (3) (2007) 367–450. [26] B. Eidel, A. Fischer, The heterogeneous multiscale finite element method for the homogenization of linear elastic solids and a comparison with the FE2 method, Comput. Methods Appl. Mech. Eng. 329 (2018) 332–368. [27] Y. Feng, K.E. Gray, XFEM-based cohesive zone approach for modeling near-wellbore hydraulic fracture complexity, Acta Geotech. 14 (2) (2019) 377–402. [28] A. Fischer, B. Eidel, Convergence and error analysis of FE-HMM/FE2 for energetically consistent micro-coupling conditions in linear elastic solids, Eur. J. Mech.- A (2019) 77. [29] G. Francfort, J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46 (1998) 1319–1342. [30] G. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46 (1998). [31] A.A. Griffith, VI. G.I.T., The phenomena of rupture and flow in solids, Philos. Trans. R. Soc. London Ser. A 221 (1921) 163–198. [32] P. Gupta, C.A. Duarte, Simulation of non-planar three-dimensional hydraulic fracture propagation, Int. J. Numer. Anal. Meth. Geomech. 38 (13) (2014) 1397–1430. [33] Y. Heider, B. Markert, A phase-field modeling approach of hydraulic fracture in saturated porous media, Mech. Res. Commun. 80 (2017) 38–46. [34] Hirshikesh, C. Jansari, K. Kannan, R. Annabattula, S. Natarajan, Adaptive phase field method for quasi-static brittle fracture using a recovery based error indicator and quadtree decomposition, Eng. Fract. Mech. 220 (2019) 106599. [35] Hirshikesh, S. Natarajan, R. Annabattula, A FEniCS implementation of the phase field method for quasi-static brittle fracture, Front. Struct. Civil Eng. 13 (2) (2019) 380–396. [36] A. Pramod, R. Annabattula, E. Ooi, C. Song, S. Natarajan, Adaptive phase-field modeling of brittle fracture using the scaled boundary finite element method, Comput. Methods Appl. Mech. Eng. 355 (2019) 284–307.

14

Computational Materials Science 176 (2020) 109519

B. He, et al. porous media using flow cohesive interface elements, Eng. Geol. 225 (2017) 68–82. [69] S. Ozaki, Y. Aoki, T. Osada, K. Takeo, W. Nakao, Finite element analysis of fracture statistics of ceramics: effects of grain size and pore size distributions, J. Am. Ceram. Soc. 101 (7) (2018) 3191–3204. [70] A. Paluszny, S. Salimzadeh, R.W. Zimmerman, Chapter 1 – finite-element modeling of the growth and interaction of hydraulic fractures in poroelastic rock formations, in: Y.-S. Wu (Ed.), Hydraulic Fracture Modeling, Gulf Professional Publishing, 2018, pp. 1–19. [71] R.H.J. Peerlings, R. de Borst, W.A.M. Brekelmans, M.G.D. Geers, Gradient-enhanced damage modelling of concrete fracture, Mech. Cohesive-Frictional Mater. 3 (4) (1998) 323–342. [72] F. Pesavento, B.A. Schrefler, G. Sciumè, Multiphase flow in deforming porous media: a review, Arch. Comput. Methods Eng. 24 (2) (2017) 423–448. [73] J.M. Sargado, E. Keilegavlen, I. Berre, J.M. Nordbotten, High-accuracy phase-field models for brittle fracture based on a new family of degradation functions, J. Mech. Phys. Solids 111 (2018) 458–489. [74] B.A. Schrefler, L. Simoni, A unified approach to the analysis of saturated-unsaturated elastoplastic porous media, in: Numerical Methods in Geomechanics, sixth ed., vol. 1, CRC Press, United States, 2017, pp. 205–212. [75] M. Sheng, G. Li, D. Sutula, S. Tian, S.P. Bordas, XFEM modeling of multistage hydraulic fracturing in anisotropic shale formations, J. Petrol. Sci. Eng. 162 (2018) 801–812. [76] L. Simoni, B.A. Schrefler, Chapter four – multi field simulation of fracture, Vol. 47 of Advances in Applied Mechanics. Elsevier (2014) 367–519. [77] A. Spencer, Continuum Mechanics, Courier Corporation, 2012. [78] B.-L. Su, C. Sanchez, X.-Y. Yang, Hierarchically Structured Porous Materials: From Nanoscience to Catalysis, Separation, Optics, Energy, and Life Science, Wiley-VCH, 2012. [79] D. Sutula, P. Kerfriden, T. van Dam, S.P. Bordas, Minimum energy multiple crack propagation. Part I: theory and state of the art review, Eng. Fract. Mech. 191 (2018) 205–224. [80] D. Sutula, P. Kerfriden, T. van Dam, S.P. Bordas, Minimum energy multiple crack propagation. Part-II: discrete solution with XFEM, Eng. Fract. Mech. 191 (2018) 225–256.

[81] D. Sutula, P. Kerfriden, T. van Dam, S.P. Bordas, Minimum energy multiple crack propagation. Part III: XFEM computer implementation and applications, Eng. Fract. Mech. 191 (2018) 257–276. [82] H. Talebi, M. Silani, S. Bordas, P. Kerfriden, T. Rabczuk, Molecular dynamics/xfem coupling by a three-dimensional extended bridging domain with applications to dynamic brittle fracture, Int. J. Multiscale Comput. Eng. 11 (6) (2013) 527–541. [83] C. Torres-Sanchez, J. McLaughlin, R. Bonallo, Effect of pore size, morphology and orientation on the bulk stiffness of a porous Ti35Nb4Sn alloy, J. Mater. Eng. Perform. 27 (6) (2018) 2899–2909. [84] Z.A. Wilson, C.M. Landis, Phase-field modeling of hydraulic fracture, J. Mech. Phys. Solids 96 (2016) 264–290. [85] J.-Y. Wu, V.P. Nguyen, C. Thanh Nguyen, D. Sutula, S. Bordas, S. Sinaie, Phase field modelling of fracture, Adv. Appl. Mech. (2019) 53. [86] S. Xiao, T. Belytschko, A bridging domain method for coupling continua with molecular dynamics, Comput. Methods Appl. Mech. Eng. 193 (17–20) (2004) 1645–1669. [87] J. Yang, H. Lian, W. Liang, V.P. Nguyen, S.P. Bordas, Model I cohesive zone models of different rank coals, Int. J. Rock Mech. Min. Sci. 115 (2019) 145–156. [88] K. Yoshioka, B. Bourdin, A variational hydraulic fracturing model coupled to a reservoir simulator, Int. J. Rock Mech. Min. Sci. 88 (2016) 137–150. [89] B.D. Zdravkov, J. Cermák, M. Sefara, J. Janku, Pore classification in the characterization of porous materials: a perspective, Cent. Eur. J. Chem. 5 (2) (2007) 385–395. [90] H. Zhang, X. Ge, H. Ye, Effective thermal conductivity of two-scale porous media, Appl. Phys. Lett. 89 (2006). [91] X. Zhang, W. Sloan, C. Vignes, D. Sheng, A modifcation of the phase-feld model for mixed mode crack propagation in rock-like materials, Comput. Methods Appl. Mech. Engrg. 322 (2017) 123–136. [92] H. Zheng, F. Liu, X. Du, Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method, Comput. Methods Appl. Mech. Eng. 295 (2015) 150–171. [93] S. Zhou, X. Zhuang, H. Zhu, T. Rabczuk, Phase field modelling of crack propagation, branching and coalescence in rocks, Theoret. Appl. Fract. Mech. 96 (2018) 174–192.

15