Regularization of continuum damage mechanics models for 3-D brittle materials using implicit gradient enhancement

Computers and Geotechnics 122 (2020) 103505

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Regularization of continuum damage mechanics models for 3-D brittle materials using implicit gradient enhancement

T

S. Mondala,b, , L.M. Olsen-Kettleb, L. Grossa ⁎

a b

The University of Queensland, School of Earth and Environmental Sciences, St Lucia, QLD 4072, Australia Swinburne University of Technology, School of Science, Mathematics Department, Hawthorn, Vic 3122, Australia

ARTICLE INFO

ABSTRACT

Keywords: Damage and fracture Brittle materials Mesh independence Finite element method Split-operator technique Non-local model

A well-known problem which arises in local damage models is that they suffer from mesh dependence and strain localization. In this study, we employ the non-local implicit gradient damage formulation for mesh sensitivity analysis of brittle materials. The basic concept of the non-local implicit gradient model is that the strain at a given point depends not only on the strain at that point but also on the nearby strain field. The local equivalent strain is replaced with a non-local equivalent strain, which is constructed by solving a system of partial differential equations (PDEs). We solve this system of PDEs using finite element method (FEM) and split-operator method which has been applied to 3-D microplane damage models. We extend this approach to continuum damage models and apply it to brittle material for three different loading experiments. In these examples, we demonstrate that the implementation of the implicit gradient damage model using the split-operator method can successfully be applied to complex, 3-D geometries requiring large-scale unstructured FEM meshes, and eliminate mesh dependence. In not all cases the simulations were able to obtain the desired fracture widths as this would require meshes at even smaller resolutions which is beyond our current computing capacity.

1. Introduction It is widely acknowledged in the literature that local continuum damage models lead to localization issues [32,23]. This means that the predictions of the finite element simulations change significantly when changing the size and orientation of the mesh used for discretization. The growth of damage tends to localize in the smallest band that can be captured by the spatial discretization resulting in damage localization at the smallest length scale set by the length of a single element in the finite element mesh. Consequently, the solutions for a case with a progressively damaged zone converge to a solution with a localization zone of zero volume when successive mesh refinements are applied. In the end, the numerical prediction fails to converge to a physically meaningful result. Therefore, the global response of the system shows a strong dependence on the spatial discretization. The causes and effects of localization have been discussed in detail in Borst et al. [4]. The introduction of non-locality in the constitutive model is one of the popular methods to regularize the localization of damage in order to eliminate mesh dependence. In a non-local damage formulation, the state of stress (or strain) at a point does not only depend on the state of strain and its history at that particular point but also depends on the state of the media in a finite neighbourhood of that point. There are two ⁎

families of non-local formulations, namely, on the one hand, integral formulations using a kernel to regularize the stress or strain [1,25,2] and on the other hand, gradient-based formulations [8,5,27,42,21]. Non-local approaches introduce a length scale to the mathematical representation of the physical system and thus prevent the localization problem. The concept of non-locality for elastic deformations was introduced by Kröner [20] and later extended to hardening plasticity by Mühlhaus and Alfantis [27]. A breakage diffusion model for strength softening rock is studied by Mühlhaus and Gross [28]. Lasry et al. [22] introduced the concept of a non-local constitutive model as a localization limiter for a strain-softening material. These localization limiters [3,19] force the damage to grow in a zone with a finite width that is independent of the spatial finite element discretization. There are two forms of non-local gradient formulation namely explicit and implicit, where the non-local equivalent strain is related to the local equivalent strain through partial differential equations [44]. Explicit gradient enhancement is be considered as weakly non-local only, as it fails to regularize the stress (or strain) under some circumstances [55]. On the other hand, implicit gradient enhancement has the advantage of being strongly non-local. It can be shown that they are largely equivalent to the integral type [43], see also in Section 2. The idea of implicit gradient enhancement is to introduce a second

Corresponding author at: The University of Queensland, School of Earth and Environmental Sciences, St Lucia, QLD 4072, Australia. E-mail address: [email protected] (S. Mondal).

https://doi.org/10.1016/j.compgeo.2020.103505 Received 31 July 2019; Received in revised form 19 February 2020; Accepted 19 February 2020 0266-352X/ © 2020 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

differential equation in addition to the quasi-static equilibrium equation [43]. The equation takes the form of a Helmholtz equation to calculate the non-local field as the counterpart of a local equivalent strain or another local internal variable. The damage evolution law now depends on the regularized equivalent strain rather than the local version. The regularized equivalent strain is obtained by solving the Helmholtz equation. Consequently, the momentum and Helmholtz equations form a non-linear, coupled system of partial differential equations (PDEs). We use the split-operator method [13,15] and small time steps to solve the two equations sequentially by alternating between them while increasing loading. Section 2.4 details the equations solved using the split-operator approach. The split-operator algorithm is a robust scheme commonly used in implementations of gradient damage [26] and phase-field models [29]. The method is also used in 2D nonlinear elastic brittle damage by Pires-Dominuges et al. [46] and for 3-D non-local microplane damage models by Zreid et al. [55]. In this paper we extend this implementation to the general 3-D case without any assumption on damage pattern. Our implementation support massive parallel computation which is crucial for the solution of problems with very large finite element meshes required to resolve the targeted spatial length scales of the damage zone. The implementation is applied to the case of a 3-D rectangular rock specimens with a preexisting surface flaw under uniaxial compression. Local damage models have been investigated to simulate this experimental setup using unstructured meshes in 3-D where results have extensively been compared with a rich set of experimental observations [31]. In order to be able to represent the flaw the finite element method with an unstructured mesh is applied. The rock heterogeneity is modeled through random sampling at a fixed grid resolution which is then interpolated to the simulation meshes with varying mesh sizes. The allows an investigation on mesh dependence for the damage model results without changing material properties for meshes with different element sizes. To our knowledge, this paper presents the first non-local damage model applied to this experimental set-up. The paper is organized as follows: In Section 2 we present the ingredients of the non-local model including the governing equations and the damage model based on the implicit gradient approach and the provision of heterogeneous material properties. Section 2.4 introduces the operator splitting and its implementation. Results of simulation of the damage propagation for three different application cases are presented in Section 3. The effect of localization length scale is discussed in Section 4. In Section 5 the mesh sensitivity analysis for all three application cases is discussed. An evaluation of the results is given in Section 6.

(1)

=0

ij, j

The lower index , j refers to the derivative with respect to spatial direction x j . Summation over double lower indices (in this case j) is applied. For the isotropic, linear elastic case the components of stress tensor ij are given by Hooke’s law as ij

(2)

= Cijkl kl,

where kl is the strain tensor. The stiffness tensor Cijkl for an isotropic medium takes the form

Cijkl =

E (1 + )(1

2 )

ij kl

+

E ( 2(1 + )

ik jl

+

jk il ).

In Eq. (2) summation over double index kl applies. The symbol ij is Kronecker delta symbol, is the Poisson’s ratio and E is the Young’s modulus. The strain tensor is given in the geometrical linearized form as kl

=

1 (uk, l + ul, k ) 2

(3)

where u = (uk ) is the displacement vector. In this paper we investigate specimens with plane top and bottom surface which are perpendicular to the x3 axis, see Fig. 1. Loading is applied at the top surface for the rectangular and cubic specimens, and along the side for the L-shaped specimen. All the specimens were fixed at the bottom surface; more details are presented in Section 3. At all other boundary surfaces natural boundary conditions are applied, i.e., the normal stress is set to zero: ij nj

(4)

=0

where n = (ni ) is outer normal field at the domain surface. This boundary condition can easily be modified to include a confining pressure. We solved the momentum Eq. (1) for the displacement vector u = (uk ) using the FEM [54]. 2.2. Damage evolution The damage model investigated is based on a single condition using a scalar equivalent strain which is defined in terms of the strain components and is defined later. We use an isotropic damage model [31] in which the elastic Young’s modulus E is reduced as damage of the rock progresses. Young’s modulus E of the damaged material is defined as

E = E0 (1

(5)

D)

where E0 is Young’s modulus prior to loading. The scalar damage variable D (0 D 1) can be interpreted as an isotropic stiffness degradation factor. The material is totally undamaged for D = 0 , and D = 1 corresponds to complete loss of material coherence. When D = 1 the damaged Young’s modulus is equal to zero according to Eq. (5), however this causes numerical issues and instead a small number, i.e. E = 10 5 Pa, was specified for the limit for the damaged elastic modulus [49]. The assumption of isotropic damage is often sufficient for isotropic rocks [49,53]. The implementation of isotropic damage models is less complex and is computationally inexpensive because the damage variable is a scalar instead of a fourth order tensor [38,37,36]. The damage variable D may increase under the influence of loading. Whether or not damage growth occurs is decided on the basis of a damage loading function related to the equivalent strain components [44]:

2. Methodology We consider a rock specimen which is incrementally loaded until failure. Fig. 1 shows the geometry of the three application cases investigated later in this paper, see Section 3. In each loading step the new displacement vector (ui ) is calculated and the resulting damage to the rock based on the new strain tensor ( ij ) is determined. Material properties are updated according to the damage model before the next loading step is applied. This iterative process has been implemented using esys.escript which is a module in python for solving partial differential equations (PDEs) using the finite element method (FEM) [47]. In all simulations presented in this work linear tetrahedral elements have been used, see Fig. 1. In this section we describe this procedure and the underlying model in more details.

f( , )=

2.1. Governing equations

,

(6)

where is the threshold variable for the initiation of damage. As long as f < 0 there is no damage and the material behaviour remains linearly elastic. Damage is initiated when f = 0 is reached at which point the damage variable D starts to increase. The evolution of damage variable D is governed by the threshold

In each loading step the new stress tensor ( ij ) needs to fulfill the conservation of momentum equation which is given in Einstein notation as 2

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 1. The rock specimen (a) a tetrahedral mesh element (b) cubic specimen with weak zone marked in read, (c) block specimen with pre-existing flaw and (d) Lshaped specimen.

variable

(

where the update of

) = 0,

0 and

is defined by three conditions

0.

(7)

where is the change of between two loading steps. This condition never decreases but is inexpresses mathematically the fact that creased to the value of if the equivalent strain exceeds its value. The update process guarantees that f ( , ) 0 at any time. The threshold variable is initialized with 0 > 0 . As a consequence damage will initiate once the equivalent strain reaches the initial strength 0 . The damage variable D is an explicit function of this threshold variable ; D = D ( ) . We use the modified power law first proposed by Geers [14] to describe the damage evolution as we found this gave realistic fracture patterns:

D=

1 1,

( 0) (

c c

0

) , if if

<

c, c,

.

(8)

In this model the parameter influences the initial rate of damage growth and determines the final softening stage, close to complete failure [14]. An example for damage evolution law is plotted in Fig. 2 for parameters given in Table 1. As initially 0 = the damage variable D is initially zero for undamaged material initially but starts to increase

Fig. 2. Damage versus equivalent strain with parameters given in Table 1.

3

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

the average of over the volume occupied by the micro-structure at x where the size of the micro-structure is determined by the radius over which is averaged. This averaged quantity ¯ is defined in each material point x by the convolution of with a smoothing kernel g:

Table 1 Parameters for the damage calculation. Parameters

Values

Damage growth influence parameter ( ) Softening parameter ( ) Threshold damage parameter ( 0 )

0.01 1.6

Ultimate damage parameter ( c )

¯ (x) =

2.0 × 10

4

1.0 × 10

3

1 1 I1 + 2 ) 2

2 (1

( (1

1) 2 2 I1 2 )2

12 J2 , (1 + ) 2

kk

and

J2 =

1 2 I1 6

1 2

ij ij.

(9)

=

0 E0 ,

=

=

0

E0.

= ,

t

=1

(x). 4

+

1 2!

2

(x).(2)

2

+

1 3!

3

(x).(3)

3

+

1 4!

4

(15)

+…

where n and .(n) are the nth order gradient operator and the nth order inner product respectively; n designates the n factor dyadic product … . Substitution of (15) into (14) and the assumption of spatial isotropy and symmetry of the kernel leads to

¯=

+c

2

+d

4

(16)

+…

where denotes the Laplacian operator and c, d and higher order coefficients are determined as moments of the weight kernel g. Notice that terms of odd powers of vanish due to the symmetry of the kernel g. Differentiating Eq. (16) twice and reordering yields 2

2

(11)

2¯

=

c

4

d

6

+…

(17)

Substitution of this identity into Eq. (16) leads to

¯

c

2¯

=

+ (d

c 2)

4

+ …..

(18)

Neglecting higher-order terms in expression (18), simplifies the equation to

¯

(12)

c

2¯

(19)

=

> 0 is the second moment of the kernel The gradient parameter, c = g. It controls the width of the local micro-structure and as a consequence the size of the damage zone. The value l is the localization length. A larger value of l produces a larger inelastic zone. The non-local model forces the damage to diffuse in an area defined by the localization length scale l. To ensure that the mesh can represent this process it is necessary that the mesh elements are fine enough to accurately represent this localization length scale. We make an assumption that we need at least two elements inside the localization length scale in order to eliminate any mesh dependence and solve the Helmholtz equation accurately. Fracture patterns from the same loading experiment may be fitted to different length scales l which indicate that this length is not just a material parameter that is correlated to material properties only [40]. It may not be constant and can be a function depending on stress and strain fields [34]. It typically decreases at larger strains. However, this does not hinder the practical application of this models, since a good calibration of the length scale may be obtained by conducting a fracture test of notched specimens as it had been shown for the microplane model [6]. The non-local strain ¯ is not given explicitly in terms of and its derivatives, but as the solution of a boundary value problem consisting of the Helmholtz Eq. (19) and appropriate boundary conditions. As it can be solved efficiently using FEM the Helmholtz equation provides a convenient tool to construct the non-local equivalent strain ¯ from the non-local equivalent strain on the same FEM mesh on which the 1 2 l 2

using the equivalent strain . When combining this condition (12) for compressive strength with condition (11) for tensile strength one obtains c

(x) + (x).(4)

(10)

where 1, 2 and 3 are the principal stresses. From Eq. (11), 0 can be seen as the tensile strain strength. Similarly for purely compressive loading ( 1 = 2 = 0 > 3 ) one sets 3 = c with compressive strength c > 0 at the point of damage initiation. This translates into the condition c

g ( )d

(x + )

The modified von Mises equivalent strain in Eq. (9) depends on the intermediate principal stress as well as the principal maximum and minimum stress. The modified von Mises equivalent strain criteria has also been recently used in several papers [33,39,48] for both 2-D and 31) is a given parameter which D studies of damage. In Eq. (9), ( conducts the sensitivity in compression relative to that in tension. In terms of principle stresses 1 > 2 > 3 tensile strength t is reached for purely tensile loading ( 1 > 2 = 3 = 0 ) when 1 = t . Using the equivalent strain the damage initiation condition (6) is expressed in the form t

1 vol ( )

in which vol ( ) is the volume of set and denotes the relative position vector pointing to the infinitesimal volume d . A potential choice for g is a Gaussian kernel where the width is related to size of the microstructure. The calculation of the convolution (14) for each FEM element in a large 3-D mesh is computationally inefficient. We use the gradient formulation instead. The gradient formulation for the non-local equivalent strain ¯ is derived from approach (14) in the following way. The local equivalent strain (x + ) is expanded into a Taylor series at x [22,41,27]:

with I1 the first invariant of the strain tensor and J2 the second invariant of the deviatoric strain tensor. These are given by:

I1 =

g ( ) (x + ) d , where

(14)

once damage is initiated and 0 < . The case D = 1 characterizes the complete loss of material consistency which is reached once has reached ultimate damage parameter c . It remains to define the equivalent strain . We follow Peerling [40] who used a strain quantity based on the von-Mises stress and Hooke’s law:

=

1 vol ( )

(13)

This shows that is the ratio the compressive and tensile strength. The detailed calculation of this relationship is described in the study by Mondal [30]. 2.3. Non-local implicit gradient model We extend the local damage model described in the previous section to be reformulated using a nonlocal description for the equivalent strain. The idea is to replace the local equivalent strain by a smoothed, non-local version, where the damage initiation criterion and evolution law remain unchanged. In this non-local damage model, is related to the non-local equivalent strain ¯ as a weighted volume average of this (local) equivalent strain [10,45]. The growth of the damage variable associated with a mesh element depends on the averaged strain across neighboring mesh elements. In terms of the elasticity based damage model the above reasoning means that the growth of damage at a point x is no longer governed by the local equivalent strain at x alone, but by 4

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

momentum Eq. (1) is solved. The resulting gradient damage formulation will be referred to as the implicit gradient formulation. The treatment of ¯ is as an independent variable, which has to satisfy the partial differential Eq. (19). A difficulty of the gradient models is the requirement to specify a suitable boundary condition in order to complete the Helmholtz Eq. (19) as a boundary value problem with a unique solution. From a mathematical point of view, it is necessary to specify either ¯ or its normal derivative at the boundaries. The simple natural boundary condition

¯, i ni = 0

heterogeneity [52] under laboratory conditions where a small spatial correlation length can be assumed. The Weibull probability distribution is given as:

f (u w ) =

(20)

2.4. Split-operator method implementation The model couples two partial differential equations, i.e. the momentum equation and the non-local Helmholtz equation defined in Eq. (21) and (22), respectively:

=0

(21)

¯

c ¯, ii =

(22)

m 1

·exp

uw u0

m

(23)

where u w is a given material parameter to be modelled (here the initial Young’s modulus E0 ) and m is the shape parameter which defines the degree of material homogeneity and is thus referred to as the homogeneity index. u 0 is the scale parameter of the distribution function which is related to the mean bulk value of the distribution. Only Young’s modulus is sampled from a Weibull distribution. High values for m correspond to an almost homogeneous structure with a uniform Young’s modulus, whereas a heterogeneous structure with a broad distribution of local Young’s modulus is associated with a relatively low m value [51]. Fig. 4 depicts the probability density function of Weibull distribution for various values of the shape parameter m. In response to an externally applied stress, the elastically heterogeneous material develops a spatially heterogeneous stress field and whenever the stress field in an element satisfies the failure criterion, the Young’s modulus begins to decrease in magnitude for that element. The Weibull distribution can be sampled over each element of the FEM mesh. That means the length scale of the heterogeneity is the element size and becomes dependent on the mesh size. The method creates realizations that are unsuitable for investigating mesh dependence as the local material parameters change when the mesh is refined. To address this problem, we have implemented the rock heterogeneity using an interpolation table based on a fixed cubic grid which is sampled from a Weibull distribution. The same interpolation table is used for all simulations when refining the mesh. In the tests presented in this paper a grid size of 5 mm3 is used.

is adopted in this paper as in other works [22,27].

ij, j

m uw u 0 u0

The local equivalent strain depends on the gradient of displacement as solution of the first PDE (21) while the elastic tensor is depends on the non-local equivalent strain ¯ via the damage evolution law. As a consequence Eqs. (21) and (22) are strongly coupled. Instead of solving the system in its fully coupled form at each loading step we apply the split-operator approach which consists of splitting the coupled problem into two sub-problems, namely (i) the quasi-static equilibrium Eq. (21) with a fixed equivalent strain ¯ for an increased loading and (ii) the Helmholtz problem (22) on a local equivalent strain calculated in (i). This method does not consider the full coupling [13,15] of the quasi-static equilibrium equation and the Helmholtz equation for the non-local equivalent strain. This simplification is made to reduce the computing time for 3-D specimens. The split-operator algorithm is a robust scheme often used in implementations of gradient damage [26] and phase-field models [29]. A similar strategy introducing an extrapolation of the threshold variable based on the two previous increments was proposed for cases with strong discontinuity [35,7,46]. We solve each equation using the FEM where we use the same mesh for both equations. For small localization length, l, the mesh needs to be very fine as at least two or more elements within l need be used to resolve the localization length scale and to obtain an accurate average value of the non-local strain. The split-operator method allows us to couple the gradient enhanced damage framework to the existing non-local finite element local damage code used in the Mondal et al. [31] in a transparent and convenient manner. The detailed solution procedure is shown in the Fig. 3 (Flow Chart). Two symmetric linear PDEs are solved sequentially in each loading step. When damage initiation has been detected ¯ 0 damage parameter D is updated based on the modified power law as discussed in Section 2.2. As the damage threshold is increasing the value for damage D may not be reduced between load increments. Finally, Young’s modulus is updated with new damage variable D and a new load increment is applied. The iterative process is terminated when the set maximum loading is reached. It is pointed out that the two solution fields – displacement and non-local equivalent strain – are slightly inconsistent but for small load increments these differences will be small and therefore quite acceptable and controllable.

2.6. Implementation The implementation of the model is based on esys.escript which is a python module for solving partial differential equations using spatial discretization techniques such as the finite element method [47,18]. It is designed to describe numerical models in the language of PDEs while using computational components implemented in C and C++ to achieve high performance for the time-intensive numerical calculations. A particular strength of esys.escript is that it supports the parallelized solution of PDEs using domain composition without additional programming work on the scripts developed and tested on a desktop computer. Parallel execution is supported through both MPI for distributed memory and OpenMP for shared memory including a hybrid mode. The efficiency of numerical simulations implemented using esys.escript when running script across large numbers of processors with large-scale FEM meshes has been demonstrated in various applications [9,17]. We use this feature of esys.escript to accelerate simulations when solving models based on very large FEM meshes required to resolve the localization length scale l for the nonlocal implicit gradient damage model on 3-D domains. The results shown in this paper are using esys.escript version 5.3 [12]. To solve the governing Eq. (1) we have written a python script using esys.escript as the PDE solver. The script has four major steps. Firstly, we define the domain representing the specimen to be loaded. For the simulation three different shaped domain ( ) are considered as shown in Fig. (1). Secondly, the PDE to be solved in each loading step is defined. Thirdly, the coefficients of the PDE based on the damage variable obtained in the previous iteration step are determined. Then the PDE is solved and the damage is updated where failed elements are detected. Finally, the termination condition is checked. For loading of the samples investigated in this paper the incremental loading is terminated once a specified displacement limit at the surface has been reached. The computational work flow of esys.escript is described in Fig. (5).

2.5. Rock heterogeneity A realistic rock heterogeneity in brittle material is modelled by sampling rock properties from a Weibull distribution. It is a continuous probability distribution and commonly used to model rock 5

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 3. Computational Work Flow for Non-Local Model.

3. Damage pattern of local and non-local model

homogeneous numerical specimen was chosen with a shape parameter m = 50 to sample the initial Young’s modulus. However, for the Lshaped specimen, Young’s modulus was set to constant value to be consistent with the previous published numerical simulations by Winkler et al. [50]. We use the modified power damage law discussed in Section 2.2 with the damage evolution law plotted in Fig. 2. As the non-local model forces the damage to diffuse in an area defined by the localization length scale l, we need to resolve the localization length scale. We make the somehow optimistic assumption that at least two elements per length scale l are sufficient to eliminate mesh dependence. All meshes are generated using the 3-D finite element mesh generator Gmsh [16]. The plots of distribution of the damage parameter D presented in following sections are showing the distribution after the maximum total displacement 1 mm was applied to the specimen.

We apply the non-local damage code to three specimens with different shapes: a cubic specimen with a weak zone in Section 3.1, a homogeneous L-shaped specimen in Section 3.2 and a rectangular specimen with a pre-existing weak flaw in Section 3.3. Materials parameters; in particular Young’s modulus; are chosen according to the values published in the relevant experimental studies [24,50]. For the three cases studied in this paper namely a cubic, L-Shaped and rectangular (with flaw) shaped specimen we have considered comparable material properties, peak stress and damage patterns to the published experimental papers by Viso et al. [11], Winkler et al. [50] and Lu et al. [24] respectively. Some validation with another numerical local damage model has already been published in Mondal et al. [31]. To match the experimental results for the cubic and rectangular specimen a fairly 6

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 4. Probability density function of Weibull distribution for different values of shape parameter m.

Fig. 6. Cubic rock specimen (a) with weak zone (b). The cube has the edge length of lzh = 100 mm. The weak zone colored in blue has a cubical shape with edge length 5 mm and is placed at the center of surface of the cube.

clamped boundary conditions by prescribing displacement in all direction at the bottom x3 = 0 and top x3 = lzh surface. To implement uniaxial compression along the x3 axis this is given in the form

u1 = u2 = u3 = 0 at x3 = 0 u1 = u2 = 0 and u3 = uBC at x3 = lzh,

(24)

where uBC (0.004 mm) is the vertical displacement increment enforced at the top of the specimen for uniaxial compression and the total displacement at the top surface increases with time. On all other faces no loading is applied which means that there is zero normal stress at the boundary of the specimens on these faces. We set the material and damage evolution parameters as listed in Table 1 and Table 2 respectively. The mesh is a very fine mesh with a total of about 40 million elements with a maximum edge length of 1 mm for the elements. The same mesh is used for both the non-local and local damage model in Fig. 7 (a) and (b). Fig. 7 plots the distribution of the damage variable D of cubic specimen with a weak zone under uniaxial compression for both local (Fig. 7 (a)) and non-local (Fig. 7 (b)) damage models. The localization length scale l is 2 mm for non-local damage in Fig. 7 (b). The maximum element size is 1 mm for both simulations. The non-local model forces the damage to diffuse at a scale determined by the localization length Table 2 Material properties for the cubic specimen with weak zone. Note that the Young’s modulus outside the weak zone gives the value of the scale parameter of the Weibull distribution used to model rock heterogeneity on a 5 mm3 grid. Fig. 5. Computational work flow for esys.escript.

3.1. Cubic specimen with weak zone In the first test the specimen has a cubic shape with edge length lzh = 100 mm. Uniaxial compression is applied along the direction x3 . Fig. 6 (a) shows the location of the weak zone with dimensions (5 mm3) in the middle of the front x2 x3 planar surface. The shape of the weak zone is shown in Fig. 6 (b). The initial Young’s modulus of the specimen is sampled from the Weibull distribution with scale parameter 36.0 GPa and shape parameter m = 50.0 as described in Section 2.5. In the weak zone the Young’s modulus is set to 18.0 GPa. We apply the following 7

Parameters

Values

Height of the specimen (lzh) Young’s modulus outside the weak zone Young’s modulus inside the weak zone Poisson’s ratio outside the weak zone ( ) Poisson’s ratio inside the weak zone Homogeneity index (m) Scale parameter (u0 ) Interpolation grid size Tensile strength ( t ) Compressive strength ( c ) Displacement increment Localization length scale (l)

100 mm 36.0 GPa 18.0 GPa 0.18 0.26 50.0 36.0 GPa 5 mm3 250.0 MPa 841.0 MPa 0.004 mm 2 mm

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 7. The damage plot of the cubic specimen with weak zone under uni-axial compression as shown in Fig. 6 using the local damage model (a) and the non-local damage model (b) with length scale parameter l = 2 mm. The color represents the value of the damage variable D. The material properties are shown in Table 2.

scale. The size of the damaged zone and the fracture width is determined by the localization length scale, and for this reason, the fracture width is larger for the non-local damage model than the local damage model. The damage propagates along shear bands emanating from the weak zone of the cubic specimen. The non-local model shows damage evolution and fracture propagation at a larger spatial scale in comparison to the local model because it has an inherent length scale prescribed by the model through the localization length scale. The fracture width is larger than the localization length in the non-local model as the damage diffuses over a scale of 2 5 times the localization length. A similar effect of localization length scale is shown for microplane damage models using an implicit gradient enhancement by Zreid and Kaliske [55]. It was not possible to reduce the localization length scale further in the non-local model because of computational resource limits for the simulations, as we require at least two finite element mesh elements to resolve the localization length.

localization length. Both the damage patterns produced by the local and non-local model are similar shape to experimental results [50] sketched in Fig. 8 (b). Again it was not possible due to computing resource limitations to reduce the localization length scale in the non-local model any further to more accurately reproduce the experimental results and fracture widths at a smaller length scale. 3.3. Block specimen with pre-existing weak flaw To replicate the experimental investigations of Lu et al. [24] we simulate a block-shaped sandstone specimen with a pre-existing flaw in the centre of the front x2 x3 planar surface of the specimen undergoing uniaxial compression. The block has dimensions 20 × 70 × 120 mm3, see Fig. 10. The flaw has a depth of 8 mm, width of 1 mm and a length of 30 mm and is positioned as angle of 75° against the x2 axis. The clamped displacement boundary conditions (24) similar to the cubic specimen are used for this block shaped specimen. A displacement increment of 0.002 mm per loading step is applied at the top face of the specimen at x3 = lzh until it reaches the final displacement (1 mm). The mechanical parameters chosen for the simulation are listed in Table 1 and Table 4. The heterogeneity of the sandstone is modeled by sampling the Young’s modulus outside the flaw from a Weibull distribution with shape parameter m = 50.0 over a regular grid of cell size 5 mm3 as described in Section 2.5. A very fine mesh is considered to calculate the damage evolution where the maximum length of an element is 1 mm and contains a total of 55.3 million elements. As observed in the cases presented previously the damage evolution and fracture propagation produced by the local model as shown in Fig. 11 (a) is very similar to the experimental results [24] (Fig. 7). Mondal et al. [31] (Fig. 8) were also able to reproduce these experiments using a local damage model based the unified strength theory. In contrast to the failure model presented here one of the main differences with the local damage model evolution laws employed in this work is that it uses separate criteria for compressive and tensile damage initiation and evolution. The fracture width emanating from the tips of the pre-existing fracture in the nonlocal model is more than double the fracture width of the local model in Fig. 11 and the experiment Lu et al. [24]. However, the fracture width in the region containing the pre-existing fracture for the non-local model is almost ten times bigger than the local model in Fig. 11 and the experiment Lu et al. [24]. This large fracture width is produced with the non-local model because there is a localization length scale present which diffuses the damage in a large area corresponding to the magnitude of the localization length. We were unable to reduce this localization length further due to memory limitations in the mesh

3.2. Homogeneous L-shaped specimen In this section we simulate a set up that was investigated in the experimental study by Winkler et al. [50]. In contrast to the other two specimens (cubic and rectangular prism) we consider a homogeneous (initially constant Young’s modulus without Weibull distribution) specimen for this case. The geometric properties and shape of the L-shaped specimen along with the boundary conditions are shown in Fig. 8 (a). The edge lengths at the top and in the recess are 500 mm and 250 mm respectively. The thickness of the sample is 100 mm. The lower horizontal edge of the vertical leg is fixed for all the directions (marked red in Fig. 8 (a)) and the vertical surface at x2 = ly is subjected to a vertical surface load (marked green in Fig. 8 (a)). The displacement boundary conditions are:

u3 = uBC at x2 = ly , u1 = u2 = u3 = 0 at x3 = 0.

(25)

The dimensions and properties of the concrete for the L-shaped specimen are given in Table 1 and Table 3. The localization length scale is 2 mm whereas the maximum element size is 1 mm for this L-shaped specimen. Total number of elements are 83.9 million for this simulation. Fig. 9 illustrates the damage for the local (a) and non-local model (b). The non-local damage in Fig. 9 (b) depicts the same pattern of damage as the local model in Fig. 9 (a). However, the fracture width in the nonlocal model is larger than in the local model. This larger damaged zone in the non-local model is expected as the localization length scale forces the damage to diffuse in a certain area which is defined by that

8

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 8. L-shaped specimen with boundary conditions and (a) and experimental results (b) redrawn from Winkler et al. [50].

mentioned already, at a large scale.

Table 3 Material properties for the L-shaped specimen. Parameters

Values

Long edge of the specimen (ly) Short edge of the specimen Scale parameter Young’s modulus (E0 ) Poisson’s ratio ( ) Tensile strength ( t ) Compressive strength ( c ) Displacement increment Localization length scale (l)

500 mm 250 mm 36.0 GPa 0.18 250.0 MPa 841.0 MPa 0.002 mm 2 mm

3.4. Summary For this Section 3 we conclude that damage patterns as indicated by the distribution of the damage parameters are qualitatively similar for both local and non-local models. As expected for the non-local model the fractured zone is smeared out over a larger region. In the next Section we will demonstrate that the size of this zone is controlled by the length parameter l. The damage pattern of the L-Shaped specimen in Fig. 9 shows qualitatively the same damaged zone as the experimental study by Winkler et al. [50]. Similarly, the damage pattern of the block specimen with a pre-existing flaw shows a qualitatively similar damaged zone to the pattern seen in experimental study by Lu et al. [24]. The results of the local model exhibit details of the damage pattern at a finer scale. This scale is set by the local element size of the FEM mesh and as a consequence the damage patterns are strongly mesh

generation. Because of computing resource limitations, the smallest localization length scale we could use was 2 mm for the non-local damage model presented in Fig. 11 (b). The same pattern of damage is observed in the non-local damage model as for the local model but, as

Fig. 9. The fracture propagation of L-shaped specimen using the local model (a) and non-local model (b) with length scale parameter l = 2 mm. Geometrical set-up is shown in Fig. 8. Material properties are listed in Table 3. 9

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

challenging to run simulations with a mesh resolution sufficient to resolve damage pattern on a realistic localization length scale. The limiting factor for the simulations is the FEM mesh generation. We have used high-performance supercomputing facilities with up to 15 compute nodes and with a total of 300 cores and were able to provide FEM solutions for very large meshes. Larger problems could be solved with our FEM solver using more cores. However, we were unable to generate 3-D meshes with more elements to resolve narrower fracture widths, due to memory limitations in the mesh generator. 4. Effect of localization length scale In the non-local model, the fracture width depends on the value of the non-local parameter c which is proportional to the square of localization length l as defined in (19). So the localization length is crucial for controlling the fracture width as the damage diffuses at the scale of the localization length. It also places limits on the maximum mesh element size that can be used so that the fracture width is resolved accurately. We make the assumption that we need at least two elements inside the localization length scale area to resolve the diffusion of damage at the scale of the localization length. This means that the mesh element size needs to be at least half the localization length scale. As in the worst case sharp changes in the damage parameter appear over a range that is 2–3 times larger than the value of localization length l this restriction on the element size in regions where damage initialization is expected sufficiently conservative in order to be able to resolve the damage patterns in the non-local model. Fig. 12 plots the damage evolution of the L-shaped specimen for different values for the localization length scales l. All simulations use meshes with a maximum element size of 1 mm. The results clearly show the effect of the localization length on the fracture width produced. In Fig. 12 (d) the localization length is 8 mm and the damaged zone is very wide. As the localization length is decreased the fracture width also decreases accordingly as demonstrated in (c), (b) and (a) of Fig. 12. For this application, it would be desirable to employ a smaller localization length in order to reproduce the experimental results of Winkler et al. [50] more accurately. However, computational resource limitations meant that we could not employ localization lengths less than 1 mm. Fig. 13 plots the stress-strain curve of the L-shaped specimen for different localization length scales l with a fine mesh with element size is 1 mm. Larger values of the localization length result in a stronger specimen which fails at a larger loads. This is caused by the redistribution of strain peaks to neighbouring locations in the averaging process. The decimation and redistribution process leads to later initiation of damage and hence a stronger specimen. However, for larger values of the localization length scale the failure threshold c is reached over a much larger area resulting in larger region of total damage (D = 1) as shown in Fig. 12 and at the same time faster deterioration of the specimen integrity once failure has reached. Thus can be identified in Fig. 12 through the steeper gradient in the stress-strain curve for larger values of l after the peak stress has passed.

Fig. 10. Rock specimen with pre-existing surface flaw (a). The flaw is disc shaped flaw (b) of a given width which penetrates the surface of the specimen at an angle. Table 4 Material properties for the block specimen with a weak flaw. Notice that the Young’s modulus outside the weak zone gives the value of the scale parameter of the Weibull distribution used to model rock heterogeneity on a 5 mm3 grid. Parameters

Values

Height of the specimen (lzh) Young’s modulus outside the weak zone Young’s modulus inside the weak zone

120 mm 13.5 GPa 1.0 × 10 5 Pa 0.26 50.0 13.5 GPa 5 mm3 250.0 MPa 841.0 MPa 0.002 mm 2 mm

Poisson’s ratio ( ) Homogeneity index (m) Scale parameter (u0 ) Interpolation grid size Tensile strength ( t ) Compressive strength ( c ) Displacement increment Localization length scale (l)

5. Mesh sensitivity analysis The aim of this section is to analyze the effect of the mesh size on the results for the non-local damage model and to show that the nonlocal model is, in fact, mesh independent. For comparison we also plot the damage plots for the local damage model as well, to show the problem of mesh dependency arising in local models. To control the computational cost and execution time for the simulations a relatively large localization length of l = 8 mm for the non-local model has been chosen. The stress-strain curves are inspected for the mesh dependence analysis. The maximum mesh element size in the meshes tested ranges from 1 to 8 mm for the three different shaped specimens investigated

dependent. In Section 5 we will demonstrate that the results of the nonlocal model converge and become mesh-independent when the mesh is refined. This solution is controlled by the value of localization length scale l, and is mesh independent. However, it is pointed out that the results shown in this section already reveal that it is computationally

10

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 11. The fracture propagation of block specimen with a pre-existing flaw with flaw angle 75°, flaw depth 8 mm, flaw length 30 mm. Distribution of damage parameter D for the local model (a) and the non-local model (b) with length scale parameter l = 2 mm. Material parameters are specified in Table 4.

Fig. 12. The damage evolution of an L-shaped specimen for different localization length scales, (a) localization length 2 mm, (b) localization length 4 mm, (c) localization length 6 mm and (d) localization length 8 mm. The element size is 1 mm for all simulations.

11

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Table 5 Peak stress of different mesh element size in local and non-local model (localization length scale 8 mm) for different specimens.

Fig. 13. The stress-strain curve for different localization length scales l with mesh element size 1 mm for L-Shaped specimen. The corresponding damage plots are shown in Fig. 12.

Cubic specimen with weak zone

Local model

Non-local model

1 mm 2 mm 4 mm

3.21e+ 07 Pa 3.25 e+ 07 Pa 3.30 e+ 07 Pa

3.34 e+ 07 Pa 3.34 e+ 07 Pa 3.35 e+ 07 Pa

L-Shaped Specimen

Local model

Non-local model

1 2 4 8

3.19 e+ 05 Pa 3.55 e+ 05 Pa 4.34 e+ 05 Pa 5.42 e+ 05 Pa

6.76 e+ 6.76 e+ 7.00 e+ 7.46 e+

Block specimen with single flaw

Local model

Non-local model

1 mm 2 mm 4 mm

1.25 e+ 07 Pa 1.30 e+ 07 Pa 1.34 e+ 07 Pa

1.33 e+ 07 Pa 1.34 e+ 07 Pa 1.35 e+ 07 Pa

mm mm mm mm

05 05 05 05

Pa Pa Pa Pa

Fig. 15. L-shaped specimen: The stress-strain curve for varying mesh resolutions for (a) the local model and (b) the non-local model with localization length of l = 8 mm.

Fig. 14. Cubic specimen with weak zone: The stress-strain curve for varying mesh resolutions for (a) the local model and (b) the non-local model with a localization length of l = 8 mm.

million respectively. The peak stresses for this simulation as for the other simulations shown in the section are summarized in Table 5. A comparison between the local model in Fig. 14 (a) and non-local model in Fig. 14 (b) reveals that the non-local model predicts a slightly stronger material in the sense that at similar peak stress of 31 MPa total failure occurs at a higher external strain. The local model shows some

earlier in Section 3. Fig. 14 plots the stress-strain curves of a cubic specimen with material properties as described in Section 3.1. In this plot, the maximum element size is varied. We use the maximum element sizes as 1, 2 and 4 mm whereas the total number of elements are 40.1, 18.9 and 9.0

12

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 16. The damage evolution of L-shaped specimen for different mesh with the local model; for a mesh with a maximal (a) element size 8 mm, (b) element size 4 mm, (c) element size 2 mm and (d) element size 1 mm. (e) shows a zoom in picture of (c); (f) shows a zoom in picture of (d). The color represents the value of the damage parameter D.

weak mesh dependence in particular after loading has passed the peak stress whereas for the non-local model the stress-strain curves for the various mesh sizes are identical indicating the mesh independence of the model. For the L-shaped specimen introduced in Section 3.2 mesh dependence of the local model is more apparent. Fig. 15 shows the stressstrain curves for different mesh element sizes. We use the maximum element sizes as 1, 2 and 4 mm with the total numbers of elements of 83.9, 12.5 and 1.7 million respectively. In addition, the result from the local model with a maximum mesh element size of 8 mm identical to the localization length of the non-local model is also presented in Fig. 15 (a). The figure demonstrates strong mesh dependence of the

stress-strain curves. On the other hand the results from the non-local model as shown in Fig. 15 (b) display stress-strain curves independent from the mesh size. Similar to the cubic specimen the non-local approach leads to material which is stronger but also more brittle (due to the large localization length employed) indicated by the sharp drop of stress after peak is reached. It is highlighted that the material behavior described by the non-local model with a localization length of l = 8 mm is not identical to the local model simulation with a mesh size of 8 mm, as indicated by the peak stresses shown in Table 5. It is also pointed out that there is a relatively large jump in the peak stress for the non-local model when reducing the element size from the 8 mm to 4 mm, see Table 5 (stress-strain curve for element size 8 mm is not shown). This

13

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

Fig. 17. The damage evolution of L-shaped specimen for different meshes in non-local model with constant localization length l = 8 mm; for a mesh with a maximal (a) element size of 8 mm, (b) element size of 4 mm, (c) element size 2 mm and (d) element size 1 mm. The color represents the value of the damage parameter D.

observation provides the justification for the recommendation that the non-local model requires a mesh with an element size half the localization length in order to be able to deliver acceptable results. The impact of element size on the damage pattern is illustrated by Fig. 16 for the local model and by Fig. 17 for the non-local model where the color represents the value of the damage parameter D. Again results are shown for different meshes with various maximum element sizes and for a fixed localization length l = 8 mm. In Fig. 17 the damage plots are shown immediately after the peak stress has reached. The patterns for the non-local model are identical where for large element sizes the zone of large damaged values D is more lacerated for larger element sizes due to the resolution of the mesh at the boundary of the damaged zone. As a comparison the damage patterns for different mesh sizes produced by the local model are shown in Fig. 16. From this figure, it is clear that the width of the damaged zone is narrower than for the nonlocal model, where the width is becoming narrower as the mesh is refined. To better highlight the effect of mesh dependence for the two very fine meshes of element size 2 and 1 mm the zone where damage localizes has been enlarged, see Fig. 16 (e) and (f) respectively. For these two cases the width of the damaged zone in the center of the arc formed by damaged zone has similar size but the thickness of the arc is thinning out towards reentrant corner and the opposite boundary. The final thickness at both ends is given by the element size of the mesh which gives evidence for the mesh dependence of the damage pattern produced by the local model which is consistent with the mesh dependency in stress-strain curves of Fig. 15 (a). Finally the mesh dependence of the simulation results for the

rectangular specimen with a pre-existing flaw as introduced in Section 3.3 is investigated. Fig. 18 shows the stress-strain curves for meshes with the element size of 1, 2 and 4 mm. The respective mesh sizes are 55.3, 9.40 and 2.0 million elements. The stress-strain curves are independent of element size for the non-local model in Fig. 18 (b). Again the results of the local model show mesh dependence which applies to the stress-strain curve as well as to the peak stress, see Fig. 18 (a). Analogously to the previous cases a reduction in the maximum element size results in reduction of strength of the material. However for this case the simulated compression shows no significant change to overall brittleness of the specimen with element refinement as well as in comparison to the non-local model. Table 5 summarizes the peak strength obtained for the three test specimens for the local and non-local model for different mesh resolutions. For the cubic mesh with a weak zone, peak stresses show similar values between the local and non-local model where the nonlocal model appears to produce slightly less brittle behaviors, see Fig. 14. For the L-shape specimen the local model shows strong mesh dependence in the peak stress and exhibits a lower degree of brittleness for smaller element sizes and in comparison to behavior shown by the non-local model, see Fig. 15. A similar trend is seen in the block specimen with a pre-existing flaw, see Fig. 18. One explanation for the Lshape specimen exhibiting the strongest mesh dependence is that the mesh involves no complex geometries such as the flaw or weak zone, and the mesh elements are more or less uniform throughout the specimen. However for the flawed specimens the mesh is usually more refined closer to the flaw boundaries. 14

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al.

examples presented in this paper, we found that even smaller localization length scales were needed to realistically reproduce the narrow fracture widths observed in the experiments. As a consequence running 3-D simulations producing mesh independent results on a physically meaningful spatial scale remains difficult to achieve even using a massively parallelized code. CRediT authorship contribution statement S. Mondal: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Validation, Visualization, Writing - original draft. L.M. Olsen-Kettle: Conceptualization, Funding acquisition, Software, Investigation, Supervision, Writing - review & editing. L. Gross: Conceptualization, Resources, Software, Investigation, Supervision, Writing - review & editing. 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. Acknowledgements This research is supported by the Australian Research Council Discovery Early Career Researcher Award DE140101398. The first author is grateful for PhD scholarship support, the tuition fee award and research assistance from this award, The University of Queensland and the Swinburne University of Technology, Australia respectively. The authors also gratefully acknowledge the software and high performance computing (HPC) support from The University of Queensland. This work uses infrastructure funded by the AuScope National Collaborative Research Infrastructure Strategy of the Australian Commonwealth. This project was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government.

Fig. 18. Block specimen with a pre-existing flaw: The stress-strain curve for varying mesh resolutions for (a) the local model and (b) the non-local model with localization length of l = 8 mm.

References

6. Conclusions

[1] Bažant ZP, Ožbolt J. Nonlocal microplane model for fracture, damage, and size effect in structures. J Eng Mech 1990;116(11):2485–505. [2] Bažant ZP, Pijaudier-Cabot G. Nonlocal continuum damage, localization instability and convergence. J Appl Mech Trans ASME 1988;55:287. [3] Belytschko T, Lasry D. A study of localization limiters for strain-softening in statics and dynamics. Comput Struct 1989;33(3):707–15. [4] Borst RD, Sluys L, Mühlhaus H, Pamin J. Fundamental issues in finite element analyses of localization of deformation. Eng Comput 1993;10(2):99–121. [5] Borst R, Pamin J, Peerlings R, Sluys L. On gradient-enhanced damage and plasticity models for failure in quasi-brittle and frictional materials. Comput Mech 1995;17(1):130–41. [6] Caner FC, Bažant ZP. Microplane model m7 for plain concrete. ii: Calibration and verification. J Eng Mech 2013;139(12):1724–35. [7] Cazes F, Meschke G, Zhou M-M. Strong discontinuity approaches: An algorithm for robust performance and comparative assessment of accuracy. Int J Solids Struct 2016;96:355–79. [8] Claudia C. Computational modelling of gradient-enhanced damage in quasi-brittle materials. Mech Cohesive-frict Mater 1999;4(1):17–36. [9] Codd A, Gross L. Electrical Resistivity Tomography using a finite element based BFGS algorithm with algebraic multigrid preconditioning. Geophys J Int 2017;212(3):2073–87. [10] de Vree J, Brekelmans W, van Gils M. Comparison of nonlocal approaches in continuum damage mechanics. Comput Struct 1995;55(4):581–8. [11] del Viso JR, Carmona JR, Ruiz G. Shape and size effects on the compressive strength of high-strength concrete. Cem Concr Res 2008;38(3):386–95. [12] Ellery A, Fenwick J, Gross L. esys-escript 5.3; 2019. https://espace.library.uq.edu. au/view/UQ:08a47ea. [13] Engelen RA, Geers MGD, Baaijens FP. Nonlocal implicit gradient-enhanced elastoplasticity for the modelling of softening behaviour. Int J Plast 2003;19(4):403–33. [14] Geers M. Experimental analysis and computational modelling of damage and fracture [PhD thesis]. Department of Mechanical Engineering; 1997. [15] Geers MGD, Ubachs RLJM, Engelen RAB. Strongly non-local gradient-enhanced finite strain elastoplasticity. Int J Numer Meth Eng 2003;56(14):2039–68.

Using the implicit gradient damage formulation a local damage model is modified with the objective to eliminate mesh dependence. The idea is to introduce regularization of the equivalent strain to eliminate mesh dependence and strain localization for the strain-softening damage model. This non-local version is constructed by solving an additional PDE – the Helmholtz equation. In contrast to the local damage model based on a dual damage initiation condition as investigated previously by Mondal et al. [31], we employ an equivalent strain condition is used for the damage evolution as it requires the regularization of a single quantity only (the equivalent strain) which is computationally less costly. The momentum and the regularization equations which are formally coupled via the damage parameter are solved sequentially using a split-operator approach. To validate the model three specimens were investigated. The simulation results show that local and non-local model produce similar damage patterns but in the non-local model the localization length introduces an additional length scale in the regularization. It determines the extent of the damaged zone and subsequent fracture propagation and fracture widths. As confirmed by the simulations smaller localization lengths lead to a narrower fracture width. The presented numerical investigations allow to conclude that the implicit gradient damage formulation in fact removes the mesh dependence from the local damage model. However, there are computational limitations for the application of the non-local model for practical cases as it is extremely computationally expensive to run simulations with smaller localization length scales. For the 15

Computers and Geotechnics 122 (2020) 103505

S. Mondal, et al. [16] Geuzaine C, Remacle J-F. Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities. Int J Numer Meth Eng 2009;79(11):1309–31. [17] Gross L, Altinay C, Shaw S. Inversion of potential field data using the finite element method on parallel computers. Comput Geosci 2015;84:61–71. [18] Gross L, Bourgouin L, Hale AJ, Muhlhaus H-B. Interface modeling in incompressible media using level sets in escript. Phys Earth Planet Inter 2007;163:23–34. [19] Jirasek M. Nonlocal models for damage and fracture: Comparison of approaches. Int J Solids Struct 1998;35(31):4133–45. [20] Kröner E. Elasticity theory of materials with long range cohesive forces. Int J Solids Struct 1967;3(5):731–42. [21] Kuhl E, Ramm E, de Borst R. An anisotropic gradient damage model for quasi-brittle materials. Comput Methods Appl Mech Eng 2000;183(1):87–103. [22] Lasry D, Belytschko T. Localization limiters in transient problems. Int J Solids Struct 1988;24(6):581–97. [23] Liu Y, Murakami S, Hayakawa K. Mesh dependence and stress singularity in local approach to creep-crack growth analysis. Trans Jpn Soc Mech Eng Ser A 1993;59(564):1811–8. [24] Lu Y, Wang L, Elsworth D. Uniaxial strength and failure in sandstone containing a pre-existing 3-d surface flaw. Int J Fract 2015;194(1):59–79. [25] Luzio GD. A symmetric over-nonlocal microplane model m4 for fracture in concrete. Int J Solids Struct 2007;44(13):4418–41. [26] Mediavilla J, Peerlings R, Geers M. An integrated continuous–discontinuous approach towards damage engineering in sheet metal forming processes. Eng Fract Mech 2006;73(7):895–916. [27] Mühlhaus H-B, Alfantis E. A variational principle for gradient plasticity. Int J Solids Struct 1991;28(7):845–57. [28] Mühlhaus H, Gross L. A breakage diffusion model for strength softening rock. Bifurcation and degradation of geomaterials with engineering applications. 2017. p. 563–9. [29] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput Methods Appl Mech Eng 2010;199(45):2765–78. [30] Mondal S. Numerical modelling of damage evolution and fracture propagation in brittle materials [PhD thesis]. School of Earth and Environmental Sciences, The University of Queensland; 2019. [31] Mondal S, Olsen-Kettle L, Gross L. Simulating damage evolution and fracture propagation in sandstone containing a preexisting 3-d surface flaw under uniaxial compression. Int J Numer Anal Meth Geomech 2019;43(7):1448–66. [32] Murakami S, Liu Y. Mesh-dependence in local approach to creep fracture. Int J Damage Mech 1995;4(3):230–50. [33] Nguyen TH, Niiranen J. A second strain gradient damage model with a numerical implementation for quasi-brittle materials with micro-architectures. Math Mech Solids 2020;25(3):515–46. [34] Ožbolt J, Bažant ZP. Numerical smeared fracture analysis: Nonlocal microcrack interaction approach. Int J Numer Meth Eng 1996;39(4):635–61. [35] Oliver J, Huespe A, Blanco S, Linero D. Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach. Comput Methods Appl Mech Eng 2006;195(52):7093–114. [36] Olsen-Kettle L. Quantifying the orthotropic damage tensor for composites undergoing damage-induced anisotropy using ultrasonic investigations. Compos Struct

2018;204:701–11. [37] Olsen-Kettle L. Using ultrasonic investigations to develop anisotropic damage models for initially transverse isotropic materials undergoing damage to remain transverse isotropic. Int J Solids Struct 2018;138:155–65. [38] Olsen-Kettle L. Bridging the macro to mesoscale: Evaluating the fourth-order anisotropic damage tensor parameters from ultrasonic measurements of an isotropic solid under triaxial stress loading. Int J Damage Mech 2019;29:219–32. [39] Paknahad A, Kucko NW, Leeuwenburgh SC, Sluys LJ. Experimental and numerical analysis on bending and tensile failure behavior of calcium phosphate cements. J Mech Behav Biomed Mater 2020;103:103565. [40] Peerlings R. Enhanced damage modelling for fracture and fatigue [PhD thesis]. Department of Mechanical Engineering, Technische Universiteit Eindhoven; 1999. [41] Peerlings R, Borst RD, Brekelmans WAM, Vree JHPD. Gradient enhanced damage for quasi-brittle materials. Int J Numer Meth Eng 1996;39(19):3391–403. [42] Peerlings R, de Borst R, Brekelmans W, Geers M. Localisation issues in local and nonlocal continuum approaches to fracture. Eur J Mech A Solids 2002;21(2):175–89. [43] Peerlings R, Geers M, de Borst R, Brekelmans W. A critical comparison of nonlocal and gradient-enhanced softening continua. Int J Solids Struct 2001;38(44):7723–46. [44] Peerlings RHJ, Brekelmans WAM, De Borst R, Geers MGD. Gradient-enhanced damage modelling of high-cycle fatigue. Int J Numer Meth Eng 2000;49(12):1547–69. [45] Peerlings RHJ, de Borst R, Brekelmans WAM, Geers MGD. Gradient-enhanced damage modelling of concrete fracture. Mech Cohesive-frict Mater 1998;3(4):323–42. [46] Pires-Domingues SM, Costa-Mattos H, Rochinha FA. Nonlinear elastic brittle damage: numerical solution by means of operator split methods. J Theor Appl Mech 1999;37(4):847–61. [47] Schaa R, Gross L, Plessis Jd. Pde-based geophysical modelling using finite elements: examples from 3d resistivity and 2d magnetotellurics. J Geophys Eng 2016;13(2):S59–73. [48] Titscher T, Oliver J, Unger JF. Implicit-explicit integration of gradient-enhanced damage models. J Eng Mech 2019;145(7):04019040. [49] Wang SY, Sloan SW, Sheng DC, Yang SQ, Tang CA. Numerical study of failure behaviour of pre-cracked rock specimens under conventional triaxial compression. Int J Solids Struct 2014;51(5):1132–48. [50] Winkler B, Hofstetter G, Niederwanger G. Experimental verification of a constitutive model for concrete cracking. Proc Inst Mech Eng Part L: J Mater Des Appl 2001;215(2):75–86. [51] Wong T-F, Wong RH, Chau K, Tang C. Microcrack statistics, Weibull distribution and micromechanical modeling of compressive failure in rock. Mech Mater 2006;38(7):664–81. [52] Zhu CW. Numerical modelling of the effect of rock heterogeneity on dynamic tensile strength. Rock Mech Rock Eng 2008;41(5):771–9. [53] Zhu CW, Tang AC. Micromechanical model for simulating the fracture process of rock. Rock Mech Rock Eng 2004;37(1):25–56. [54] Zienkiewicz OC, Taylor RL, Zhu JZ. The finite element method: its basis and fundamentals. 7th ed. Elsevier; 2013. [55] Zreid I, Kaliske M. Regularization of microplane damage models using an implicit gradient enhancement. Int J Solids Struct 2014;51(19):3480–9.

16