An FFT-based spectral solver for interface decohesion modelling using a gradient damage approach

This work presents a fast Fourier transform (FFT) based method that can be used to model interface decohesion. The inability of an FFT solver to deal with sharp interfaces discards the use of conventional cohesive zones to model the interfacial mechanical behaviour within this framework. This limitation is overcome by approximating sharp interfaces (e.g. grain/phase boundaries) with an interphase. Within the interphase, the background plastic constitutive behaviour is inherited from the respective adjacent grains. The anisotropic kinematics of the decohesion process is modelled using a damage deformation gradient that is constructed by mapping the opening strains (in normal and tangential modes) to the associated projection tensors. The degradation (damage) of the interfacial opening resistances is modelled via a dimensionless nonlocal damage variable that prevents localisation of damage within the interphase. This nonlocal variable results from the solution of a gradient damage based regularisation equation within the interphase subdomain. The damage field is constrained to the interphase region by applying a relatively large penalisation on the damage gradients inside the interphase. The extent of nonlocality ensures that the damage is largely uniform in the direction perpendicular to the interphase, thus making its thickness the theoretical lengthscale for dissipation. To achieve model predictions that are objective with respect to the interphase thickness, scaling relations of the model parameters are proposed. The numerical performance is shown for a uniform opening case and then for a propagating crack. Finally, an application to an artificial polycrystal is shown.


Introduction
The industrial requirement of high strength and tough materials necessitates a clear understanding of these properties in relation to the role of the underlying microstructure. Computational tools that link microstructures with their macroscopic response are of great importance for the multiscale modelling of materials [9,23,24]. During the last two decades, the fast Fourier transform (FFT)-based spectral method has emerged as a useful tool to serve this purpose. The method was originally introduced to study the homogenised response of composites [18]. Lebensohn et al. [14] extended the algorithm to incorporate viscoplasticity to study polycrystalline microstructures at small strains. An extension to finite strain crystal plasticity was presented in [7]. Several modified formulations [17,34, and references therein] and improved solution schemes [6,26,35, and references therein] have been developed to deal with the high stiffness contrast in volume elements (VEs). The method is specially suited to perform micromechanical simulations on experimentally acquired images. An FFT-based crystal plasticity solver was used by [32] to understand the shear banding in pure magnesium. Geus et al. and Tasan et al. [10,31] used it to explore the stress-strain partitioning in DP-steels.
Multifield problems take direct benefit of the computational efficiency offered by the FFT-based solvers. The matrix-free solution scheme in particular lends the advantage of scalability when the system size, i.e. the grid density increases. Recent works on coupling of the mechanical equilibrium problem with damage or thermal fields [4,29] or martensitic phase transformation fields [13] are examples of the increasing popularity of the FFT solver for micromechanical problems. A necessary extension of these FFT-based damage modelling is the interface decohesion. The challenges associated with it and more generally with modelling of any other interfacial phenomena, need to be addressed.
The modelling of crack propagation is generally done by using either the cohesive zone (CZ) technique or the continuum damage method. The former requires the use of interface elements to model the displacement discontinuity across the surfaces representative of the crack faces. The evolution of this discontinuity is governed by a tractionseparation law (TSL). For an overview of cohesive surface modelling, please see [19] and the references therein. CZ approaches are typically finite element (FE) based. It is not clear how such interface elements, providing displacement discontinuity, can be embedded in FFT solvers. The continuum damage approaches that approximate the stress relaxation across the separating surfaces by degradation of a material volume have been successful in FE context, for example [28], and also FFT context [4,29], without being specific to inter/intragranular fracture. The volumetric approximation allows for a conventional constitutive law for the softening behaviour. This approach also has the advantage that a complete three-dimensional state of stress is available. Stress triaxiality within the damage process zone can influence the crack evolution by in-plane stretching, crack splitting etc. Motivated by this phenomenon, [22] contributed towards extending the classical cohesive zone approach by incorporating in-plane strain components in the formulation. This also necessitated the definition of a finite thickness region to normalise the related jump components.
Motivated by the use of cohesive zones for grain boundary sliding and separation in nanocrystalline FCC metals [33], analysis on TSL parameters influencing polycrystalline interfacial crack paths [25], it is worthwhile to explore methods such that FFT solvers can also be used for similar studies.
This research contribution introduces a method to include traction separation like behaviour for the interfaces in FFT solvers. Though a strongly simplified setting is considered, this contribution discusses the possible shortcomings and ways to address them. The volumetric approach is exploited by a priori identifying the grid points near the non-resolvable (due to grid discretisation) interfaces as interphase, see Fig. 1.
Anisotropic separation kinematics can then be employed within these interphases to reflect the decohesion activity of arbitrarily oriented interfaces to loading. Aslan et al. [1] incorporated this anisotropy in a small strain setting to model crystallographic fracture in Nickel-based super alloys. Subsequent works [15,28] extended it to a finite strain setting. The former applied it to model the damage initiation at lath martensite boundaries using a viscous regularisation model while the latter modelled the crystallographic cleavage using a phase field damage formulation. This interface as interphase approximation introduces a lengthscale (the interphase thickness, l), the role of which needs to be understood, see Fig. 2. Within the interphase, the aim is to have multiple Fourier points, such that for arbitrarily oriented interface the discretization artefacts are minimized .
It is well known that the continuum damage method must involve some sort of nonlocality to restore the well-posedness and mesh objectivity of the governing equations in the localisation regime. Therefore, a strongly nonlocal regularisation technique needs to be used. In a broad sense, this regularisation can be achieved using two approaches-the gradient based differential equation approach [8,20] and an integral equation approach [3]. In this contribution, we shall explore the former.
The organisation of this paper is as follows. After this introduction motivating the work, the theory (Sect. 2) of the constitutive model (Sect. 2.1) with some parametric studies to understand the homogeneous response of the model (Sect. 2.2) and the field problem (Sect. 2.3) are discussed. The numerical implementation is outlined subsequently in Sect. 3. The theoretical scaling of the damage parameter in order to obtain an interphase thickness insensitive work of damage dissipation is discussed in a 1D setting (Sect. 4.1). The scaling is then numerical assessed on the 1D bar (Sect. 4.2), with particular emphasis on the influence of the interphase subdomain boundary condition and the defect on the damage dissipation energy. A similar study is then performed on a more complex problem, i.e. crack propagation (Sect. 4.3). Finally, in Sect. 4.4, a polycrystalline application of this model to a cluster of grains is presented to show the potential of the FFT-based modelling for interface decohesion, followed by conclusion in Sect. 5.

Theory
The governing formulation consisting of the constitutive and kinematics formulations (Sect. 2.1), its response (Sect. 2.2) and the field formulation (Sect. 2.3) of the model are presented in this section.

Constitutive formulation and kinematics
The constitutive formulation consist of the elastic, plastic and damage models and their corresponding contribution to the kinematics. Exploiting the anisotropic damage models within a crystal plasticity framework [28], the current work on interfaces is also modelled in DAMASK [24] is used. F is the total deformation gradient.

Elasticity
The elasticity is based on the classical Hooke's law. A Green-Lagrange strain measure E e in the unloaded (damage) configuration Ω d is defined in the same sense as [11], as where F e is the elastic deformation gradient and I is the identity tensor. Using the elastic stiffness tensor 4 C, the elastic constitutive relation, yielding the second Piola-Kirchhoff stress S I in damage configuration reads, Equivalent definitions can be described in the plastic configuration after applying appropriate pull back operations.

Plasticity
To account for the plastic anisotropy, a plastic slip variable γ α is defined for each slip system α. Following the rate dependent description, for example described in [21,23], the evolution of γ α is given aṡ whereγ 0 is the reference shear rate, m is the rate sensitivity, τ α is the resolved shear stress on the slip system with direction d α and normal n α , obtained by projecting the Mandel stress M p on the projection tensor d α ⊗ n α as τ α = M p : ( d α ⊗ n α ). The critical resolved shear stress g α is transient and evolves as a function of the plastic slip according tȯ where h 0 is a hardening constant and g 0 and g ∞ are the initial and saturation values of the critical resolved shear stress, respectively. h αβ is the hardening matrix whose diagonal entries, i.e. for α = β, contribute to self-hardening and the off-diagonal entries (α = β) provide cross-hardening. a is the hardening exponent. The scalar contributions from each of the slip systems contribute to the rate of the plastic deformation gradient F p througḣ where L p is defined as,

Damage
Interfaces are prone to damage due to the heterogeneity imparted by differences in crystalline orientations (grain boundaries) or different phases (phase boundaries). Moreover, they are preferable sites for impurities to segregate and weaken them. On these accounts, upon loading, the interfaces may start to initiate damage by micro-void nucleation. During the loading processes these voids may grow, thus accommodating the applied deformation. This strain accommodation can be decomposed in a normal and a shear part. To represent them, scalar variables, namely the normal opening damage strain δ n and the shear damage strain δ s , are defined. The equations for their evolution are given as, for the normal mode anḋ for the shear mode. They are respectively driven by the normal (t n ) and shear (t s ) components of the resolved traction vector (M I · n I ) in the direction normal to the interface plane ( n I ) and parallel to the plane ( d I ), respectively. The vector d I is the normalised direction of the shear stress obtained after subtracting the normal component from the interface resolved traction vector. t n 0 and t s 0 are the critical stresses of the interfaces in normal and shear modes respectively. The term within the Macaulay brackets, (.) = ((.) + |(.)|)/2, represents the damage limit surface. In order to model the softening behaviour, the critical stresses are reduced by φ q nl as the individual opening strains (normal and shear) evolve. φ nl is a nonlocal damage field (or synonymously, a material integrity field) obtained from regularisation of φ l . A value of φ nl = 1 represents an undamaged material point which drops towards 0 as the opening strains evolve. The softening exponent q controls the softening slope.δ n 0 andδ s 0 are the rate constants scaling the < · > confined terms in Eqs. 8 and 9, and allow to model the rate dependent traction-separation behaviour, stemming from the underlying micro-mechanical activity like the micro-void growth, coalescence etc. The dependence of the local measure of damage on the normal and shear opening strain, calculated from Eqs. 8 and 9 respectively, is defined using a history variable δ H as where δ H i and δ H f are related to the opening strains (normal and shear combined) at which the stress-softening initiates (i.e φ l starts to drop from 1) and reaches complete failure, respectively. The evolution of the history variable at a material point is defined aṡ with the initial condition δ H = 0. The Macaulay brackets oṅ δ n enforce that φ l does not evolve when the normal opening strain decreases. The absolute value onδ s , implies that the formulation does not distinguish the role of positive and negative shear in sliding mode on degradation. It is noted that this model does not consider a full cyclic loading. The contribution of these openings to the interface damage part of the total deformation gradient, is obtained usinġ where L d consists of the normal and shear opening strain rates multiplied with their respective projection tensors as,

Response of the constitutive model
The response of the constitutive model is analysed next to understand how it works in conjunction with plasticity.
Depending on the orientation of the abutting crystals, their misorientation and the damage parameters, the overall softening response of a material point in interphase may be preceded by significant, little or no plasticity. Figure 3 depicts the different responses that the material points within the interphase region may give. Depending on the interaction between plasticity and damage, the total dissipation within the interphase can be calculated. Since the interphase zone inherits the plasticity of its surrounding crystals, the dissipation within this zone also occurs through plastic deformation. It is noted that the zone of plastic dissipation in not confined, unlike the damage model. In this work, we restrict our analysis to the damage dissipation. Another characteristic property of the model is its rate dependency, the extent of which at a straining material point, let's say, at a uniaxial tensile strain rateε =ε 11 can be controlled by the model parameterδ n 0 . Figure 4 shows the stress-strain response for different values of the reference opening strain rateδ n 0 , with E as the Young's modulus. For non-zero loading rates, smaller values ofδ n 0 , entail a (viscous) overstress surpassing the critical stress defined in the constitutive model. The peak stress always overshoots the rate-independent solution, which is approximated for larger values ofδ n 0 in Eq. 8. At final failure, all the responses converge. Hence, the model can sufficiently approximate the rate independent limit, and also highly rate dependent material behaviour, if needed. Similar discussion also applies to the shear opening strain evolution described in Eq. 9. The other important damage constitutive model parameter is the degradation function, φ q nl . Figure 5 shows that the exponent q controls the initial and final softening slope. Higher values give a steep softening at initiation and a gradual slope towards final failure. In the context of ductile damage, the degradation function controls the microscopic void growth and coalescence.

Field problem
At the boundary value problem level, the solution procedure consists in finding two primary unknown fields, the deformation gradient field F and the nonlocal damage field φ nl , in the reference configuration Ω. F(X) is retrieved as the solution of the quasi-static mechanical equilibrium equation where P is the first Piola-Krichhoff stress tensor, which is a nonlinear function of the field variables-F and φ nl -that are not known a priori. The FFT-based spectral implementation enforces periodic boundary conditions on the solution variable, i.e. the field of F satisfies, F + = F − , across the domain boundaries ∂Ω approaching from either side (represented here by + and −). The tractions across the boundaries are anti-periodic, implying, P · N | − = −P · N | + , where N represent the normals to the domain boundaries.  The second unknown field, the nonlocal damage φ nl , is obtained from the following regularisation equation where, φ l is the local damage variable as defined in Eq. 10 while φ nl is its regularized counterpart. The extent of regularisation is governed by the regularisation tensor D. The fields φ nl and φ l concern interface damage and hence ideally should only exist in sub-domain Ω I . Accordingly, the Eq. 15 should also be solved in Ω I only, with the flux-free condition at the edges ∂Ω I of the interphase with normal n I . However, in FFT-based spectral solvers, the interpolatary shape functions encompass the entire domain Ω. Hence, a numerical approximation of Eq. 15 can not be solved only in the subdomain Ω I . Under this restriction, a method to approximate the boundary condition Eq. 16 is proposed. This is done by adopting a non-uniform regularisation tensor D, defined as with a piecewise constant λ, having a value λ in inside the interphase Ω I and a value λ out → 0, outside it i.e. Ω \ Ω I . Accordingly, defining D in and D out using Eq. 17, the boundary condition at the edges of the interphase become which in the limit λ in /λ out → ∞, approximates Eq. 16. Equations 14 and 15 are coupled nonlinearly through Eqs. 8, 9, 13, 12 and 1 in a non linear manner.

Numerical implementation
The numerical implementation of Eq. 14 is discussed in [26], while in this work the implementation of Eq. 15 using an FFT-based spectral method is discussed. Similar to the idea of introducing a homogeneous reference media for the mechanical problem, an average diffusion tensorD is introduced. Splitting the regularisation tensor D into its homogeneous (average) and inhomogeneous part as, allows transforming Eq. 15 into, The above equation contains differential operators that are nonlocal. The Fourier transform provides an easy way to calculate them, as the nonlocal differential operations become local in Fourier space. Defining the Fourier transform and its inverse as F and F −1 , they operate on a function f as F f =f and its Fourier transform as F −1f = f , respectively. Defining the property, F f = 2πι ζf = ι kf , where ζ represents the frequency vector, f the gradient of f and ι as pure unitary imaginary number. Applying the Fourier transform on Eq. 20, is obtained. The inhomogeneity ofD makes the Fourier transform of the second term on the left hand side more involved than that of the remaining terms. It is simplified by first taking the derivative in Fourier space and then transforming the result back to real space to point-wise multiply the obtained vector field with the inhomogeneous field,D. Rearranging the terms and taking the inverse Fourier transform, the fixed point form, is obtained. A solution t n φ nl at time t n is required. Starting with an initial guess at the first iteration t n φ i=0 nl = t n−1 φ nl the solution t n φ nl can be found by solving the roots of the resulting residual, using a Jacobian-free Newton-Krylov scheme [12] as implemented in the PETSc library [2]. The solution of the coupled system, Eq. 23 and the discretised form of Eq. 14, is obtained using a staggered scheme, similar to [27,28].

Results and discussion
This section assesses the proposed interphase based interface decohesion model. The energy dissipated by the damage model, represented by G is used to examine the model behaviour.

A theoretical scaling of work of dissipation
A theoretical relationship between the damage dissipation energy and the interphase thickness is sought. A bar of length L with a sharp interface at X = 0 is considered in Fig. 6, loaded in pure tension monotonically. Considering this onedimensional case, the resolved normal opening stress t n is also the overall stress σ . This allows rewriting Eq. 8 as In order to compare Eq. 24 to a cohesive zone like tractionseparation law, a link to displacements needs to be made.
Integrating the opening strain distribution δ n and its ratė δ n along the thickness direction of the interphase X = [−l/2, l/2] gives a normal displacement separation Δ n and its rateΔ n which only approach the total displacement and displacement rate, respectively at the completion of softening. Depending upon the material properties or defects, the distribution of both δ n andδ n can be heterogeneous along X . Δ n andΔ n can be calculated by integrating them through the thickness, respectively. Motivated by the fracture energy trick [5], where damage element properties were scaled, a limiting case, such that both δ n andδ n are uniformly distributed (i.e. no localisation, φ nl = φ l ) is considered. Considering the application to the metallic systems and thus neglecting elastic deformation, within the interphase, Eq. 24 can be written as For the cohesive response given by Eq. 25 to be invariant of the interphase thickness, the conditions anḋ δ n 0 l =Δ n 0 (28) must be satisfied. Δ n i and Δ n f are defined as initial and final critical opening displacement at which damage initiates and final failure is attained.Δ n 0 can be regarded as the displacement rate constant. A similar derivation can be made for the shear traction-separation law. At final failure all the energy stored in the bar will be dissipated through the interphase. Thus, the area under the traction separation law obtained by scaling the material parameter as given in Eqs. 26, 27 and 28 altogether, should be independent of interphase thickness l.

Numerical analysis: 1D
Assuming isotropic elasticity and uniform bulk properties i.e., E 1 = E 2 = 10 11 Pa and ν 1 = ν 2 = 0, numerical experiments are performed on the 1D bar (Fig. 6), by applying a pure Two different cases, where the bar is discretised with a Fourier grid of 1000×1×1 and 2000×1×1, are considered to evaluate the effect of the discretisation. Within the interphase, the regularisation lengthscale λ in = l is applied, while λ out is set such that λ out = 0.01λ in is maintained. This ratio is maintained throughout the analysis. In an FFT based solution for problems involving heterogeneity of model parameters, the numerical convergence issues are expected. [30] provides a systematic analysis of the influence of this contrast of the mechanical response. Following the combined static and kinetic scaling, i.e. Eqs. 26, 27 and 28, Δ n i /L = 0, Δ n f /L = 5 × 10 −5 andΔ n 0 /L = 1s −1 are considered. The result of this scaling in terms of overall tensile response, true stress -true strain component σ 11 − ε 11 , is shown in Fig. 7. Ideally, the responses should match perfectly but some deviations can be clearly observed. These deviations decrease for the finer discretisation, Fig. 7b. Better convergence with a larger interphase thickness is also observed. Thinner interphases reveal somewhat less softening due to the larger influence of damage outside the interphase. This influence becomes smaller for finer discretisation, Fig. 8b. The curve l/L = 0.01 for both resolutions also reveals a less smooth softening. This could be attributed to the large value ofδ n 0 , leading to rapid growth of opening strain rates, as depicted in Fig. 10. Although the damage field seems to be uniformly regularised, the opening strain δ n in Fig. 9 and opening strain ratesδ n in Fig. 10, are still localised. There is a considerable gradient in the interphase region for local field variables.

Fig. 7 Overall response for interphases with no defect
The fracture energy is then calculated by integrating σδ n dt over the time of deformation depicted in Fig. 7. Figure 11 shows the plot for different l/L combinations. The ideal scaling would be a horizontal line, which is only approximated in the numerical simulations. An improved scaling behaviour is observed at higher resolutions.
In the simulations presented in preceding discussion, the properties within the interphase region were uniform. The heterogeneity in the resulting damage and strain fields is arising due to the approximation made in the boundary condition Eq. 18 at the interphase edges. For more realistic interphases, the constituting crystals may have different mechanical properties. This may induce asymmetry in the resulting damage and strain fields within the interphase. In order to examine this, the interphase simulations with l/L = 0.04 are reconsidered. A small defect spanning a thickness 0.2l from the left edge is placed. In this defect zone the critical stress is slightly reduced such that the critical stress becomes 0.995t n 0 . The overall response and the  Fig. 12. Different regularisations λ in = λ = {1.0l, 1.5l, 2.0l} are applied. Localisation is clearly visible in these results, which reduces for larger values of λ in , Fig. 12b. Although the damage fields are regularised significantly, the opening strain rates Fig. 12c, are localised. The effect of this is also visible in the overall response Fig. 12a and accordingly the influence on the fracture energy can also be expected. To reduce this influence of defect, stronger nonlocal regularisation can be explored.

Numerical analysis: 2D crack propagation
The response of the cohesive parameter scaling in a onedimensional setting is now extended to a more complex two-dimensional setting. Instead of examining damage evolution in a region with predefined dimensions approximating a sharp crack surface, analysis is now performed for propagating crack surfaces. The moving crack surfaces are approximated by a volume with smooth moving boundaries in the propagation direction, parallel to interface. Such a process zone is also associated with the classical cohesive zone model, where its size depends on the fracture parameters. For a given loading, fixed specimen dimensions and fracture properties, the overall peak stress can be determined. There is no ambiguity on the crack energy release rate in the propagation regime since the crack is modelled as physically interpreted, i.e., as a surface. In the nonlocal continuum damage based approaches, the process zone size is controlled by the additional regularisation parameter λ which affects the dissipation. Introducing a damage dissipative mechanism in   an interphase region with predefined finite thickness l, necessitates a clearer understanding of the relationship between l and G. As the interphase thickness dimension is purely numerically motivated, the model's response in the propagation regime should be insensitive to it.
For this purpose, a volume element of size L X × L Y resolved by a grid of 111 × 101 points as depicted in Fig. 13 is used. Since nucleation of a heterogeneous crack is required before the damage dissipation energy release rate is calculated, a stress concentrator is introduced in the left part of the specimen. Due to the fact that the FFT method can only deal with a regular grid and overall domains, this stress concentrator is modelled by a relatively very compliant phase. This phase also provides shielding of the periodicity related stress fields across X = 0, L X . The interface emanating from the notch is approximated by an interphase of thickness l. Different interphase thickness values, i.e. l = {0.0297L, 0.0485L, 0.0693L} are considered for the same geometry. The model parameters for the simulation are given in Table 1. For simplicity, the material is initially elastic and damage is the only available dissipative mechanism allowed to be active. A tensile loading rate such thaṫ F 11 = 10 −3 s −1 is applied upto an overall true strain of ε 22 = 0.818σ c /c 11 in vertical direction.
The overall response in terms of loading direction component of true stress and true strain (σ 22 −ε 22 ), is depicted in Fig. 14.
For different interphase thicknesses, the notch geometry remains the same, thus the stress prior to damage initiation is not affected. Since the opening criterion is stress based, the damage evolution (indicated by the start of deviation from the elastic response in Fig. 14) at the notch tip starts at the  Fig. 14 The overall stress-strain response for different interphase thicknesses same moment. For a thicker interphase, peak load is somewhat higher. This observation is in contradiction with the 1D result, where thinner interphases showed more stress owing to a larger rate activity. This can be attributed to the nonlocal effect, i.e. the averaging of stresses from a bigger zone with comparatively lower stresses. For a larger process zone, the critical damage is attained at a larger overall deformation. After elastic loading all three cases show a quasi-brittle softening branch with a steep slope that starts at similar overall deformation levels. The slope becomes steeper for thicker interphases. This can be attributed to lower opening strain rates thus lower overstress. On the same account, a larger strain to failure is also observed for thinner interphases. The peak stress values 0.56t n 0 − 0.59t n 0 at an overall strain value of approximatelyε 22 = 0.7σ c /c 11 . This insensitivity of peak stress is a good measure to assess the interface as interphase approach in comparison to cohesive zones where the interfaces are sharp. Simulations without scaling of the kinetic contribution (Eq. 28), show lack of correlation in terms of the final failure strain. Focus is next put on the incremental damage dissipation as the crack grows, which allows to quantitatively scrutinize the scaling response of the model considered and also the nonlocal continuum damage approach for fracture in general. The damage dissipation energy release as the function of the incremental crack growth Δa, requires a definition for a as the crack length. In the present model, where the process zone is significant, a unique definition of the crack length does not exist. A crack zone consists of material points which have lost their integrity. Similar to [16] where an equivalence between damage and sharp fracture phenomenon was sought, a definition is therefore proposed for the crack length a as where l Z is the size of the interphase in the out-of-plane direction. Figure 15 shows the evolution of the crack length a as a function of time. Very small values of a can already be identified before the peak load, in particular for the thinnest interphase which grows more steadily. The thicker interphases eventually start to grow faster and traverse the entire length of the interphase much quicker. Having defined a and moving on to the quantity of interest, the normalised (by static component G 0 ) damage dissipation energy release rate G with respect to crack length a, is plotted as a function of the crack length in Fig. 16. Ideally, the curves for all three cases should be identical in a fully developed crack regime. This is not exactly the case but still a remarkable similarity for thicker cracks is observed. In Fig. 16a, the thinnest interphase crack shows only a slightly higher amount of damage dissipation per unit crack growth. The results accounting for the scaling of the kinetic part also (Fig. 16b), indicate that (b) (a) Fig. 15 Evolution of the crack length with overall strain this deviation is largely due to the viscous contribution to the cohesive law given in Eq. 25. Moreover, the curves in Fig. 16b seem to saturate and cluster around a value that is related to the inherent fracture property G imparted by model parameters.

Polycrystal example
The one-dimensional study in Sect. 4.2 showed that the accuracy in the calculated work of dissipation improves as the resolution increases. In spite of the matrix-free advantage provided by the FFT method, the regular and equispaced grid is a matter of concern in obtaining high resolution results in regions of interest. Using any other grid entails loss of efficiency. Under these restrictions and having quantified the  (Table 2) and ferrite phase crystal plasticity parameters (Table 3) such that low angle grain boundaries are obtained. The compliant phase follows isotropic elastic response with the elastic constants (c 11 = 10 8 Pa and c 12 = 100 Pa). Due to this interruption of periodicity, effectively the grain count equal to six should be considered. The interphase approximating the interface has a thickness l = 0.081L X . A predominantly horizontal loading (in terms of conventionally used averaged quantities, the  Table 2 Orientations of the grains of the volume element shown in Fig. 17 in terms of their Bunge-Euler triplets (ϕ 1 , ϕ, ϕ 2 ) The overall response (Fig. 18), depicts an elasto-plastic response followed by softening. As the crack propagates, stress decrease is visible in Fig. 18 with comparatively slower drops corresponding to when the crack traverses more inclined (with respect to the loading) interfaces. The sharp stress drop is a result of brittle parameters chosen for this simulation. The damage field is shown in Fig. 19. Due to large resolved stresses, crack nucleates along an interface between grain 1 and grain 4. The damage on this interface (between grain 1 and grain 4) is predominantly in normal opening mode. The crack then propagates on to a less inclined interface (between grain 1 and grain 2) wherein the damage evolves in both modes. Similar behaviour is observed upon further propagation along interfaces between grain 3 and grain 2, and grain 3 and grain 1. This example demonstrates the capability of the model to simulate crack branching at the triple junctions without any additional criterion with interface damage confined significantly within the interphase region (Table 4).

Conclusions
In this contribution, a method has been proposed that enables FFT-based simulations of interface decohesion using a volume-based material degradation model. As a localisation limiter, a gradient damage approach within the interfacial (interphase) subdomain was adopted. An approximate method to solve the gradient damage differential equation within this subdomain was presented. The effectiveness of the approach is tested in a onedimensional setting with an evolving damage variable in an interphase. The significantly confined damage activity within the interphase allows applicability of a theoretical approach to damage model parameter scaling. The scaling allows to exploit the interphase thickness as a purely numerical parameter.
The damage dissipation was found not to be very sensitive to the thickness parameter even in a more complex setting of straight propagating crack. Interestingly, the peak stress in the overall response also showed a similar insensitivity. The response of the model is therefore comparable with the more established FE-based methodology for interfaces, i.e. the cohesive zones. The application of the model to a grain cluster with inclined interfaces, shows that the model can work in synergy with other mechanisms like plasticity.
Some challenges still remain. For example, strong inhomogeneities within the interphase may require a larger regularisation. This would further assess the effectiveness of lengthscale contrast trick, i.e. λ in >> λ out such that λ out → 0. The requirements to confine the interfacial damage inside the interphase and the related scaling, may need a large contrast of the lengthscales in the gradient damage equation, thereby making its numerical solution more difficult even with flexible iterative schemes. Other regularisation schemes may therefore be explored and compared as well.
is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.