Analysis of a model of field crack mechanics for brittle materials

A computational model for arbitrary brittle crack propagation, in a fault-like layer within a 3-d elastic domain, and its associated quasi-static and dynamic fields is developed and analyzed. It uses an FFT-based solver for the balance of linear momentum and a Godunov-type projectionevolution method for the crack evolution equation. As applications, we explore the questions of equilibria and irreversibility for crack propagation with and without surface energy, existence of strength and toughness criteria, crack propagation under quasi-static and dynamic conditions, including Modes I, II and III, as well as multiaxial compressive loadings.


Introduction
This paper develops a computational strategy and uses it to critically evaluate the capabilities of a recently proposed continuum mechanical model of fracture that we shall refer to as Field Crack Mechanics (FCM) (Acharya, 2018(Acharya, , 2020, restricted here to brittle bulk response and to evolution of arbitrary crack patterns in a fixed plane. Beyond basic verification, the evaluation is carried out in the context of analyzing several illustrative initial-boundary-value problems in nonlinear fracture involving elasto-statics/dynamics in three space dimensions. The model has been placed in the context of existing literature in the aforementioned references, in particular classical fracture mechanics (e.g. the monograph by Freund (1998), Salvadori and Fantoni (2016)), cohesive zone approaches to fracture (e.g. Barenblatt (1962); Dugdale (1960); Needleman (1990)), the eigenerosion approach to fracture (e.g. Pandolfi and Ortiz (2012)), phase field methods for fracture (e.g. Bourdin et al. (2008Bourdin et al. ( , 2000; Tanné et al. (2018); Staroselsky et al. (2019)) regularizing the variational approach to fracture (Francfort and Marigo (1998); Dal Maso and Toader (2002)) and also independently conceived from the Ginzburg-Landau paradigm (e.g. Hakim and Karma (2009)), peridynamics (e.g. Silling and Lehoucq (2010); Lipton (2014)), and dislocation based approaches to fracture (e.g. Weertman (1996); Bilby et al. (1963); Fressengeas and Taupin (2014)); of note is also the work of Smyshlyaev and Willis (1994) on accounting for crack-face contact constrains in dynamic fracture, and the extension of the variational approach to fracture to elastodynamics (Dal Maso et al. (2016)).
Our computational strategy adapts the prior work of Morin et al. (2019) for the computation of the Field Dislocation Mechanics (FDM) model of plasticity to FCM, with extension to the entire range of crack evolution under quasi-static to dynamic balance of forces. Notably, to our knowledge, our work involves the first fully dynamic implementation of the FFT method for solid mechanics due to Moulinec and Suquet (1998). In this connection, we develop a general set of sufficient conditions that guarantee formal uniqueness of solutions of the linear, spatio-temporally heterogeneous, 'periodic' elastodynamic system relevant to the class of problems we analyze. The FFT method for 'stress analysis' is combined with a central-upwind scheme for a Hamilton-Jacobi equation representing the evolution of the crack field.
Two distinguishing features of our physical model are the representation of crack patterns by a vector field whose direction at any field point represents the direction of the crack-surface normal, with its magnitude representing the extent of the elastic modulus degradation, i.e. brittle damage, at that location. This uniquely equips the model to distinguish between compressional and extensional loadings normal to the crack surface in a physically transparent, direct manner, without reliance on volumetric-deviatoric (Amor et al. (2009)) or eigen-decomposition (Ortiz (1985); Miehe et al. (2010); Borden et al. (2012)) splits of the strain tensor that do not work uniformly (consider the possibility of deviatoric-strain energy driven fracture propagation of a planar crack under transverse compression utilizing the volumetric split and see, e.g., Agrawal (2016), for a simple example related to simple shearing vis-a-vis the eigen-decomposition split); an attempt to define a crack-surface normal field, and the difficulties in doing so in a physically appropriate manner, within a 'scalar' phase field model of fracture can be appreciated from Agrawal (2016). The evolution of our vectorial crack field is a direct consequence of a conservation law for topological charge associated with the crack-tip field, constrained by the second law of thermodynamics. Crack propagation in the presence of material inertia is a natural consequence of the model, as is the representation of zero-surface energy cracks leading to irreversibility in crack evolution under unloading, without relying on any irreversibility criteria put in 'by hand.' This latter fact is the second distinguishing feature of our model: allowing for the possibility of crack growth without any surface energy cost, while maintaining classical expectations on prediction of strength and toughness. Furthermore, while lending itself naturally to non-singular crack-tip stress fields and geometric representation of cracks by a field in contrast to classical fracture mechanics, it is a manifest feature of the theory that crack evolution is localized around the crack-tip field, regardless of levels of stress or energy, or elastic wave propagation in the body, much like physical notions of classical fracture. This paper is organized as follows. In Section 2, the Field Crack Mechanics theory is recalled and a planar ansatz, in which the crack field is constrained in a thin layer, is derived. A numerical formulation of the model is proposed in Section 3, relying on a FFT-based solver for the balance of linear momentum and Godunov-type projectionevolution solver for the crack evolution equation. In Section 4, several problems including Mode I, shear and compressive loadings are investigated in quasi-static conditions. Crack evolution in dynamic loadings is then considered in Section 5. Section 6 demonstrates the capability of our framework in modeling progressive fracture in materials with a random initial (small) crack microstructure. Section 7 contains concluding remarks, and Appendix A contains arguments for formally deducing uniqueness of solutions in spatio-temporally heterogeneous, periodic, elastodynamics.
2 The Field Crack Mechanics model of brittle fracture

General case
The problem we are addressing is the numerical simulation of brittle crack propagation in an elastic solid. The approach that we follow relies on the use of a crack field as an internal variable (Acharya, 2018(Acharya, , 2020. This requires solving (i) for mechanical force balance (with or without inertia) for a given crack field and some applied macroscopic loading in an heterogeneous anisotropic elastic media and (ii) an evolution problem, consisting of the evolution of the crack field involving the local mechanical and crack-related energetic fields constrained by kinematics related to conservation of topological charge of the crack-tip field.
The objective of this work is to study the local and overall mechanical fields of brittle materials with complex microstructures containing cracks that can propagate. Thus, this work focuses on the effect of cracks within some heterogeneous material and is in line with previous studies devoted to the simulation of materials with heterogeneous microstructures for which it is well established that FFT-based methods are adapted to solve the mechanical problem (Moulinec and Suquet, 1998). Indeed, it is a powerful tool to solve the elasticity equations with heterogeneous fields that are periodic, such as in polycrystalline or composite materials. In the case of cracked media, heterogeneity is also built-in as the crack field modifies the local elastic properties, which makes this method appealing. In the case of a phase-field fracture model, the FFT-based approach has proved to provide valuable insights in brittle crack propagation problems (Chen et al., 2019;Ernesti et al., 2020;Ma and Sun, 2020). Periodic boundary conditions, intrinsic of the method 1 , are of course slightly restrictive in crack propagation problems, but the benefits of the method outweigh this limitation in several practical cases as, e.g. the solution to arbitrary static crack distributions (Gasnier et al., 2018), problems with random distributions of cracks, and in general any large scale problems with complex microstructures.
With in mind the use of FFT-based methods, we are thus concerned with a medium with a periodic microstructure so that the crack field vector c(x, t) and the (possibly damaged) elastic moduli tensor C(x, c(x, t)) are periodic fields. The unit-cell is a 3D cubic domain Ω = [−L, L] 3 with 2L the period of the microstructure. Tensorial components refer to a system of Cartesian coordinates (e 1 ; e 2 ; e 3 ). In quasi-statics, the problem consists in finding, the local stress, strain and displacement fields (respectively denoted by σ(x, t), ε(x, t) and u(x, t)) and the crack field c(x, t), for given initial conditions on c and some prescribed, generally time-dependent boundary conditions. In dynamics, initial conditions on the displacement u and velocityu are required, in addition. These fields satisfy (in the absence of body force) where ρ is the time-independent mass density field (corresponding to the reference configuration of the body from which all displacements are measured) and V is the crack velocity field. In the above, the local strain ε(x) is split into its average ε = ε (where · denotes the spatial average over the unit-cell Ω) and a fluctuation term ε(u * (x)): The fluctuating part u * (x) is periodic (notation: u * (x)#) and the traction σ.n is antiperiodic on the boundary between two neighboring cells with n the outer normal along the boundary ∂Ω of Ω (notation: σ.n − #).
It should be noted that the crack evolution equation can be modified in order to account for some irreversibility condition of crack propagation in order to prevent crack healing; this condition would reaḋ The free energy function ψ is supposed to be of the form where ψ E is the elastic strain energy density of the material (with its elastic modulus degraded) and ϕ is an energy function which should account notably for an energy barrier to damage from an undamaged state. The positivity of the mechanical dissipation D implies that the crack velocity is of the form (Acharya, 2018(Acharya, , 2020) where M is a symmetric, positive definite tensor of crack mobility; in the isotropic case it is given by where I is the identity matrix and B is a (positive) scalar drag coefficient.
We assume that the elastic strain energy density is of the form In this equation, C i denotes the elasticity tensor of the intact material and C d corresponds to the damaged moduli (depending on the crack field). The function H(x) is the Heaviside step function which is given by The last term in equation (7) permits to introduce a different behavior in the compressive regime in the direction normal to the crack surface in order to prevent the interpenetration of crack flanks (corresponding to an undamaged behavior or some prescribed contact stiffness), which corresponds to the case ε n < 0 where ε n is given by In the case of a non-damaged behavior in compression, the constant α would read From the definition of the elastic strain energy, the stress tensor is given by where δ is the Dirac delta function. Assuming the term δ(ε n )ε n to be negligible (since, in the sense of generalized functions (cf. Bracewell (1986)), the product of δ(x) with x is equal to 0), the stress tensor reduces to This allows the expression of the elasticity tensor: Typical examples of intact and damaged elastic behaviors C i and C d in the isotropic case (in the strain argument) are given by where λ and µ are Lamé's coefficients of the intact material,λ(|c|) andμ(|c|) are Lamé's coefficients of the damaged material, I is the fourth-order identity tensor and I is the second-order identity tensor. Simple forms for the damaged elastic moduliλ andμ are monotonically decreasing convex functions of the magnitude of cracking: where f is a degradation function of a non-negative scalar argument of the form where f m 1 is the coefficient describing the maximal degradation of the elastic modulus and d c is a threshold value of |c| for inducing this maximal degradation in elastic stiffness.
The simplest form of the crack energy density ϕ may be expected to be of the form where η(|c|) is a non-convex local crack energy density function representing an energy barrier to crack growth and the other term is the crack-tip energy density representing a lowest-integer-order approximation of any smooth function that assigns an energy cost to the formation of a crack-tip. In this paper we will deal with two candidate functions sketched in Figure 1. The functions correspond to qualitatively different mechanisms of energy storage in cracked regions. The function η 1 represents a 'Griffith' type energy cost to cracked regions, endowing a permanent energy cost to the formation of a crack. The function η 2 , on the other hand, while presenting an energy barrier that has to be surmounted for a material point to be cracked leaves no residual stored energy (beyond a small amount of strain energy density, regulated by the residual stiffness related to the coefficient f m ) once the material is at or beyond the damage threshold ω 2 . Before proceeding further, we nondimensionalize our model.
In the presence of quasi-static loading (i.e., rates of loading asymptotically slower than the time scales of elastic wave propagation in the body or the relaxation dictated by the drag B) the dimensional equations to be solved are div(σ) = 0 and −curl(c) × V = 0 (considerations of the irreversibility condition do not change the essential argument here).
Non-dimensionalizing physical quantities with dimensions of length by L a linear extent of the domain, stresses and energy densities by the shear modulus µ (we consider c as dimensionless, cf. Acharya (2020)), assuming the crack mobility tensor to be of the form (6) and denoting all dimensionless quantities by a 'circle superscript', the resulting nondimensional system of equations is The rates of evolution are entirely controlled by the time-rates of the applied boundary conditions.
For rates of loading (perhaps 'slow') where inertial effects cannot be ignored, we again nondimensionalize lengths by L, stresses and energy densities by µ, and time by L vs , where v s = µ ρ is the elastic shear wave speed, to obtain the following non-dimensional system: (20) In both cases the irreversibility condition is appended when used.
In the rest of the paper we deal with the above nondimensional systems and, for notational convenience, all 'circle exponents' will be dropped, except when explicitly noted to the contrary.

A planar ansatz
Hereafter we consider a simplified planar ansatz in which the crack field is confined in a thin layer (see Figure 2). This constitutes a first step toward the simulation of fully 3d crack propagation in heterogeneous materials and follows developments in Morin et al. (2019) and Zhang et al. (2015) for the case of dislocation mechanics.
The two outer regions outside of the layer are supposed to be purely elastic isotropic linear with material constants µ and λ. This problem will allow to reduce the crack evolution equation to a scalar degenerate parabolic equation with dominant nonlinear wave transport behavior for which efficient numerical solvers can be deployed.
In this planar problem, the following ansatz is assumed: (1) The crack surface is supposed to be the plane x 3 = 0. The crack field c and the crack-tip velocity field V are supposed to be constrained in a layer of finite thickness h around the crack surface, and to be of the form Non-damaged elastic regions

Planar layer
Crack field x 1 x 3 x 2 (2) The elastic energy (7) is given by 2ϕ E (ε, c 3 ) = C 1111 (ε 2 11 +ε 2 22 )+C 3333 ε 2 33 +2C 1122 (ε 11 ε 22 +ε 11 ε 33 +ε 22 ε 33 )+2C 1212 (ε 2 12 +ε 2 13 +ε 2 23 ), (21) or alternatively where I 1 and I 2 are functions of some symmetric second-order tensor A given by The non-negative components of the stiffness tensor C, accounting for the tensioncompression asymmetry according to the crack normal e 3 , are given by The elasticity behavior is thus heterogeneous due to the crack field c 3 , non-linear and tension-compression asymmetric with respect to ε 33 due to the compressive behavior of the crack flanks. In a case of a non-damaged compressive behavior of the crack flanks, the constant α reads In the following we will only consider this case.
(3) The degradation function appearing in the damaged moduli (15) reduces to (4) The evolution equation for the crack field consistent with (24) and this ansatz iṡ where · layer represents an average of · over the height of the layer. Since the crack is defined only in the planar layer, the driving force related to the strain energy density is averaged within the layer (this is motivated from a similar result in the FDM model (Zhang et al. (2015)) where it is shown that such a fault-thickness averaged elastic driving 'force' is adapted to the assumed layer-kinematics ansatz, and exactly satisfies thermodynamics). In the case of a non-damaged compressive behavior where α is given by equation (25), one has where h is the thickness of the layer.

General procedure
The resolution will consist in finding the mechanical state S n+1 = {ε n+1 , u n+1 , σ n+1 , c n+1 } at time t n+1 , knowing the previous mechanical state S n = {ε n , u n , σ n , c n } at time t n and considering boundary conditions (typically prescribed overall strain ε = ε or stress σ = σ ). We follow the strategy adopted in Morin et al. (2019) for FDM problems, which consists in treating separately the elasticity problem and the evolution problem through a staggered strategy: • The crack evolution equation will be solved using Godunov-type solvers; • Two distinct solvers will be developed for the elasticity problem, a non-linear FFTbased solver in the absence of inertial effects (elastostatic problem) and a dynamic explicit solver in the presence of inertial effects (elastodynamic problem).

Elastostatic problem
The static problem will be solved using an FFT-based solver. This class of solvers is particularly interesting in heterogeneous elasticity problems (Moulinec and Suquet, 1998), as it is the case with cracks (Gasnier et al., 2018) (see also (Rovinelli et al., 2020)). The present case differs slightly from that considered by Moulinec and Suquet (1998) and Gasnier et al. (2018) because the elastic behavior is non-linear due to the condition of non-interpenetration of crack flanks (the components of the elasticity tensor depend on ε).
Assuming that the present crack field c n+1 is known, the static problem is first written as the so-called auxiliary problem (Moulinec and Suquet, 1998) where C 0 is the stiffness tensor of some homogeneous medium and τ (x) is the polarization tensor given by In practice, the reference stiffness tensor C 0 can be taken isotropic with Lamé's coefficients µ 0 and λ 0 generally given by (Moulinec and Silva, 2014) The principle of FFT-based solvers is to iterate in order to determine the correct polarization tensor solution of the problem. The basic scheme of Moulinec and Suquet (1998) permits to solve linear (elastic) as well as non-linear (visco-plastic) problems; here this scheme is slightly adapted in order to account for the non-linearity of the stiffness tensor with respect to the strain field. In fact it is actually straightforward to consider this non-linearity since, with this method, the stress field is deduced from the strain field by the constitutive law. The following fixed point algorithm is thus considered in order to determine τ (x): In step (d), Γ 0 is the classical (continuous) Green operator related to the reference medium C 0 (Moulinec and Suquet, 1998). Convergence is reached when σ i n+1 is in equilibrium. In practice, the iterative procedure is stopped when the error serving to check convergence (Moulinec and Suquet, 1998) is smaller than a prescribed (small) value (typically 10 −4 in practice).
It should be noted that other derivatives rules (such as the finite differences rule proposed by Willot et al. (2014)) or the rotated finite difference scheme of Willot (2015) could be used to improve the numerical results. In particular, the rotated finite difference scheme of Willot (2015) would permit to handle infinite contrasts using the fixed point method as shown by Djaka et al. (2017). Nevertheless, we have checked that the continuous operator together with the fixed point were sufficient to handle cracks with reasonable (finite) values for the degradation function (roughly f m = 0.01).

Elastodynamic problem
We develop a fully dynamic, FFT based explicit solver for the balance of linear momentum including the effects of material inertia. In Appendix A we deduce sufficient conditions for formal uniqueness of solutions in spatio-temporally heterogeneous and periodic elastodynamics that indicate an appropriate set of supplementary conditions that allow for a well-set problem in this 'periodic' setting. The local acceleration field is assumed to be of the formü where the fluctuation of the accelerationü * (x) is a periodic field. The macroscopic acceleration term,ε, is assumed to vanish in the algorithm. This is actually not restrictive because piecewise constant strain rate loadings are admissible within this setting, as argued in Appendix A.
The following algorithm (explicit central difference in time) is thus considered: Step time n + 1 u * n , u * n , ü * n , σ n and c n+1 being known In practice, the dynamic explicit solver was found to be quite sensitive to numerical artifacts (oscillations), emerging from pseudospectral differentiation in the presence of traction jumps (Morin et al., 2021). Thus, to improve the computation of the local field and remove spurious oscillations, the pseudospectral differentiation rule is replaced by finite differences. Following Willot et al. (2014) (see also Berbenni et al. (2014)), the strain field is computed from the displacement field using a backward finite difference rule; in step (b) the term ξ i is replaced by k i given by The acceleration is, on the other hand, computed from the stress field using a forward finite difference rule; in step (f) the term ξ i is replaced by k i given by We have verified our implementation by the 'method of manufactured solutions'.

Crack evolution problem
Following Morin et al. (2019), we consider the central scheme of Kurganov et al. (2001) for the resolution of the crack evolution problem, because it is a nonoscillatory high-resolution Godunov-type projection-evolution method. The crack evolution problem given by (27) can be formally written as where the energetic driving force, e 0 , is given by This term e 0 will be evaluated at the beginning of the time step; following Kurganov et al. (2001), the Laplacian (diffusion) term in equation (38) is discretized by a fourthorder central difference scheme.
A uniform grid (defined at the voxels' nodes) is chosen with the following notations: x j = j∆x, y k = k∆y, t n = n∆t and c n jk = c(x j , y k , t n ), where ∆x, ∆y and ∆t are respectively the spatial step sizes and the time step. Provided that the value of c is known at time t n , we are looking for the point values of c at time t n+1 . The algorithm is briefly detailed below: Step 1: construction of a continuous piecewise interpolant. A continuous piecewise quadratic interpolantc(x, y, t n ) is constructed in order to avoid spurious oscillations (see Kurganov and Tadmor (2000) for the full detail). A nonlinear limiter is used in the definition of the quadratic interpolant in order ensure the nonoscillatory behavior of the central scheme. In practice, a one-parameter family of the minmod limiters is considered (Morin et al., 2019).
Step 2: estimation of the one-sided local speed of propagation. The one-sided local speeds of propagation in the x− and y−directions are evaluated at the grid point (x j , y k ): where c ± x =c x (x j ± 0, y k , t n ) and c ± y =c y (x j , y k ± 0, t n ) are the right and the left derivatives in the x− and y−direction, deduced from the quadratic interpolant (Kurganov and Tadmor, 2000).
Step 3: approximate solution of the Hamilton-Jacobi equation at intermediate grid points. The Hamilton-Jacobi equation (36) is then solved at the intermediate points (x n j± = x j + a ± jk ∆t, y n k± = y k + b ± jk ∆t). This leads to the approximate Riemann solver c n+1 j±,k± =c(x n j± , y n k± , t n ) − ∆t H c x (x n j± , y n k± , t n ),c y (x n j± , y n k± , t n ) .
Step 4: projection of the intermediate solution onto the original grid. The solution previously obtained at the intermediate points x j± and y k± is projected onto the original grid which leads to the fully discrete scheme In problems of dislocation motion with constant velocity, the predictions of this solver almost coincide with the exact solution (Morin et al., 2019).
4 Crack evolution in quasi-static loadings

Description of the simulations
Quasi-static loadings can be performed by combining the FFT-based solver for the static problem with the Godunov-type algorithm for the evolution problem. In practice, we apply successive macroscopic strain increments ∆ε. In order to ensure the quasi-static evolution, we let the crack field evolve for a given macroscopic strain increment. When the crack field is converged, we move to the next macroscopic strain increment and so on.
The algorithm related to the quasi-static runs is provided in Figure 3.
Evolution problem (Central upwind scheme) Knowing the previous strain and crack fields: Update of macroscopic loading ε n+1 = ε n + ∆ε No update of macroscopic loading ε n+1 = ε n
Since FFT algorithms require data sampled on a regular grid, we consider throughout the paper a [−1, 1] × [−1, 1] × [−1, 1] regular unit-cell discretized with 128 × 128 × 128 voxels (the mechanical fields are evaluated at the center of each voxel) so the spatial step is ∆x = 0.015625 and we have L = 1. The planar layer is supposed to be centered in the plane x 3 = 0 and a constant thickness h = 0.046875 is considered for the layer, which corresponds to a number of 3 voxels 2 . The following values are taken for the elastic properties: µ = 1 and λ = 2 (which corresponds to a Poisson ratio ν = 1/3, typical of metallic materials). Finally the values e c = 10 −5 and e i = 10 −4 have been chosen for the convergence criteria, where e i serves in the convergence criterion of the iterative FFTbased solver (see equation (32)) while e c serves in the convergence test of the evolution problem (see Figure 3). The (artificial) time step ∆t which enters in the evolution equation is chosen in order to check the following stability constraint In equation (42), v l is the maximal local speed of propagation given by where a + ij ,a − ij , b + ij and b − ij are the local speeds of propagation given by equation (41). In practice we take the value It was checked that decreasing the time step does not modify significantly the results.
It should be noted that in the evolution problems studied (subsections 4.3 to 4.6), the irreversibility condition for crack propagation given by equation (3) will always be considered. The possible irreversibility of crack propagation without the condition (3) will be investigated separately in subsection 4.7.

Stress distribution at the tip of a penny-shaped crack
We begin with the study of the stress distribution at the tip of a penny-shaped crack without crack evolution. This (static) problem is of interest since it admits an analytical solution which will allow to assess the FFT-based algorithm developed in subsection 3.2.
Let us consider an infinite (isotropic) body containing a penny-shaped circular crack with radius R. We denote by r the radial coordinate and the crack plane is (x 3 = 0). The infinite body is loaded by a uniform remote stress σ ∞ at infinity in the longitudinal x 3 direction, which corresponds to Mode I conditions. The analytical solution for this circular penny-shaped crack of radius R dates back to Sneddon (1946) (see also Broberg (1999) and Cornetti and Sapora (2019)) and the normal stress (σ 33 ) in the crack plane (x 3 = 0) is given by The stress field created by a penny-shaped circular crack in an infinite isotropic medium can be envisaged in a periodic framework with a sufficiently large unit-cell with a prescribed strain tensor of the form For the simulations, an initial penny-shaped circular crack of radius R = 0.25 is considered (see Figure 4). The values d c = 1 and f m = 0 have been chosen so that the material properties are totally degraded within the crack region. First, the normalized stress field σ 33 /σ ∞ calculated numerically (in the simulations σ ∞ corresponds to σ 33 ) is represented in Figure 5 in the planes x 2 = 0 (corresponding to a plane orthogonal to the crack plane) and x 3 = 0 (corresponding to the crack plane). Some numerical artifacts can be observed on the solution close to the crack tip. Those artifacts corresponds to (non-physical) spurious oscillations that arise in FFT-based solutions due to the pseudo-spectral differentiation operator in presence of discontinuous material fields (Morin et al., 2021). This can be easily circumvented by adding smoothness to the material properties field (Morin et al., 2021) (which corresponds in this case to a spreading of the crack tip) or by using other differentiation rules (Willot et al., 2014). In what follows in this paper, the crack-tip is always spread out since t/a = 0 and hence this is not a concern in static.
The numerical solution is then compared in Figure 6a to the analytical solution (45). Since in our case the crack is spread onto 3 voxels, the average stress component σ 33 layer is considered for the comparison. A very good agreement is observed between the analytical solution and the numerical results. It should be noted that the case f m = 0 induces an infinite contrast between the elastic properties of the crack and the bulk, which implies that it is challenging for the FFT-based solver as it requires a large number of iterations with the basic scheme of Moulinec and Suquet (1998) (which is considered in this work). Thus, it is of interest to consider non-null values for f m leading to a finite contrast between the elastic phases. The distribution of the average stress component σ 33 layer is represented in Figures 6b and 6c for the the cases f m = 0.001 and f m = 0.01 respectively. For these values, the effect of the parameter f m is negligible on the stress concentration at the crack tip. Consequently in the following, the case f m = 0.01 will be considered for the sake of practical expedience related to the adopted FFT-based scheme.

Initial crack equilibrium
The previous case was of interest since it has permitted to assess the FFT-based solver. However the initial crack distribution considered previously may not be initially at equilibrium. Indeed, even in absence of macroscopic loading (inducing a microscopic strain field), the energetic driving force e 0 (given by equation (38)) appearing in the crack evolution equation (36) is not null due to the presence of the non-convex local crack energy and crack-tip energy densities, so an initial driving force (5) may initially exist.
First we consider the non-convex local crack energy density η 2 with ω 2 = 1. Since in quasistatic evolution the drag coefficient B does not play any role, only the ratio t/a is relevant. Several values are considered: (t/a)/L 2 = [0, 10 −3 , 1, ∞]. As explained previously, we let the crack field evolve until it reaches an equilibrium; such equilibria are represented in Figure 7 for the cases (t/a)/L 2 = [0, 10 −3 , 1]. In the absence of the non-convex local crack energy term (t/a = ∞) the initial crack completely expands due to the smoothing from the Laplacian term arising from the cracktip energy density, even in the absence of any loading (hence t/a = ∞ is not presented in Fig. 7). Conversely, in the absence of the crack-tip energy density (t/a = 0), the initial crack is already at equilibrium due to the absence of any elastic driving force due to loading, and thus it does not move. In the other intermediate cases (a = 0 and t = 0), the effect of the crack-tip energy density consists in a regularization of the crack tip which spreads on a finite domain. This corresponds to some interphase between the damaged and undamaged areas (the crack-tip region) whose size depends on the ratio t/a. As explained previously, a spreading of the crack tip also improves the calculation of the local stress field by removing numerical oscillations produced by the FFT-based solver (Morin et al., 2021).
The influence of the relative strength of the non-convex, local crack energy density is now investigated in Figure 8, in the case (t/a)/L 2 = 10 −3 . In each case an equilibrium is reached and the shapes obtained are quite close. The sizes of the crack-tips between the damaged and undamaged zones are similar in all cases.

Evolution of a circular crack in mode I
We first investigate the response of a crack under a mode-I tensile loading (normal to the crack surface). An initially equilibrated penny-shaped crack is thus subjected to an increasing macroscopic strain ε = ε 33 (t)e 3 ⊗ e 3 (with ε 33 (t) > 0); a strain increment ∆ε 33 = 5 × 10 −6 was chosen and it was checked to be small enough to capture precisely the onset of failure. The non-convex local crack energy density η 2 has been chosen and the parameters considered for the simulation are given in Table 1.
Elastic degradation functions Local crack energy density Crack-tip energy density Table 1 Parameters considered for the mode-I tensile loading.
The evolution of the (normalized) macroscopic tensile stress σ 33 /µ is represented in Figure  9 and shows three domains: • A linear regime without any evolution of the crack field (between the instants 0 ○ and 1 ○); • A quasi-linear evolution with a stable crack propagation (between the instants 1 ○ and 2 ○). This regime corresponds to an alternation of propagation and equilibrium of the crack tip. The crack slightly propagates after some strain but then stops propagating due to the local elastic unloading: this corresponds to a discontinuous increase of the crack field mimicking a stick-slip behavior.
• A steep softening evolution associated with an unstable crack propagation (after the instant 2 ○). After some critical strain, the crack continuously propagates until it reaches the cell boundary. This defines a threshold after which the crack evolution appears to propagates 'instantly'. To illustrate the different regimes observed in the macroscopic strain-stress curve, the distributions of the crack field c 3 and the normalized stresses σ 33 /µ and σ 13 /µ are respectively represented at several instants in Figures 10, 11 and 12. Furthermore, the (normalized) norm of the FCM crack evolution rate H /max( H ) as well as the quantity |e 0 |/max(|e 0 |) given by equation (38) (the purely energetic crack evolution rate) are represented in Figure  13 (the max are over the spatial domain). Some comments are in order: • The crack field distribution just before the stable crack propagation (see Figure 10a) corresponds exactly to the initial crack field after equilibrium (corresponding to the snapshot 0 ○) which is represented in Figure 7b. This confirms that during the linear regime there is no crack evolution. Thus in this regime, the driving force from elastic energy is balanced by that from the non-convex local crack energy which results in a nil net driving force in the crack evolution equation. This points to the possible existence of a threshold stress to crack motion in the model, much like a Peierls stress for dislocation motion. • After the stable crack propagation regime the crack field shows slight growth (see Figure 10b). The stable regime is characterized by a continuous balance of strain increase (due to the loading conditions) and softening (due to dissipation associated with crack propagation). • The unstable part is characterized by a growth of the crack without an increase of the macroscopic strain (see Figure 10c). At the end of the calculation, the whole crack plane is totally damaged due to periodic boundary conditions. • The distribution of the normal stress follows a similar profile as that investigated in Section 4.2 in the static case. It is worth noting that the stress field at the crack tip is smooth and shows no numerical artefacts, due to the crack-tip energy density term which induces a regularization of the crack front (through a Laplacian in this layer ansatz, with a quasilinear, degenerate 'diffusion coefficient'). It is interesting to note that the crack propagation is associated with the expansion of a large domain at a very low stress level (normal to the crack surface), which results in the softening observed on the macroscopic stress-strain curve. • The distribution of the FCM crack evolution rate H /max( H ) is highly localized at the crack tip (see Figure 13), regardless of the overall situation being close to or strongly out-of-equilibrium. This is consistent with classical notions of crack evolution, while requiring no singularity in stresses or the geometric representation of the crack. This is a direct structural consequence of the theory; the evolution of the crack field is not only of purely energetic/thermodynamic origin, but is additionally constrained by the special kinematics obeyed by the crack tip field which allows evolution only at locations where there is a crack-tip (represented by a non-vanishing value of curl(c)). As shown by the comparison between the two rows of Fig. 13, the purely energetic evolution-rate (as would be operational in a phase field/(L 2 ) gradient flow model) cannot discern the crack-tip in out-of-equilibrium/strongly driven situations. The implications of this feature of dynamic crack propagation in our model, vis-a-vis modeling by phase field models as explained in Agrawal and Dayal (2017), remains to be explored.

Preliminary strength and toughness predictions of the model
We investigate predictions of toughness and local crack-tip stress in the model for onset of unstable crack propagation. For this problem, the stress intensity factor is given by (Cornetti and Sapora, 2019) A set of (initial) crack radii are considered: R/L = [0.15, 0.25, 0.5, 0.75], where L = 1 is the half-period of the microstructure. The macroscopic normal stress σ 33 /µ, the maximal local normal stress σ 33 /µ, and the normalized stress intensity factor K I /(µ t/a) at the onset of unstable propagation are represented in Figure 14 for initially circular cracks (with differing radii).
The macroscopic stress promoting crack propagation increases as the crack radius decreases. However the maximum local stress (located at the crack tip) is almost constant for all crack radii, which suggests that the onset of unstable crack propagation is driven by a local stress criterion (the so-called σ c value). In the case R/L = 0.15 it slightly increases; it should be noted that in this case, the ratio R/h (where h is the thickness of the planar layer) is more important than in the other cases. It appears that smaller R/h ratios beyond a threshold exhibit a different scaling of local strength against crack propagation, but this needs to be studied further. It is also worth noting that the stress intensity factor is bounded by a critical constant, which seems to imply that an apparent toughness (K Ic ) is predicted by the model: crack propagation is triggered when the stress intensity factor reaches a critical value (at least for crack radii above a threshold), which would correspond to some toughness criterion.
These preliminary results of our model appear to be in agreement with the strengthtoughness paradox in brittle materials (Leguillon, 2002): when fracture occurs, both criteria (strength and toughness) appear to be approximately simultaneously fulfilled.

Influence of the ratio a/µ
Since the driving force from the local crack energy density η balances the elastic energy release in the definition of the crack velocity, the ratio a/µ should have an important effect on the onset of unstable propagation. Thus we now study the influence of this ratio. The values a/µ = [10 −8 , 10 −6 , 10 −4 ] are considered and the corresponding evolution of the macroscopic stress σ 33 /µ is represented in Figure 15a. Again, the stress-stain curves are composed of three regimes, (i) a linear regime with no crack propagation, (ii) a stable crack propagation and (iii) an unstable crack propagation. An increase of the ratio a/µ delays the unstable crack propagation occurrence because a greater strain energy is required to exceed the energy barrier from the non-convex local crack energy density: the ratio a/µ thus permits to define the threshold of crack propagation.
It is interesting to note that the maximum value of the stress σ 33 /µ is approximately scaled by a function of the ratio a/µ. Indeed the quantity is approximately constant for the ratios a/µ considered (see Figure 15b).

Crack evolution in shear loadings; Modes II and III
We now consider the response of a circular crack under a shear loading (cf. Lancioni and Royer-Carfagni (2009)). The macroscopic strain is thus of the form ε = ε 13 (t)(e 1 ⊗ e 3 + e 3 ⊗ e 1 ), with ε 13 (t) > 0. Again a strain increment ∆ε 13 = 5 × 10 −6 was chosen and it was checked to be small enough to capture precisely the onset of failure. The parameters given in Table  1 have been chosen.
The evolution of the shear stress σ 13 /µ is represented in Figure 16. As in mode-I tensile, the overall behavior shows three domains: (i) a linear regime without evolution of the crack field (before the snapshot 1 ○), (ii) a quasi-linear evolution with stable crack propagation (between the snapshots 1 ○ and 2 ○) and (iii) an unstable crack propagation (after the snapshot 2 ○).
The distribution of the crack field (in the plane x 3 = 0) and the shear stress σ 13 /µ (in the plane x 2 = 0) at several instants is represented in Figure 17. This permits to highlight the effect of heterogeneity on the stress that is induced by the presence of the crack. An important stress concentration is observed at the crack tip while the domain inside the crack is at zero stress. Furthermore, it is worth noting that the stress concentration (at the crack tip) increases when the crack propagates which likely results in an acceleration of the crack propagation.

Crack evolution in (multi)axial compressive loadings
We consider crack propagation in compressive loadings. There is ample experimental evidence of brittle crack propagation and 'axial splitting' under compressive, uniaxial and longitudinal (Haeri et al. (2014); Li and Wong (2013)), biaxial (Argon, A. et al. (1983); Hoek and Bieniawski (1965); Jaeger and Hoskins (1966); Horii and Nemat-Nasser (1985)), and triaxial (Jaeger and Hoskins (1966)) loadings; also see Wastiels (1979) and Ortiz (1985). On the other hand, it is natural to expect that under transverse uniaxial compressive loading a crack should not propagate. Apart from this latter case and the case of longitudinal, uniaxial, compressive loading where transverse tensile strain on the crack may be expected to be generated, an intuitive explanation for crack propagation does not seem to be apparent. In fact, even in the case of longitudinal, uniaxial (or 'polyaxial' (cf. Jaeger et al. (2009)) compressive loading, say in a cuboidal domain, longitudinal crack propagation resulting in axial splitting does not fit within the tenets of classical fracture mechanics, even though seen to occur in nature, as explained in Lajtai et al. (1990). This is because for such loading all of the Modes I, II, III stress intensity factors vanish (either because of the absence of far-field stresses or due to them being compressive). Further-more, in laboratory experiments on brittle solids, single wing cracks, progressing from a single notch like flaw oriented slightly away from the direction of major compressive stress, stop propagating after becoming parallel to the direction of major compressive stress and traveling distances of ∼ 10 − 15 times the size of the initial crack (Hoek and Bieniawski (1965); Argon, A. et al. (1983); Horii and Nemat-Nasser (1985)). Notwithstanding this, axial splitting under compressive load is an observed phenomenon even if its exact cause is still not understood clearly (to our knowledge), being attributed to collective effects of many microscopic cracks (e.g. Horii and Nemat-Nasser (1985)) as one example, or composite material effects Ortiz (1985). In the following results in this section, we demonstrate stable crack propagation leading to axial splitting with increasing longitudinal compressive load in our simple model, without commitment to microscopic details of whether our crack is a single 'Griffith' crack or a notch of finite width incorporating, in a gross fashion, the interaction of many axial microcracks.
We begin with providing some approximate reasoning of why it may be within the realm of possibility in our model to predict such observed behavior for cracks under compressive loadings, based purely on energetic considerations. We then treat each case in turn, and qualitatively correlate our results with experimental observations.
From the definition of the elastic energy (21), it appears that crack propagation could be triggered with a local triaxial compressive loading. If we exclude, for simplicity, the contribution of the local crack and crack-tip energy densities on crack propagation, one has to study the local change of elastic energy at the crack tip. Let us consider near the crack tip a local strain tensor of the form where we assume that the crack flanks are in compression or unstrained in the transverse direction (i.e., ε 33 ≤ 0). For this loading, the strain energy reduces to 2ϕ E = 2(C 1111 + C 1122 )ε 2 11 + 4C 1122 ε 11 ε 33 + C 3333 ε 2 33 .
Assuming that (i) elastic properties are totally degraded within the crack and (ii) λ = 2µ, the change in elastic energy between cracked and non-cracked region is solely driven by the quantity 3ε 2 11 + 2ε 11 ε 33 .
(52) Thus a crude first estimate of energetic advantage for the crack to grow (corresponding to a greater elastic energy density in intact material next to crack tip to be rendered vanishing by crack propagation) is 3ε 2 11 + 2ε 11 ε 33 > 0 so that an approximate necessary propagation condition would reduce to ε 11 < 0 and ε 11 < −2ε 33 /3.
Therefore it may be expected that a fully compressive loading (ε 11 < 0 and ε 33 ≤ 0) can promote crack propagation. It is also easy to see by essentially the same argument that crack propagation may be expected in the compressive uniaxial and longitudinal, and biaxial loading cases as well.
Uniaxial, transverse compressive loading (mode I compression). We consider a circular crack under a uniaxial compressive macroscopic strain ε = ε 33 (t)e 3 ⊗ e 3 (with ε 33 (t) < 0), no crack propagation is observed, even for very high values of ε 33 . In that case, the normal stress σ 33 is uniform in the whole domain (including the crack) since the crack compressive behavior reduces to the intact material due to the non-linearity of the elastic behavior considered. Indeed for this loading, there is no crack propagation because there is no energetic advantage for the crack to grow since there is the same elastic energy in the crack and outside of the crack.
'(FFT) Brazilian test' with a confined crack. The Brazilian test is a widely used testing method to obtain the tensile strength of rocks and concrete (see the review by (Li and Wong, 2013), Haeri et al. (2014)). It consists in the uniaxial, diametral compression of a circular disk. Failure in these specimens inevitably results from the splitting of the disk about the diametral loading plane, most likely by the propagation of an interior extensional crack in the direction of the (compressive) loading; this is related to tensile strains which develop in the Brazilian specimen due to the Poisson effect. While our employed FFT based methodology is not suited to model this test exactly, we consider an analogous loading protocol of the form of a cuboidal domain with a pre-existing interior crack in the x 3 = 0 plane. The macroscopic transverse, normal stress σ 33 is close to zero; this also implies that there cannot be a Mode I stress intensity at the crack tips in this problem. The evolution of the macroscopic stresses, which are represented in Figure 18, confirm that σ 33 is zero. Stable crack propagation occurs between the instants 1 ○ and 2 ○. The line plot of the local stress σ 33 is shown in Figure 18b at the beginning of the stable crack propagation (instant 1 ○) in the line x 2 = 0 of the plane x 3 = 0. Small oscillations are observed which are due to the continuous Green operator but it was checked that the maximal and minimal values were correctly predicted (using a finer spatial discretization). Near the crack tip, the local normal stress σ 33 is positive which means that locally the specimen is subjected to a tensile loading. This confirms that local tension is the mechanism promoting crack propagation in the Brazilian test (Li and Wong, 2013).
Through crack under longitudinal compression. Next we consider the case of a through crack subjected to the same loading as in the previous case. This numerical experiment is motivated by an experimental results of (Haeri, 2015, Fig. 5(a)). The specific experiment of interest to us involves a through-notch in the plane spanned by a diameter and a generator of a cylindrical specimen of Portland Pozzolana cement; the cylinder is compressed in the axial direction and crack propagation from the notch is observed to progress in the plane of the notch.
In our calculation, the crack is supposed to be infinite in the direction x 2 and its width in the direction x 1 is equal to 0.5.
The evolution of the longitudinal stress σ 11 /µ is represented in Figure 19. In addition, the distribution of the crack field is represented in Figure 20. As in the previous Brazilian case, crack expands due to the presence of tensile stress within the crack layer.
Triaxial compressive loading. This numerical experiment is motivated by experimental results on fracture under biaxial and triaxial compressive loadings (Argon, A. et  al. (1983); Hoek and Bieniawski (1965); Jaeger and Hoskins (1966); Lajtai et al. (1990)). Restricting crack motion to a plane, in-plane crack propagation under all-round, average, compressive stressing is demonstrated. Such crack propagation has been observed for specimens with pre-existing through cracks under biaxial loading (see e.g., Argon, A. et al. (1983); Hoek and Bieniawski (1965), and failure, of externally intact rock cylinders/disks (presumably containing pre-existing cracks in their interior), in the triaxial loading experiments of Jaeger and Hoskins (1966). Such behavior apparently also abounds in nature (Lajtai et al. (1990)).
The evolution of the longitudinal and normal stresses (σ 11 /µ and σ 33 /µ) are represented in Figure 21. It should be noted that for ε 11 (t) = 0 the macroscopic stresses are not null due to the axial strains ε 22 and ε 33 . In that case the evolution of the macroscopic normal stress σ 33 /µ does not show any softening. It is interesting to note that in this case, the crack propagates (see Figure 22) but this does not lead to a brutal softening in the macroscopic stress-strain curves. In the direction normal to the crack, the loading is compressive (ε 33 < 0), which implies that the stress σ 33 is almost uniform in the whole domain including the crack area. The distribution of the longitudinal stress, which is represented in Figure 23, is heterogeneous. In this case, the crack has an energetic advantage to propagate due to the elastic energy produced by the longitudinal strain. However this crack evolution is not associated with a softening because the compressive loading is parallel to the crack surface resulting in very small loss in load carrying capacity, even for a through crack. Of course, in the transverse direction, even the through crack responds in an 'intact' manner under a compressive strain field.

Crack irreversibility
We investigate crack propagation and its (possible) irreversibility without using the irreversibility condition (3). When (3) is dropped, crack healing is not a priori prevented in the evolution equation (1) 4 . Indeed, consider a situation where there is no external loading in a body and it has a crack specified by the initial condition. If the local crack energy density function is η 1 , endowing the cracked region with 'Griffith' type surface energy, then, purely on energetic grounds, the evolution should heal the crack to rid the body of its energy content unless prevented from doing so by the ad-hoc imposition of the irreversibility condition (we ignore the effects of the crack-tip energy density for the essential argument here which only regularizes the crack-tip). On the other hand, for a local crack energy density described by η 2 a completely damaged crack region incurs no energetic cost and there is no energetic impetus to heal. We illustrate these features in the following numerical experiments.
Equilibrium solution for η 1 . We begin with the 'Griffith' type local crack energy density η 1 (with ω 1 = 0.5). In that case no equilibrium position can be reached for all ratios t/a. When the latter ratio is non-vanishing, the crack-tips tend to smear out while the local energy density tries to cover as little area of the layer by a crack. The overall result is the diminution of the magnitude of the crack field with time in the layer. When t/a = 0, the process is governed entirely by the energetic healing discussed above; the evolution of the crack field (of the initially penny-shaped circular crack considered in subsection 4.3) during this process of equilibration is represented in Figure 24, at several instants before its complete healing/annihilation. Unloaded crack equilibria with η 2 . next we study the local crack energy density function η 2 , again for an initially penny-shaped circular crack. In the absence of the irreversibility condition, equilibrium positions can be obtained in contrast with the Griffith energy. As expected, the equilibrium profiles are sensitive to the ratio t/a: spatially localized, crack-like equilibria are attained for (t/a)/L 2 < 10 −4 . For values (t/a)/L 2 > 10 −4 the regularizing 'diffusion' from the crack-tip energy density overwhelms the localizing influence of the non-convex local crack energy density, leading to a complete spreading of the crack field in the layer. For illustrative purposes, the equilibrium position in the case (t/a)/L 2 = 10 −5 is represented in Figure 26a.
Equilibrium after loading-unloading. A loading-unloading cycle is now considered in order to investigate crack irreversibility. Thus, we consider a Mode-I tensile loading without the irreversibility condition (3). The function η 2 is chosen because it allows to define an initial equilibrium position as shown previously; the parameters considered for the simulation are given in Table 2. In order to study the possible irreversibility of crack propagation, we consider a progressive loading until the crack slightly propagates inside the stable domain, and then an unloading is applied until the macroscopic stress reaches zero (see Figure 25).
Elastic degradation functions Local crack energy density Crack-tip energy density f m d c a/µ ω 2 (t/a)/L 2 0.01 1 10 −4 1 10 −5 Table 2 Parameters considered for the mode-I tensile loading.
The distribution of the crack field is represented at several instants in Figure 26. As expected the crack has slightly growth after the loading (Figure 26b). Interestingly, the crack does not close up during unloading (Figure 26c) which suggests that with this model the evolution of the crack is irreversible (up to some local fluctuations).

Crack evolution in dynamic loadings
We demonstrate the dynamic FFT algorithm for elastodynamics of Section 3.3 by undertaking a preliminary exploration of the interaction of rate of loading, the effect of material inertia, and the crack mobility in FCM. The physical problem that most closely resembles our simulations is that of a region containing a fault in its interior with a pre-existing rup-ture zone, with the whole system being 'slowly' loaded by a spatially uniform strain-rate, as may be envisaged in the loading of tectonic plates (admittedly, our Mode I loading is not entirely realistic for this specific physical setting). Here, the rupture zone and its possible propagation is treated as a crack with an emergent slip-weakening response not put in 'by hand.', without getting into a discussion of the appropriateness of such a ruptureas-a-crack assumption (cf. (Zhang et al., 2015, Sec. 8). Depending on the nondimensional parameter B • characterizing the mobility of the crack tip, interesting crack propagation behavior in the fault spanning sonic thresholds of the adjoining elastic material are observed. Unlike crack propagation driven by dynamic loading applied to the boundary of the domain, the present situation provides a large reservoir of stored elastic energy almost everywhere in the domain which lends itself to dissipative release by crack propagation. Due to this available supply of energy in front and behind the crack-tips (as opposed to energy supply arising from a stress pulse brought to the crack-tip via loading from a boundary), crack motion beyond sonic speeds in our results is perhaps not so surprising.

Description of the simulations
Simulations with dynamic loadings can be performed by combining the fully dynamic explicit solver for the equilibrium equations and the Godunov-type solver for the crack propagation equation. In this work, we will consider only loadings with uniform macroscopic strain ratesε (so that the macroscopic accelerationε is null). In practice, we just apply successive macroscopic strain increments ∆ε which are given by where ∆t is the time step used in central scheme for solving the crack evolution problem. It is chosen in order to check the following stability constraint where v l is the speed of propagation given by equation (43) and v d is the dilatational wave speed given by In practice we take the following value: The procedure used to evolve the crack field during quasi-static calculations (see Figure  3) is not needed in the dynamic case because, in the dynamic setting, the crack field is updated for each strain increment. Crack propagation is expected to be slowed down when inertia effects are present.
In the following, we consider the same conditions than previously. A [−1, 1] × [−1, 1] × [−1, 1] unit-cell discretized with 128×128×128 pixels so the spatial step is ∆x = 0.015625. The planar layer is again supposed to be centered in the plane x 3 = 0 and a constant thickness h = 0.046875 is considered for the layer, which corresponds to a number of 3 pixels. An initial small square-shaped crack is considered in all the simulations in order to investigate crack propagation in a large domain (see Figure 27). The following value is taken for the elastic properties: λ/µ = 2.

Overall response in mode-I dynamic loadings
We focus on a mode-I dynamic loading so that the macroscopic imposed strain is of the form ε = ε 33 (t)e 3 ⊗ e 3 (with ε 33 > 0). The parameters related to the non-convex local crack energy and crack-tip energy densities are given in Table 1.
In these dynamic simulations, different behaviors are expected to manifest according to the value of the drag parameter B. Indeed, (and explicitly resorting to circled and uncircled quantities representing corresponding non-dimensional and dimensional ones, respectively) the order of magnitude of the crack velocity v c is about so that the ratio between the crack velocity and the elastic shear wave velocity We will consider two values for the nondimensional drag parameter B • = (BL)/ √ µρ = [10 −4 , 10 −1 ], which approximately lead to a ratio v c /v s ∼ [31.6, 0.0316]. Thus the value B • = 10 −4 is expected to correspond to supersonic crack behavior (noting λ/µ = 2), while the value B • = 10 −1 would correspond to subsonic crack behavior.
In order to investigate the effect of dynamics, we will consider several values for the macroscopic strain rate. The evolution of the (normalized) overall tensile stress σ 33 /µ (where σ 33 = σ 33 , · being the spatial average over the unit-cell Ω) is represented in Figure 28.
The effect of inertia is, as expected, a retardation of the softening occurrence. When the strain rate is high, the crack propagation is limited by 'viscous' effects emerging from the reciprocal crack mobility, B. If the intrinsic time scale set by B is large compared to that set by the loading rate and elastic wave propagation, then there is not sufficient time available for stress/energy relaxation by crack propagation to occur. Computationally, for a given strain increment, the crack is not evolved until it reaches some equilibrium position as in the quasi-static case; it is simply evolved at each strain increment (without internal iterations to achieve equilibrium of the crack). Thus the crack 'takes more time' to propagate in contrast to the strain.This implies that the evolution of the macroscopic stress follows only two regimes: • A linear evolution with no evolution of the crack field. In both cases, this corresponds to the evolution before the snapshot 1 ○; • A non-linear evolution with the propagation of the crack. This can lead to an apparent quasi-linear behavior when the loading increases more rapidly than the crack, but softening eventually takes place due to a build-up of elastic energy that results in a large driving force which overcomes the relatively larger drag, at which point the crack rapidly propagates through the entire domain.

Subsonic behavior: B • = 10 −1
In order to understand the effect of strain rate on the crack propagation, we study at several snapshots the distributions of the crack field c 3 , the evolution rates, the hydrostatic stress σ h /µ (where σ h = Tr(σ)/3) and the shear stress σ 13 /µ for the strain rateε 33 = 2 × 10 −4 , which are respectively represented in Figures 29, 30, 31 and 32. It must be noted that the the distributions of the crack field c 3 , the hydrostatic stress σ h /µ and the shear stress σ 13 /µ have not been represented after the snapshot 6 ○ because of the periodic boundary conditions. During propagation, the initial square-shape crack becomes a circular crack. As in the quasi-static case, the distribution of FCM crack evolution rate H /max( H ) is highly localized at the crack tip (Fig. 30), in contrast to the energetic evolution rate |e 0 |/max(|e 0 |) which is significantly spread out -the latter would result in crack evolution behind the tip in dynamic and out-of-equilibrium situations, in a model where this was the 'driving force'. Both elastic shear and pressure stress waves propagate ahead of, and behind, the moving crack-tip. The stress level inside the crack tip is lower than in the bulk, but it is not nil, this being an artifact of the value of the parameter f m = 0.01. It is expected that the stress level inside the crack can be systematically reduced with a decrease of f m (bounded away from 0). Next we study the distributions of the crack field c 3 , the hydrostatic stress σ h /µ and the shear stress σ 13 /µ for the strain rateε 33 = 2 × 10 −4 and a significantly lower value of the nondimensional drag, which are respectively represented in Figures 33, 34 and 35. Before crack propagation, the distributions of the normal and shear stress (at the snapshot 1 ○) have the same features as those of the quasi-static case. At the very beginning of the crack propagation (snapshots 2 ○ and 3 ○), spherical stress waves are produced at the crack tip and propagate in front and behind the crack tip. This suggests that, during the early crack propagation, the crack velocity is in a subsonic regime. As the overall strain increases, the stress waves ahead of the crack-tip gradually disappear and totally vanish. Mach cones of hydrostatic stress, including 'stress-release' waves in the bulk arising from propagating damage in the fault, can be seen corresponding to snapshots 5 ○ and 6 ○ in Fig. 34, characterizing supersonic behavior. Mach cones in shear stress are also observed in the field plots corresponding to snapshots 5 ○ and 6 ○ in Fig. 35, with the shear stress within the propagating crack layer small/vanishing because of the damage due to tensile loading. The shear stress Mach 'wings' are of opposite sign across the crack, as dictated by symmetry.
6 Towards the fracture of materials with random microstructures As a final example, we study the possibilities of the FFT-FCM framework to handle more general situations with complex distribution of the crack field. Such capability is of practical interest because the fracture properties of real materials are expected to depend significantly on the microstructure. Indeed, in brittle materials the fracture properties are usually dominated by the size and spatial distributions of flaws or microcracks.
We begin with the generation of a microstructure with random pixels (comprising ∼ 2.2% of the crack plane surface) having c 3 = 1 , which is represented in Figure 36a. We investigate the response of this crack field under a mode-I loading with simulation parameters given in Table 1.
First the equilibrium crack field starting from the above initial microstructure is computed without applied load and is represented in Figure 36b. Several micro-cracks have been created in this process due to the coalescence of several initially close 'voxels'. The surface fraction of damage is ∼ 7.7% after equilibration.
Then the equilibrated, no-load crack microstructure is subjected to an increasing macroscopic strain ε = ε 33 (t)e 3 ⊗ e 3 (with ε 33 (t) > 0). The macroscopic stress-strain curve is represented in Figure 37. Overall it has, for this specific realization, the same features as that of the case of a body with a single crack, that is (i) a linear regime with no crack evolution (before snapshot 1 ○), (ii) a stable crack propagation regime (between snapshots 1 ○ and 2 ○) and (iii) an unstable crack propagation regime (after snapshot 2 ○). It must be noted that the stable crack propagation regime, although it has less influence on the macroscopic curve, comes quite early compared to the unstable regime.
The distribution of the crack field and normal stress are represented, at the snapshots 2 ○ to 4 ○, in Figures 38 and 39, respectively.
At the end of the stable regime (see Figure 38a), new clusters have been formed which suggests that the stable regime corresponds to local organization of the micro-cracks similar to that observed during no-load equilibration. As shown in Figure 39a, the (normalized) local stress level during this regime is close to 8 × 10 −3 which is of the same order of magnitude as the results of Section 4.4. The unstable crack propagation is characterized by a growth and coalescence of the biggest micro-cracks (see Figure 38b-c) which is associated with a decrease of the local stress (see Figure 39b-c).

Concluding remarks
We have developed and demonstrated a pde-based tool for the modeling and analysis of (dynamic) brittle fracture, restricted here to crack propagation in a single fault layer. Crack-path selection within the fault as well as strength and toughness criteria are pri-mary emergent features of the model. Various other features of brittle fracture have been demonstrated, indicating preliminary promise as a general-purpose tool for analysis, and hopefully design, in fracture mechanics. The tool has the flexibility of accommodating experimental input from observed crack velocities and shear and normal contact stiffnesses within cracked regions, if so desired. We have intentionally focused on zero-surface energy cracks where all elastic energy released on crack propagation is dissipated (except for what is stored at crack-tips), but the capability allows for situations where substantially larger amounts of released elastic energy can be stored in cracked regions behind the crack-tip (as in a 'Griffith crack').
Beyond this first analysis, much remains to be done in realizing the full potential of our theoretical model. This includes, first and foremost, demonstrating a computational capability without the planar fault restriction, allowing for general crack pattern evolution in 3-d. This is a challenging endeavor, because of the novel 3 , multi-dimensional, essentially first-order Hamilton-Jacobi transport system for crack evolution, singularly perturbed by a degenerate, quasilinear, second-order, parabolic 'regularizing' term, all entirely dictated by the mechanics and physics of the model. Other comparatively straightforward extensions correspond to incorporating ductile fracture and the mechanics of finite deformation. Finally, we mention the recent work of Kumar et al. (2020) related to crack nucleation in phase field models. This is essentially achieved in the said model by incorporating a phenomenological yield-like criterion reflecting the 'strength' of the material. Our work, restricted to small deformations is silent on the matter of nucleation; however, as observed in Acharya (2018), based on the results of Garg et al. (2015) and the similarities of the FCM model to FDM, especially in the nature of its evolution equation for the crack-tip field to that of the dislocation density field, it appears to be a reasonable hypothesis that crack nucleation may be a natural instability in the finite deformation version of FCM, involving the limits to elastic strength as encoded in the material velocity field through the solution of the relevant nonlinear elastic (i)bvp. We leave such questions for future work.
(2) it is assumed that the elastic modulus C(x, t) has major and minor symmetries and is positive semi-definite.
The goal is to determine a set of conditions that guarantee uniqueness of solutions within a class of functions of the form where H is a specified second-order-tensor-valued function of time alone, x 0 ∈ Ω is an arbitrarily chosen fixed point, and u * belongs to a class of functions on Ω × [0, T ] whose properties will be subsequently specified as part of our analysis.
Consider two solutions u 1 and u 2 of (62) in the class (63) and denote the difference u 1 − u 2 =: u and similarly the difference strain ε := (∇u 1 ) sym − (∇u 2 ) sym . We then have that where t := t 1 − t 2 = (C (ε * 1 − ε * 2 )) n, u :=u 1 −u 2 =u * 1 −u * 2 are the difference of the tractions produced by the two solutions u 1 and u 2 on ∂Ω and the difference of the corresponding velocities, respectively. If it can be ensured (3) that for any two solutions in the class defined by (63), the first term on the righthand-side of (65) vanishes; (4) Ω (·) :Ċ(··) dx, viewed as a bilinear operator on symmetric second-order tensor fields on Ω is negative semi-definite; (5) Ωρ (·) · (··) dx, viewed as a bilinear operator on vector fields on Ω is negative semidefinite; and (6)u(x, 0) = u(x, 0) = 0, then E(t) ≤ E(0) = 0 for all t ∈ [0, T ] 4 . But by the positive semi-definiteness of C(x, t) and the positiveness of ρ(x, t) in Ω × [0, T ], 0 ≤ E(t) for all t. Thus we have E(t) = 0 for all t ∈ [0, T ]. Due to the positiveness properties of the mass density and elastic modulus fields, the kinetic energy and the strain energy terms in the definition of E(t) in (64) have to individually vanish. Then, assuming that the kinetic energy density 1 2 ρu ·u is a continuous function, we have thatu(x, t) = 0 for all (x, t) ∈ Ω × [0, T ]. But then the initial condition u(x, 0) = 0 implies that u 1 = u 2 for all (x, t) ∈ Ω × [0, t], and we have uniqueness.
With reference to Fig. 40, we now consider Ω to be a domain whose external boundary can be decomposed into pair-wise opposing surfaces which are parallel-translates of each other (a rectangle or a cuboid fits this definition). We now think of the functionsu * , .. C, and ∇u * sym = ε * being constrained to attain equal values at points that are parallel translates of each other on these pairwise-parallel surfaces comprising ∂Ω for all times (a cuboidal domain with u * represented by a Fourier series on it satisfies such conditions). These requirements ensure that the stress field on Ω σ(x, t) = C(x, t) (ε * (x, t) + H sym (t)) also satisfies such 'periodic boundary conditions' on ∂Ω and therefore the traction field on ∂Ω satisfies 'anti-periodic boundary conditions' (because of the parallel but oppositely pointing normals at each pair of parallel-translated points on each pair of parallel surfaces of ∂Ω). It is clear then that under these stipulations, condition (3) is satisfied. This is so because the tractions arising from two solutions are individually 'antiperiodic' and therefore their difference is as well, and while the individual velocities (given byu i (x, t) = u * i (x, t) +Ḣ(t)(x − x 0 ), i = 1, 2) are not 'periodic', their difference is. Furthermore, requiring any solution to satisfy a specified displacement initial condition that can be represented in the form u * 0 (x) + H(0)(x − x 0 ) and specified velocity initial condition of the form v * 0 (x) +Ḣ(0)(x − x 0 ) where u * 0 and v * 0 satisfy 'periodic boundary conditions,' ensures the satisfaction of condition (6). Thus for problems defined by a given function of time H and the structural conditions (1), (2), (5), and (4) defined on a domain of the type just discussed, solutions belonging to the class (63) satisfying initial conditions of the type just discussed, are unique.
Corollaries: When H = 0, we have purely periodic solutions, driven solely by the initial conditions. Elastodynamics allows infinitesimally rigid unique solutions; thus, skew symmetric H are allowed. 'Periodic' elastostatics is recovered for ρ ≡ 0.
Remark 1 As shown in this paper in Secs. 3.3 and 5, at least for the finite dimensional approximation defined by (62) and (63) with u * represented by a Fourier series with a finite number of terms, solutions can be constructed whose values at the final time T are consistent with (63), providedḦ = 0 in (0, T ). This suggests that a simulation spanning the intervals [0, T 1 ] and [T 1 , T 2 ] with two different 'forcing' functions H 1 (defined on [0, T 1 ]) and H 2 (defined on [0, T 2 − T 1 ]) can be conducted with our technique provided H 1 (T 1 ) = H 2 (0) withḢ 1 andḢ 2 taking arbitrary constant values in the intervals [0, T 1 ] and [T 1 , T 2 ], respectively. The final values of the displacement field and the velocity field from the interval [0, T 1 ] serve as the initial conditions for the next interval; there are velocity discontinuities at the 'joins.' Such ideas can be used to approximate 'arbitrary' loading histories within this periodic-in-space framework.
Remark 2 In reality, the space-time dependence of the elastic moduli C arises from a dependence on the evolving displacement field through the coupling to the crack evolution. This suggests the study of conditions of well-posedness for the fully coupled problem, which is beyond the scope of the present paper.