FFT-based interface decohesion modelling by a nonlocal interphase

In this paper, two nonlocal approaches to incorporate interface damage in fast Fourier transform (FFT) based spectral methods are analysed. In FFT based methods, the discretisation is generally non-conforming to the interfaces and hence interface elements cannot be used. This limitation is remedied using the interfacial band concept, i.e., an interphase region of a finite thickness is used to capture the response of a physical sharp interface. Mesh dependency due to localisation in the softening interphase is avoided by applying established regularisation strategies, integral based nonlocal averaging or gradient based nonlocal damage, which render the interphase nonlocal. Application of these regularisation techniques within the interphase sub-domain in a one dimensional FFT framework is explored. The effectiveness of both approaches in terms of capturing the physical fracture energy, computational aspects and ease of implementation is evaluated. The integral model is found to give more regularised solutions and thus a better approximation of the fracture energy.


Introduction
Polycrystalline materials at a microscopic level show clear heterogeneous deformation patterns. This heterogeneity arises from the locally fluctuating mechanical properties of different phases and differences in lattice orientations between different grains.
FFT based spectral solver was originally introduced to model the mechanical behaviour of composite microstructures [1]. Since then it has emerged as a promising tool for modelling the micromechanical response of polycrystalline materials [2,3]. A comparison of different FFT formulations and solution approaches in a crystal plasticity constitutive framework [4] was presented in [5]. Recently, an FEM perspective on an FFT based spectral formulation for small strain non-linear material behaviour was given in [6] and its extension to a finite strain setting was presented in [7]. Alongside such improvements, much effort has gone into making the method suitable for various applications. The computational efficiency of FFT methods makes them attractive to solve multi-field problems, for e.g. a nonlocal crystal plasticity formulation [8], ferroelectric switching [9], etc.
Polycrystalline microstructures show a combination of intergranular and intragranular damage. Gradient based nonlocal damage simulations in a polycrystalline material were performed using an FFT approach in [10] without differentiating between bulk damage and interface decohesion. Interfaces like phase boundaries and grain boundaries are . a 1 shows a polycrystal with interface elements embedded at sharp boundaries. a 2 depicts an interface conformal finite element mesh with cohesive zones that open during loading according to the traction (σ )-separation (δ) law shown in a 3 with characteristic fracture energy G c . b 1 shows the same polycrystal with interfacial bands representing the sharp interfaces. An example constitutive behaviour for the material points (within the volumetric interphase) assigned dissipation density g c is shown in b 3 particularly susceptible to damage initiation and propagation. For example, damage nucleation in Dual Phase steel has been attributed to interface decohesion at ferrite-martensite phase boundaries and martensite cracking along prior austenite grain boundaries [11]. In our contribution, we focus on the computational treatment of interface decohesion and its interaction with the (undamaged) bulk behaviour in FFT solvers.
In finite element (FE) based solution schemes, such problems are generally tackled using cohesive zones [12,13], see Fig. 1 (a 1 , a 2 , a 3 ). Interfaces are identified a priori and interface elements are introduced at these physical locations. This is possible since the FE discretisation is made to conform to the sharp interfaces. The opening behaviour of these interface elements is governed by a traction-separation law (TSL) that is assigned to them. FFT based spectral methods, however, rely on a regular grid and a method to include sub-dimensional (e.g. planar in 3D) elements cannot be used. Composite pixel based approaches, for example [14], which make use of the interface normal to assign homogenised mechanical properties to grid points near the interfaces, are a step in that direction. But at this point it is not clear how they will perform in situations of damaging interfaces that are inclined or curved. The high mechanical contrast combined with inherent pixelation effects of FFT may cause stresses relayed across damaging interfaces, which is undesirable for the modelling of cracks. This limitation is dealt with using the idea of an interphase (volumetric) band, as depicted in Fig. 1 (b 1 , b 2 , b 3 ). The material points in the vicinity of the sharp interface are identified and furnished with a damage constitutive behaviour (softening) and corresponding kinematics-in addition to the deformation mechanisms attributed to them as a part of the respective grains they belong to. We aim to have multiple points across the thickness of the interphase. We expect this to reduce the pixelation effect and ensure that the two crack faces are fully decoupled. This approach allows capturing the interfacial mechanics and still benefits from the computational (memory and speed) and implementation related advantages of FFT based methods. Such an approach to represent interfaces as interphase has been used previously in the literature. Hsueh-Hung et al. [15] used it to assign different material behaviour to grain boundary regions than the bulk. Their approach was physically motivated and directed towards understanding metal plasticity in the nanocrystalline regime. Clayton et al. [16] used this approach within a finite element implementation of phase field damage to understand the competition between intergranular and intragranular damage.
In volumetric dissipation based models, one specifies a dissipation density and an internal length scale to the model-which for the present case is the band thickness, l. The total dissipation resulting from the volumetric damage process must equal the physical dissipation of the real sharp interface (or cohesive zone model). Accordingly, the dissipation density has to scale inversely with the band thickness. A straightforward inversely proportional relationship results if the entire interphase band damages uniformly. This is not guaranteed in a conventional local softening model due to localisation and lacking objectivity with respect to the grid spacing. Various regularisation strategies have been used in the literature to remedy this. They can broadly be classified into two categories-integral based averaging [17] and gradient damage based regularisation [18][19][20][21][22].
The objective of this paper is to analyse the applicability of these two regularisation approaches in an FFT method to model the softening interphases. We compare the effectiveness of the two approaches for the idealized case of one dimensional (1D) bar with a single interphase band. The constitutive response of the underlying (bulk) material is assumed to be linear elastic. This simple case allows us to systematically study the effect that the nonlocality has on the damage evolution-with emphasis on the desired objectivity of the predicted fracture energy with respect to the interphase band thickness. The organisation of this paper is as follows. In "Problem statement" section the 1D problem considered in our study is stated and the constitutive model and kinematics are described. Next, the gradient damage and integral averaging approaches are discussed in "Gradient based nonlocal damage" and "Integral based nonlocal damage" sections, respectively. In "Results and discussion" section, results of both the approaches are discussed. A conclusion and outlook is presented in "Conclusion and outlook" section.

Problem statement
Consider a 1D bar of length L as shown in Fig. 2. Let us assume a sharp interface at the center (x = 0) of the bar. To capture the response of this interface using an FFT based spectral scheme, the sharp interface is substituted by a volumetric interphase band of thickness l spanning the region from x 1 = − l 2 to x 2 = l 2 . The elastic properties are taken to be homogeneous, with Young's modulus defined as E. To model the kinematics of decohesion, a damage strain field, is the sub-domain representing the interphase region in the total domain Ω = (−L/2, L/2). Since only the interphase region accommodates the damage strain, the damage constitutive model is only used in Ω i . In Ω i it is superimposed on the elastic strain ε e (which exists on Ω) as, A local damage field φ l (x), for x ∈ Ω i , is introduced. This field represents the material integrity. It is initialised to unity, representing a perfectly bonded state. When the damage evolves it drops towards zero, thereby reducing the capacity of a material point to transmit stresses. A damage constitutive model exploiting only the local field φ l is ill-posed and provides non-physical results due to a pathological tendency to localisation [17]. Therefore, a nonlocal counterpart field, φ nl , is defined also for x ∈ Ω i . The bar is subjected to an overall monotonic uniaxial tensile loading at a constant ratė ε. This loading rate is related to the local strain rate field aṡ The stress at equilibrium has to satisfy,

Constitutive model
The constitutive model relating the stress and the elastic strain for all material points at The evolution law for ε d is taken aṡ where σ c φ nl 2 is the effective opening resistance, with σ c the critical stress of the model; ε d0 is the reference opening rate. The Macaulay brackets, defined as · = (· + |·|) /2, make sure that the damage strain does not evolve until σ reaches the initial critical stress σ c . The nonlocal damage φ nl , in Eq. (5) is obtained by appropriately regularizing the local damage field φ l defined as The critical damage strain ε f governs the softening slope. It is noted that the rate dependent form depicted in Eq. (5) is motivated from the sub-scale mechanics (micro-void formation, plasticity etc.).

Analytical solution for the rate-independent case
Since the applied loading is monotonic tension, one can rewrite Eq. (5) beyond the damage threshold as In order to retrieve an analytical solution, a special case is considered. It is assumed that the nonlocal field coincides with the local field, which is furthermore assumed to be uniform within the interphase band. In this case we have, also using Eq. (6), In the rate-independent limit, implyinġ Equation (8) yields Due to the equilibrium condition in this 1D setting, a spatially constant stress field results. The average strainε for any instant after the damage has initiated can be obtained by inverting Eq. (10) to give ε d as a function of σ and integrating the elastic and damage contributions to the strain field: Thus, the overall stress-strain response for this rate-independent uniform damage case can be written in piecewise form as The dissipation per unit volume due to the damaging process can be computed analytically as where use has been made of Eq. (10).
Relating the dissipation throughout the interphase volume to the expected interfacial dissipation G c , we obtain Accordingly, the product of the damage constitutive parameters σ c and ε f should depend inversely on the interphase band thickness. This equality, together with Eq. (13), indicates how the strain to failure needs to be scaled with respect to the interphase band thickness l, in order to recover the strength and dissipation of a cohesive zone, for arbitrary l:

Gradient based nonlocal damage
In the differential form of nonlocality, the nonlocal damage field is determined as the solution of a differential equation. Here, the classical gradient based damage equation is considered to retrieve the nonlocal counterpart of φ l , which reads in one dimension: where λ(x) is the nonlocal length scale. In principle, the damage variables φ nl and φ l are meaningful only for x ∈ Ω i , and therefore Eq. (16) for the nonlocal field φ nl should be solved on the sub-domain Ω i only, with flux free boundary conditions (also referred as interface condition) on ∂Ω i , i.e.
where ∂Ω i = {x 1 , x 2 }. Since FFT solvers are based on global and periodic shape functions with a regular discretisation grid, the problem has to be solved on the full domain instead, still respecting Eq. (17) at least approximately. This may be achieved by extending Eq. (16) to the entire domain Ω, but with a piecewise constant λ, equal to the desired value λ in inside the band and a much smaller value λ out outside the band. In that case, the interface condition at the edges of the interphase reads: which corresponds with Eq. (16). In the limit of λ in /λ out → ∞ this condition enforces Eq. (17). In practice, a finite, but large, contrast λ in /λ out is used and Eq. (17), the flux free boundary condition, is satisfied only approximately. Phase field approaches also provide a differential equation based continuum solution for evolution of damage. In principle, a phase field model based on a Ginzburg-Landau approach can be used as well, for example [21,23]. The resulting governing equations are similar to models based on the micromorphic approach, with an evolving length scale [24]. In fact, these approaches were introduced to avoid spurious spreading of the damage field (away from the nucleation point), which was observed when using a linear phase field potential in the former and a decreasing length scale in the latter.
In the present approach, based on the classical gradient based damage, the issue of compact support of damage field is avoided by using the length scale contrast to restrict the damage field to the interphase band. In order to recover physically meaningful results, which are objective with respect to the thickness of the interphase band, the interphase band needs to damage uniformly, even for bands that are not extremely thin. From this perspective, the classical gradient damage model provides better regularisation characteristics than these phase field approaches. If a phase field approach would be considered instead, models based on quadratic or double-well potential would suit best for the good regularisation characteristics needed in our application.
The numerical solution of the boundary value problem requires solving Eqs. (3) and (16) simultaneously. Note that these equations are coupled through the constitutive Eq. (5). Both differential equations are discretised using the FFT scheme. The coupled system of equations is solved using a staggered iterative scheme as detailed in an FEM context in [21]. The discretisation and residual evaluation of the mechanical equilibrium differential equation is well documented in [5]. Here, we only outline the residual formulation of the damage Eq. (16).
In order to avoid taking the Fourier transform of the length scale field, the non-linear residual is constructed by utilizing the fixed-point concept. We start by splitting λ 2 (x) into a homogeneous and a fluctuation field as Substituting Eq. (19) in Eq. (16) and taking the Fourier transform (F ) we obtain which can also be written as, Now using the property of the Fourier transform applied to differential operators, i.e, F ∂ n f ∂ x n = (ι k) nf , the above equation can be rewritten aŝ After some rearrangements and taking the inverse Fourier transform (F −1 ), Eq. (22) is rewritten in residual form as Initialising φ nl ≡ 1, Eq. (23) is solved using the matrix free Newton solver available in the PETSc library [25]. The frequencies used in the discrete Fourier transforms in this work correspond to the modifications presented in [26].

Integral based nonlocal damage
For the integral based nonlocal damage case, a weighted moving average is used to calculate the nonlocal field φ nl from the local damage φ l , i.e. [17,27].
where y is the position of the infinitesimal volume dΩ ∈ Ω i . The weighting function Ψ (y; x) that is commonly used is the Gaussian distribution, which reads where x is the position vector of the point at which the distribution is centred and λ is the characteristic length, which determines the distance along which Ψ decays to zero. The denominator in Eq. (24) normalizes the weighting function, which ensures that for a homogeneous φ l (x), the nonlocal field φ nl calculated using Eq. (24) coincides with its local counterpart. In our numerical setting, Eq. (24) is approximated by a cell averaging method, which can be written as where i and j represent material points in the interphase sub-domain and Δ x is the grid spacing. Equation (26) can be implemented by looping over all the points in Ω i . A more efficient way is to store the normalised weights in Eq. (26) in a matrix and implement the convolution via a matrix-vector product.

Results and discussion
In this section, we present and discuss the results of two numerical studies. The first is on the effect of the length scale contrast method for the gradient damage model. Next, we compare this method with the integral averaging model. The comparison is made based on their performance in providing a delocalised damage field and thus giving a numerically calculated fracture energy G n c in close agreement with two ideal solutions: the rate-independent analytical solution derived above (denoted by 'Rate indep.') and a ratedependent uniform damage strain solution denoted by 'Rate dep.'. The latter is easy to obtain in our 1D numerical setting: it only requires running simulations of the underlying

Model parameters used
The numerical studies are performed on the 1D bar of Fig. 2, discretised by 1000 uniformly spaced Fourier grid points. Young's modulus of the material is E = 100 GPa and a critical stress σ c = 1 MPa is used in all the simulations performed. The elastic properties throughout the bar are assumed to be same. The bar is loaded at an average strain ratė ε = 10 −5 s −1 for 5 s in time steps of Δt = 0.01 s. Since the considered properties are uniform, we introduce a small imperfection by means of a 1% reduction of the critical stress σ c for the grid points in the interval x ∈ [−l/20, l/20] within the interphase band, in order to trigger localisation.

Effect of length scale contrast in the gradient based model
In this part, we assess the performance of enforcing the interface condition (Eq. (17)) by the length scale contrast method. For this study, apart from the parameters described above, an interphase band thickness of l = 0.05L, reference strain rateε d0 = 0.1 s −1 and the final damage strain ε f = 0.001 were used. The nonlocal length scale used inside the band was λ in = l while outside it, the value λ out were varied to investigate the influence of the contrast (λ in /λ out ). Figure 3 depicts the distribution of nonlocal damage φ nl inside the interphase band and the region outside, to which it has spread significantly. Due to symmetry, only half of the band (and its surroundings) is shown. The profiles are plotted for the cases (λ in /λ out ) 2 = 10 1 , 10 2 , 10 3 , 10 4 . The nonlocal damage variable field extends more than 1.5 times the interphase band thickness into the bulk region on either side of the interphase for (λ in /λ out ) 2 = 10 1 . This is due to the poor approximation of Eq. (17) by Eq. (18) with this value of length scale contrast. However, this extension of the interphase damage into the bulk drops systematically for higher contrasts. The maximum contrast that the numerical solver for the damage equation could handle was (λ in /λ out ) 2 = 10 4 . For this contrast, we observe that the zone of spread of interphase damage into the bulk has reduced to less than 0.2 l. Figure 3 also shows that for high contrast ratios, the non-uniformity of φ nl caused by the imperfection decreases and the damage profile becomes increasingly uniform inside the band. Consistent with Eq. (18), the slope of φ nl at the two edges of the interphase tends to zero as the length scale contrast is increased. For the lower ratios, the interaction with the undamaged region next to the band results in a slower evolution of φ nl near the edge of the band compared to its centre.
The effect of this slowing down is observable in the averaged stress-strain response, see Fig. 4. As a result of it, the initial softening slope is less steep. This in turn results in an over-prediction of the amount of dissipation during the damage process as compared with that for the cases with a sufficient contrast, i.e. (λ in /λ out ) 2 = 10 3 -10 4 . Figure 5 (left) shows the calculated fracture energy G n c relative to the one predicted by Eq. (14). Remember that G c was derived for the rate-independent limit. Hence, we should expect the rate-dependent, computed G n c to approach G c , but not necessarily to asymptote to it. A much closer approximation of the rate dependent ideal solution (dashed curve) can be expected, as it captures the additional dissipation due to rate effects. Indeed, the ratio G n c /G c drops to less than 1.1 for the highest length scale contrast applied. This higher accuracy however entails significantly more iterations (Fig. 5 (right)) required to solve Eq. (23): this number goes up by a factor of 10 as the contrast is increased from (λ in /λ out ) 2 = 10 1 to 10 4 . This is essentially due to the

Comparison of gradient and integral nonlocal damage models
The analysis in "Analytical solution for the rate-independent case" section showed that for a uniform damage strain distribution in the rate-independent limit, choosing ε f according to Eq. (15), a constant fracture energy can be obtained, irrespective of the interphase band thickness l. In this section, we want to test the same hypothesis but in a numerical setting, where localisation is triggered by the presence of an imperfection (as discussed above). This is a test of the effectiveness of the two types of nonlocality in regularizing the damage (and damage strain rates) and thus facilitating the use of simple scaling Eq. (15) for the interphase band. Different interphase band thicknesses (l/L = 0.010, 0.015, 0.020, 0.025, 0.030) are studied.ε d0 = 1.0 s −1 and ε f = 5 × 10 −5 L/l in combination with the parameters in "Model parameters used" section are used. For both models, three different values of the nonlocal length scale are used: λ/l = 0.5, 1.0, and 1.5, where for the gradient damage based approach λ in = λ and λ out is varied such that a constant length scale contrast (λ in /λ out ) 2 = 10 4 is obtained. Figures 6 and 7 show the damage profiles obtained towards the end of the damage process for the gradient and integral model, respectively. The damage distribution obtained for the gradient damage model, shown in Fig. 6, is increasingly uniform inside the interphase band as the length scale λ is increased, for all the band thicknesses l. However, a slight gradient remains even for λ/l = 1.5. This is due to the imperfect insulation of the interphase damage from the surrounding (undamaged) bulk, as discussed in "Effect of length scale contrast in the gradient based model" section. The integral nonlocal model, in Fig. 7, shows a nearly constant damage inside the band for λ/l ≥ 1.0-again for all l/L. For the smallest length scale considered, λ/l = 0.5, some of the non-uniformity introduced by the imperfection remains.
A more discriminating result is the damage strain rateε d , which is shown in Figs. 8 and 9-again towards the end of the damage process and for the gradient and integral approach respectively. Here we do observe a significant amount of non-uniformity in all cases, demonstrating that the damage strain does localise late in the damage evolution process, triggered by the imperfection. The effect is however significantly milder for the integral model ( Fig. 9) compared with the gradient model (Fig. 8). Particularly for the integral  Fig. 10 The overall stress-strain response of the gradient damage model for different interphase band thicknesses and nonlocal length scales the applied average strainε for all cases considered. Also shown (in black) is the underlying rate-independent response as given (analytically) by Eq. (12). For both types of nonlocality, and for all l/L shown, the computed stress-average strain responses for l/L = 1.0 and 1.5 practically coincide. They furthermore deviate only slightly from the analytical solution. This shows that for these values of the nonlocal length scale both regularisation approaches are effective in rendering the global response objective with respect to the interphase band thickness, despite the non-uniformity observed in φ nl (Figs. 6, 7) and, particularly, the damage strain rate (Figs. 8, 9). The smallest value of λ, λ = 0.5 l, clearly is insufficient, as the softening responses computed for it are systematically steeper than for larger values (and than the analytical solution).
For a more quantitative assessment of the objectivity of the computed response with respect to the interphase band thickness adopted, the variation in the calculated fracture energy G n c , normalised by the ideal input fracture energy G c , is shown for the gradient damage method and the integral averaging method in Fig. 12. In the rate-independent limit, one would like to observe this ratio to be constant at unity. Given the fact that the numerical model is rate-dependent, however, a slight deviation should be anticipatedwhich ideally would be constant. To illustrate the effect of the rate dependence, numerical solutions for the rate dependent case with piecewise constant damage are included in both diagrams as black dashed curves. Note that these are slightly above G n c /G c = 1 and show a downward trend, which is due to the fact that post-peak the strain rate in the interphase band varies with the interphase band thickness. The fracture energy computed by the gradient damage approach show a discrepancy with both reference values (i.e. rate-independent and rate-dependent) on the order of 10%. For λ/l = 0.5, the fracture energy is under-estimated, whereas the larger values of λ result in an over-estimation. In each of these cases, the deviation furthermore varies with l/L. The integral nonlocal model performs much better. For λ/l = 0.5, which we already observed to be too small to obtain a uniform damage in Fig. 7, it under-estimates G c by up to 10% (Fig. 12). But for λ/l ≥ 1.0, it over-estimates G c by ≤ 5%, in a way which is consistent with the rate-dependent reference solution: the difference between the regularised solutions for λ/l = 1.0 and 1.5 and the piecewise constant reference solution is less than 1%. For practical purposes, this is more than satisfactory and the slight trend with l/L, due to the rate-sensitivity, should also pose no problem.

Conclusion and outlook
In this paper, the issue of nonlocality associated with a method to incorporate interface decohesion at polycrystalline interfaces approximated by an FFT based spectral method was discussed. Interfaces were approximated as interphase bands. The softening nature of decohesion required the use of nonlocality within the interphase domain. The applicability and performance of integral and gradient damage based nonlocal averaging methods have been discussed.
For the gradient damage approach, it was found that the mechanical response and dissipation depends on the accuracy with which the flux boundary condition at the edges of the interphase domain can be enforced. In order to restrict the damage to the interphase only, a flux free condition at the interphase boundaries was used. Since FFT solvers do not allow solving an equation on an irregular domain, in the gradient approach, a contrast in the nonlocal length scale was used to approximate the flux free boundary condition at the interphase edges. The error in the fracture energy reduces upon increasing this length scale contrast. However, a high contrast renders the gradient based approach highly heterogeneous and entails numerical problems.
The implementation of the integral approach on the other hand was straightforward and effective. It only requires storage of the nonlocal weights in a matrix that can be optimized using sparse storage. Regularisation length scale values equal to or greater than the interphase band thickness were found to give accurate predictions for the fracture energy, largely independent of the (arbitrary) interphase band thickness l. The slight (≤ 1%) variation which remains is due to the fact that the scaling of ε f according to Eq. (15) does not take into account the strain rate sensitivity of the damage model. From the current study it is very clear that integral approach offers more advantage from the computational efficiency point of view. We expect the same advantages to carry over in multidimensional cases on periodic microstructures. Nevertheless, this still remains to be tested.
We wish to emphasize that, although it enables a rigorous and transparent comparison, the 1D problem considered here may not reveal all complexities that could be encountered in two or three dimensions. For instance, care should be taken that the nonlocality introduced to homogenise the damage across the interphase band thickness does not affect the propagation of damage along the band in an unrealistic manner. One possible way to avoid this is by introducing anisotropy in the nonlocal averaging. In the gradient damage approach this can be achieved by having a tensorial form for the nonlocal length scale, while for the integral approach an orientation dependent averaging kernel can be used. Modelling decohesion of polycrystalline interfaces will require a proper treatment of triple junctions to avoid the issue of non-unique interface normals. Furthermore, a method to couple interface damage with the bulk damage in voxel based models [28] can also be explored to model more complex crack patterns-kinking and branching of cracks into the bulk. These issues are currently under investigation and will comprise future works.

Authors' contributions
LS carried out most of the study and drafted the manuscript. LS and PS developed the methodology and implemented it. RHJP, FR and MGDG conceived of the study and also participated in its design and coordination in addition to critically reviewing the manuscript. All authors read and approved the final manuscript.