Extending intergranular normal-stress distributions using symmetries of linear-elastic polycrystalline materials

Intergranular normal stresses (INS) are critical in the initiation and evolution of grain boundary damage in polycrystalline materials. To model the effects of such microstructural damage on a macroscopic scale, knowledge of INS is usually required statistically at each representative volume element subjected to various loading conditions. However, calculating INS distributions for different stress states can be cumbersome and time-consuming. This study proposes a new method to extend the existing INS distributions to arbitrary loading conditions using the symmetries of a polycrystalline material composed of randomly oriented linear-elastic grains with arbitrary lattice symmetry. The method relies on a fact that INS distributions can be accurately reproduced from the first (typically) ten statistical moments, which depend trivially on just three stress invariants and a few material invariants due to assumed isotropy and material linearity of the polycrystalline model. While these material invariants are complex averages, they can be extracted numerically from a few existing INS distributions and tabulated for later use. Practically, only three such INS distributions at properly selected loadings are required to provide all relevant material invariants for the first 11 statistical moments, which can then be used to reconstruct the INS distribution for arbitrary loading conditions. The proposed approach is demonstrated to be accurate and feasible for an arbitrarily selected linear-elastic material under various loading conditions.


I. INTRODUCTION
The ability to predict damage initiation and progression in structural materials heavily relies on understanding the local mechanical stresses present in the material.For certain aging processes, damage to the material begins at grain boundaries (GBs) where intergranular microcracks can form.Over time, these microcracks may grow along the GBs and merge into larger macroscopic cracks, ultimately threatening the structural integrity of the entire component under load.Therefore, it is essential to develop accurate models that can predict GB damage initiation based on sufficient knowledge of stresses present at GBs.An example of material aging mechanism that affects the GBs is an intergranular brittle fracture, which can occur in structural steels if the thermal history of a component leads to the segregation of phosphorus at the GBs [1].
Various models have been developed to address damage initiation and propagation on different length and time scales.Among these, probabilistic modeling of GB damage initiation has emerged as a valuable tool for predicting and mitigating damage in components under specific loading conditions, and for identifying critical locations where damage is most likely to occur.This method involves quantifying the probability of GB damage initiation as a function of local (GB normal) stresses, material properties and microstructure (grain morphology and crystallographic texture), and takes into account the variability and statistical distribution of these properties.One important input parameter for this type of modeling is the intergranular normal stress (INS) distribution, * samir.elshawish@ijs.siwhich needs to be estimated locally (i.e., on a scale of several hundreds of grains) for a given set of macroscopic loading conditions and material properties.
An example of a probabilistic modeling approach using INS distributions as an input parameter is the intergranular microstructure-informed brittle fracture (IG-MIBF) model, which is based on the original MIBF model [2] developed in the framework of the local approach to fracture defined by [3] and used since then to predict the probability for brittle fracture of low alloy steels [4].The MIBF model is used to capture the effects of some microstructure characteristics such as the scatter of critical stresses induced by the carbide size distribution and the incompatibility stresses arising from the polycrystalline structure under applied deformation.Although the original MIBF model considers maximum principal stresses, the IG-MIBF model accounts for INS distributions, which typically require numerical evaluation at different loadings (defined by, e.g., stress triaxiality) using techniques such as crystal plasticity finite element (FE) simulations [5][6][7][8][9][10] or crystal plasticity fast Fourier transform simulations [11][12][13].While these simulations can provide accurate results, they can become impractical and time-consuming when applied to numerous different loading cases.As a result, efforts have been made to identify the explicit role of external loading on INS distributions [14,15].Despite recent developments, only a limited number of models exists in literature for predicting INS distributions in polycrystalline materials under general loading conditions.For instance, in [15] we proposed a simple semi-analytical model for elastic polycrystalline materials with arbitrary-shaped and randomly-oriented grains under uniform loading.This model provides algebraic expressions for local INS and corresponding uncertainties, which explicitly include the loading stress components.It was derived in a perturbative manner and systematically validated against FE simulations of various Voronoi microstructures with cubic crystal lattices.However, the model has not yet been calibrated for lower (non-cubic) lattice symmetries.
The primary goal of this study is to develop a new alternative method for INS distribution modeling, based on the assumed properties of an isotropic linear-elastic polycrystalline material.Here, the isotropy of a polycrystal is understood in statistical sense, emerging from a sufficiently large number of anisotropic linear-elastic grains of arbitrary lattice symmetry, which are randomly oriented (providing zero crystallographic texture) and randomly shaped (providing zero morphological texture) 1 .The idea is to leverage these symmetries to establish an efficient and practical approach for extending the INS distributions from a limited number of loading conditions to any arbitrary uniform macroscopic stress state 2 .To achieve this, we derive a general expression for central statistical moments that explicitly incorporates the effect of loading stresses, enabling us to construct accurate INS distributions for various loading scenarios.
The paper is structured as follows.We begin in Sec.II by developing a general expression for the central mo-ments of the INS distribution, and then present the strategy for extending the existing INS distributions to arbitrary loading stress.In Sec.III, we demonstrate the feasibility of our method using the FE simulations.Finally, in Sec.IV, we provide some concluding remarks.All technical details are given in the Appendices.

A. Central moments
The purpose of this Section is to identify the role of the loading stress tensor Σ on the behavior of INS distribution and its respective statistical (central) moments in a general isotropic linear-elastic polycrystalline material subjected to Σ.This will allow us later to develop a method to extend the existing INS distributions (calculated at specific Σs) to arbitrary loading conditions.
Since a polycrystalline aggregate under loading stress Σ is assumed to be isotropic 3 , the INS distributions should be invariant to any aggregate rotation R. Equivalently, the INS distributions should be invariant to any loading stress rotation RΣR T .Selecting such R that RΣR T becomes diagonal, we may proceed with the most general loading stress where Σ i are loading-stress eigenvalues and I i unitloading tensile stresses along the i direction.If a local INS solution on an arbitrarily selected GB is σ nn for applied I i , then a local INS solution for applied Σ is following the linearity principle (of the assumed generalized (3D) Hooke's law).Since all three directions i are equivalent in an isotropic polycrystalline aggregate, the 1 In principle, grain shapes and crystal lattice orientations can exhibit a correlation, provided that the isotropic response of the polycrystal is preserved. 2A non-uniform loading can be approximated as a quasi-uniform stress state within a sufficiently small material domain.The applicability of this study can be extended to such a domain when its volume exceeds that of a representative volume element (e.g., several thousand grains). 3Although the grains are assumed anisotropic (with arbitrary crystal lattice symmetry), a polycrystalline aggregate becomes isotropic in the limit of a very large number of randomly shaped grains and randomly oriented crystal lattices.
average INS scales with the trace of loading (see also [15]), where ⟨. ..⟩ represents the averaging over all GBs or a subset of GBs of specific type4 in a polycrystalline aggregate and C ≡ σ (3) nn is a constant, which depends on elastic material parameters and correlation strength between grain shapes and grain lattice orientations [15].If the grains possess cubic lattice symmetry, then C = 1/3 (see Appendix A).Here we assumed that the grains are composed of the same material and they differ only by their crystallographic orientations.
The m-th central moment of the INS distribution is defined as for the applied stress Σ.It is convenient to decompose Σ into hydrostatic (Σ hyd ) and deviatoric parts (Σ dev ), Σ = Σ hyd + Σ dev , and analyze both contributions to µ m separately.Since Σ hyd is fully symmetric, its effect on µ m is significantly simplified.

Effect of Σ dev
We begin by deriving the central moments under the assumption of purely deviatoric loading.In fact, this loading is the sole contributor to µ m when the grains posses cubic lattice symmetry 5 .For non-cubic crystal lattices, however, the effect of hydrostatic loading is nonzero.This aspect is discussed later in Secs.II A 2 and II A 3.
The evaluation of the central moments can be reduced to when the applied stress Σ is replaced by the deviatoric stress Σ dev , which is traceless by definition, In this respect, only two (out of three) deviatoric stress invariants, J 2 ≡ tr(Σ 2 dev )/2 and J 3 ≡ det(Σ dev ), control the shape of the INS distribution.In the case of cubic crystal lattices, any hydrostatic contribution enclosed in the first stress invariant I 1 ≡ tr(Σ) = tr(Σ hyd ), manifests only as a shift of the INS distribution 6 .
Using Eqs.(2), ( 5) and ( 6), the m-th central moment can be expressed as Applying again the statistical equivalence of three spatial directions, the material average term (i.e., independent of loading) in Eq. ( 7) becomes invariant to the permutation of the three indices 1,2 and 3, Thus, a 3-parameter average M (i, j, k) (with i, j, k = 0, 1, 2, . . .and i + j + k = m) can be reduced to a 2parameter average (using a commutative property of the multiplication of three indices), This allows us to rewrite the m-th central moment in a symmetrized form where the averaging function M m depends solely on the (linear-elastic) properties of the grains for the assumed zero crystallographic and zero morphological texture of the aggregate.Moreover, if Σ1 , Σ2 and Σ3 in Eq. ( 10) are expressed in terms of the two (commonly chosen) deviatoric stress invariants J 2 and J 3 , the expressions for µ m dev can be further simplified, Here, a new parameter M i,j is introduced to collect all the material prefactors M m (k) in front of each J i 2 J j 3 loading term (with i, j = 0, 1, 2, . . .and 2i + 3j = m).In this way, the µ m dev finally simplifies to The first 20 central moments are presented in Tab.I. Finding a general expression for µ m dev under purely deviatoric stress loading (note, however, that µ m = µ m dev for cubic crystal lattices) by decoupling the loading (J 2 , J 3 ) and material (M ) contributions, is the first main result of this study.
It is relevant to note that the simplicity of µ m dev in Eq. ( 12) results from the following assumed symmetries: (i) the isotropy of the polycrystalline material, (ii) the linearity of the employed Hooke's law and (iii) the absence of hydrostatic loading (or its zero effect on the shape of INS distribution in the case of cubic lattice symmetry).Moreover, in the expression for µ m dev , the two deviatoric stress invariants J 2 , J 3 are independent parameters which are, however, bounded7 by J 2 ≥ 0 and 27. Also, since odd central moments µ m dev (m odd) scale with only odd powers of J 3 , the sign inversion Having a relatively simple expression for µ m dev for a general isotropic linear-elastic polycrystalline material subjected to purely deviatoric loading stress Σ dev (i.e., J 2 and J 3 ), allows to extract the unknown material parameters M i,j from N existing INS distributions obtained at few specific Σ dev .For example, if two (N = 2) such INS distributions are available (e.g., from FE simulations or literature), along with their central moments, one can identify all M i,j of the first 11 (and also 13-th) central moments (see Tab. I), assuming that J 2 ̸ = 0 and J 3 ̸ = 0 of the N = 2 existing INS distributions 8 .With this, one can calculate the new central moments µ m dev for m ≤ 11 and arbitrary Σ dev (i.e., J 2 and J 3 ).As demonstrated in Sec.III, m ≤ 11 is usually sufficient to accurately (re)construct the INS distributions.

Effect of Σ hyd
Grains with lower (non-cubic) crystal lattice symmetries yield non-trivial local responses even under purely hydrostatic (symmetric) loading conditions, Σ = Σ hyd .This indicates that the corresponding central moments are non-zero, µ m ̸ = 0 Σ hyd .Consequently, in the derivation of µ m , it is essential to consider a complete loading, Σ = Σ hyd + Σ dev .
Following the same derivation steps as in Sec.II A 1, the final expression for µ m , defined in Eq. ( 4) for arbitrary loading Σ, simplifies to (m > 1) where all three stress invariants, I 1 ≡ tr(Σ), I 2 ≡ tr(Σ) 2 − tr(Σ 2 ) /2 and I 3 ≡ det(Σ), multiply the unknown material prefactors M i,j,k .Since several combinations (i, j, k) satisfy i + 2j + 3k = m for a given m-th central moment, many more (N ≫ 1) existing INS distributions are required to extract M i,j,k even for relatively low orders m, For example, if four (N = 4) existing INS distributions (obtained at different loadings Σ with I i ̸ = 0) are available, it is possible to determine the corresponding M i,j,k for only the first four central moments.Clearly, the method outlined in Sec.II A 1 for extracting M i,j from µ m dev proves impractical here for any realistic application.
To circumvent this, we propose an alternative path.The core idea is that on any given GB the two INS contributions, σ nn = σ dev nn + σ hyd nn , denoted as σ dev nn for purely deviatoric loading Σ dev and σ hyd nn for purely hydrostatic loading Σ hyd , can be approximated as being independent of each other.This follows from the notion that σ dev nn is strongly dependent on GB normal orientation n, while σ hyd nn is not (or much less) since there is no preferred direction in a (random) aggregate under purely symmetric (hydrostatic) loading 9 .Under this approximation, σ hyd nn can be modeled as stochastic variable with assumed normal probability distribution 10 centered around the mean 9 The assumed independence of σ dev nn and σ hyd nn can be easily tested in the FE aggregate by calculating a linear correlation between the two values on each GB using Pearson correlation coefficient r.For various Σ dev and Σ hyd used in Sec.III B, we indeed found very low vales, |r| < 0.06, which supports the assumption. 10The absence of a preferred direction under purely hydrostatic loading, coupled with the concept of a random crystallographic texture and random morphological texture, implies that PDF(σ hyd nn ) should be a symmetric function around its mean.From this perspective, a normal probability distribution appears to be an appropriate choice, as also supported by numerical results in Fig. 3. value for isotropic (or cubic crystal lattice) grains I 1 /3, and with variance µ 2 hyd , Since the response of grains is linear with respect to (hydrostatic) loading, the variance scales with I 2 1 so that µ 2 hyd = I 2 1 M 2,0,0 , where M 2,0,0 is unknown material parameter 11 to be extracted from a single existing INS distribution obtained at hydrostatic loading Σ hyd = 1  3 I 1 1 3×3 .With this, all higher (even) central moments of normal distribution follow straightforwardly, 3. Effect of Σ hyd + Σ dev Since the (exact) approach outlined in Eq. ( 13) is not feasible, the final (approximate) expression for the central moment µ m is obtained from the probability theory, which states that the probability distribution of the sum of two or more independent random variables is the convolution of their individual distributions.In this respect, where µ m dev is the (exact) m-th central moment of INS distribution evaluated for purely deviatoric loading, Eq. ( 12), while µ m hyd is the (approximate) m-th central moment of INS distribution evaluated for purely hydrostatic loading, Eq. ( 16).In the limit of vanishing hydrostatic contribution (i.e., when grains posses cubic lattice symmetry or I 1 = 0), the expression in Eq. ( 17) correctly reduces to µ m = µ m dev .
To summarize, we outline the following methodology for computing central moments for arbitrary crystal lattice symmetry and external loading Σ, utilizing a limited number (N + 1) of existing INS distributions evaluated at specific Σs: 1. determine the unknown material parameters M i,j in µ m dev from N existing INS distributions obtained at N specific purely deviatoric loadings Σ dev (with I 1 = 0), 2. determine the unknown material parameter M 2,0,0 from a single existing INS distribution obtained at purely hydrostatic loading Σ hyd = 1 3 I 1 1 3×3 (with 3. calculate the central moments µ m from obtained prefactors M i,j and M 2,0,0 and arbitrary loading Σ (any I 1 , J 2 , J 3 ) using Eqs.( 12), ( 16), (17) and Tab.I.
The proposed method is exact (limited by the accuracy of the inputs) for isotropic linear-elastic polycrystals composed of grains that exhibit (i) cubic lattice symmetry under any external loading Σ, or (ii) non-cubic lattice symmetry under purely deviatoric external loading Σ dev .In both cases, µ m = µ m dev .As validated in Sec.III, the method performs very accurately also for most general crystal lattices and loadings.

B. Construction of INS distribution
The next step is to construct the INS distribution PDF(σ nn ) from first K central moments µ m (m ≤ K).There are various strategies how to get the spectra most representative for K → ∞ (see, e.g., [16]), expecting a smooth function PDF(σ nn ).We follow here the procedure proposed by [17] to obtain the following expression (for details see Appendix B) PDF(σ nn − ⟨σ nn ⟩ ; K, P, λ) ≈ − 1 π Im lim ϵ→0 P/P S(ξ;K,λ) ξ→ (σnn−iϵ)± √ where P/P S(ξ;K,λ) is a diagonal Padé approximant (of order P , typically P ∼ 6) to the modified moment series about the point ξ → ∞.The modified moments G m (λ), defined as are linear combinations of central moments µ i with prefactors λ m−i (with 0 ≤ i ≤ m) listed in Tab.III 12 .In addition to the two newly introduced (integer) parameters K and P in Eq. ( 18), the parameter λ is a real number, which represents approximately the lower and upper bounds of the (centered) INS distribution as In practice, λ should be as small as possible (while still obeying Eq. ( 21)), which is usually achieved by setting λ to a fraction of von Mises loading stress Σ mis = √ 3J 2 , e.g., λ/Σ mis ∼ 0.7.This scaling (λ ∼ Σ mis ) results from the fact that standard deviation of INS distribution scales as µ 2 ∼ √ J 2 ∼ Σ mis (see Tab. I).However, the optimal values for λ as well as P depend (also) on the number of available central moments K.For a given number N (and thus K) of available INS distributions, the fine-tuning of λ and P can be thus achieved by back-fitting of the N reconstructed INS distributions obtained at specific loadings.The optimized λ and P can then be used to predict the INS distributions for other loadings.
In Eq. ( 18), the off-shift by ⟨σ nn ⟩ in the argument of PDF is required because the central moments are used in G m (λ) instead of raw moments (see Appendix B).Also, the "±" sign should be understood as "+" for σ nn ≥ 0 and "−" for σ nn < 0, while a reasonably small value of ϵ (e.g., ϵ/Σ mis ∼ 0.001) should be taken instead of lim ϵ→0 to provide smooth distributions in the range of interest, Eq. ( 21).

III. VALIDATION OF THE METHOD
The proposed method is validated against polycrystalline FE simulations of an artificial quasi-isotropic aggregate model (4000 randomly shaped and oriented grains) under various (uniform) loading conditions Σ.More details about the FE simulations can be found in Appendix C. For demonstration purposes, two materials were selected for the grains: gamma iron (γ-Fe) with cubic crystal lattice symmetry (see Sec. III A) and calcium sulfate (CaSO 4 ) with orthorhombic crystal lattice symmetry (see Sec. III B).The corresponding material properties are listed in Tab.IV.

A. Cubic crystal lattice
Figure 1 demonstrates the feasibility of the method presented in Sec.II B for the reconstruction of PDF(σ nn ) from available first K central moments µ m = µ m dev (m ≤ K).The moments µ m were calculated numerically from the FE simulation of γ-Fe polycrystal under tensile loading (Σ = I 1 ) and then used to reconstruct the PDF by employing Eq. ( 18) with different parameter sets (K, P, λ).The reconstructed PDFs show (Fig. 1(a)) a rapid convergence with increasing K, reaching wellconverged results already for K ∼ 9. Also, the best performance of the method is obtained for P ∼ 6 (e.g., P = 5, 6 or 7 provide the same result, see Fig. 1(b)) and λ/Σ mis ∼ 0.7 (Fig. 1(c)).The derived optimal values, K ≥ 9, P ∼ 6 and λ/Σ mis ∼ 0.7 should not change significantly when using other material properties 13 .
In the following, we present the feasibility of the method introduced in Sec.II A for the calculation of central moments (and corresponding PDFs) at arbitrary loadings Σ given the central moments at few (N ∼ 1) specific Σs.Two cases are considered in Fig. 2, one with N = 1 and one with N = 2 available input PDFs.Having only one such distribution at hand, restricts us to identifying M i,j of the first K = 5 central moments (see Tab. I), which is obviously insufficient for accurate reconstruction and prediction (see Fig. 2(a)).However, the situation is greatly improved when two PDFs become available (see Fig. 2(b)).It this case, K = 11 central moments are clearly sufficient to correctly reproduce and extend the existing two INS distributions to arbitrary loading conditions (defined by J 2 = 1 and varying J 3 , see next paragraph).This holds true for any pair of input PDFs provided that their |J 3 |s are different and nonzero.As practically identical results were obtained for N = 3 (and K = 17), they were omitted from Fig. 2 for brevity. 13The most relevant material parameter here is (elastic) anisotropy of grains.More anisotropic materials produce wider INS distributions, thus larger λ/Σ mis is anticipated there.To exhaust different loading possibilities in a condensed manner, σ nn was off-shifted in Fig. 2 by ⟨σ nn ⟩ ∼ I 1 and normalized by von Mises stress Σ mis .Since the latter scales with J 2 , it is actually sufficient to fix J 2 to an arbitrary (positive) value and vary only J 3 .In Fig. 2, therefore, various possible PDF shapes (for a given material γ-Fe) were reconstructed by setting J 2 = 1 and varying |J 3 | ≤ 2/ √ 27 ∼ 0.385.This approach allows for the behavior of the PDF((σ nn − ⟨σ nn ⟩)/Σ mis ) to be controlled by a single dimensionless parameter J 3 /J 3/2 2 , as opposed to the PDF(σ nn ), which requires three parameters (I 1 , J 2 and J 3 ) to be specified.
To summarize, the second key result of this study, illustrated in Fig. 2, shows that for cubic crystal lattices we can accurately extend the INS distributions from specific to arbitrary loading conditions using only two input distributions evaluated at two specific loading conditions.

B. Orthorhombic crystal lattice
In this section, the methodology for INS distribution reconstruction is validated for orthorhombic crystal lattice system CaSO 4 at various loading conditions.The emphasis is put on hydrostatic stress, which (contrary to cubic crystal lattices) provides non-zero contribution to central moments µ m .As shown in Fig. 3, the effect of purely hydrostatic stress, Σ hyd = 1  3 I 1 1 3×3 , to INS distribution is indeed non-zero and can be accurately approximated 14 by a normal distribution centered at I 1 /3 and 14 Small deviations observed at higher stresses are attributed to finite size effects, which emerge due to insufficient aggregate size Applied external stresses Σn used in Fig. 4. In the labeling of input-PDF stresses, one additional letter is used ("d" for purely deviatoric, and "h" for purely hydrostatic stress).
with variance scaling as I 2 1 .This supports the modeling assumption made in Sec.II A 2.
The effect of general applied stress is analyzed in Fig. 4, where two approaches are compared.While in Fig. 4(a) two PDFs are considered as input (N = 2), both evaluated at purely deviatoric stresses (I 1 = 0), in Fig. 4(b) one additional input PDF (N = 1 + 2) is taken at purely hydrostatic stress (J 2 = J 3 = 0).In contrast to Fig. 2, the resulting INS distributions are presented in absolute scale to facilitate the effect of I 1 , J 2 and J 3 on PDF position and shape.
The results of Fig. 4 demonstrate that INS distributions can be accurately predicted for both deviatoric (I 1 = 0) and non-deviatoric applied stresses (I 1 ̸ = 0).Best accuracy is achieved by employing the strategy outlined in Sec.II A 3 for the estimation of central moments under both Σ hyd and Σ dev (see Fig. 4(b)).However, when the effect of Σ hyd is neglected (i.e., by omitting (4000 grains) and/or too coarse FE mesh (∼5 million elements).the Σh1 input), the agreement between FE simulations and predictions may reduce, especially at larger values |I 1 | (see top curves in Fig. 4(a)).This validates the proposed methodology also for non-cubic crystal lattices.Similar to the γ-Fe case shown in Fig. 2, also for CaSO 4 the K = 11 central moments are sufficient to correctly reproduce and extend the existing INS distributions to arbitrary loading conditions.As practically identical results were obtained for N = 1 + 3 (and K = 17), they were omitted from Fig. 4 for brevity.
In summary, for isotropic linear-elastic polycrystals with non-cubic crystal lattices three input INS distributions evaluated at three specific loading conditions are sufficient to accurately extend the INS distributions to arbitrary loads.

C. Limitations on the applicability of the method
The limitations of the proposed method, such as the assumption of linear elasticity of the grains and the absence of crystallographic and morphological texture in the aggregate, restrict its applicability to studying INS (distributions) in untextured aggregates under predominantly elastic loads.This could be particularly relevant for predicting the onset of GB damage in various ma-terials.Many failure mechanisms, which result in GB crack initiation, often begin in the elastic regime of the grains.Examples include intergranular stress corrosion cracking, hydrogen embrittlement, radiation-induced embrittlement, low-temperature brittle fracture, and intergranular corrosion.In all these cases, damage (and consequently cracking) can occur at GBs without significant plastic deformation.By examining GB normal stresses in this regime, it is possible to predict the onset of these failures, which is crucial for ensuring material reliability and safety.
Although the method requires both crystallographic and morphological textures to be absent for an isotropic polycrystal response, it does not limit the potential correlation between grain shape and the orientation of the crystal lattice.This allows the method to be used for studying INS (distributions) not only on random GBs but also on special GBs, where the two neighboring grains have fixed crystal lattice orientations relative to the GB plane.Such local arrangements of lattice orientations typically provide a certain strength to a GB, which can be compared to the corresponding INS for damage modeling.While many interesting results have been obtained on cubic crystal lattices [14,15], the responses on such GBs in less symmetric crystals are less known.In this regard, the proposed method could prove to be useful.

IV. CONCLUSION
In this study, we have derived a general expression for the m-th central moment of the intergranular-normalstress (INS) distribution when evaluated in an isotropic linear-elastic polycrystalline material under uniform applied stress.These central moments depend on three stress invariants and only a few (generally unknown) material invariants, due to material linearity and rotational symmetry of a polycrystal.Based on these findings, we have developed a new strategy to extend existing INS distributions from specific stress loadings to arbitrary loading conditions.The method relies on the fact that INS distributions can be accurately (re)constructed using only few first central moments (typically ten) whose material invariants can be extracted numerically from the existing INS distributions.Our study has demonstrated the accuracy and feasibility of this approach for two selected linear-elastic materials with different crystal lattice symmetries, showing that in a general case only three input INS distributions evaluated at three specific applied stresses are sufficient for accurately extending them to arbitrary loading conditions.These results have important implications for probabilistic modeling of grain boundary damage initiation on a component scale, particularly in the presence of complex stress states that can arise in front of propagating macroscopic discontinuities such as cracks.
the above integrand about z → ∞ and integrating over dω ′ , we get where ⟨ω m ⟩ is the m-th (raw) moment of F (ω).However, if we assume in Eq. (B3) that z = ω − iϵ (with both ω and ϵ > 0 real), then as ϵ → 0 we obtain (using Sokhotski-Plemelj theorem) A direct approximation of the moment series in Eq. (B6) by, for example, Padé approximants is doomed to failure as one obtains only a set of poles which represents the branch cut along the Re(z) axis.The imaginary piece of F (ω) is simply a set of δ functions rather than a continuous function which we wish to obtain.In order to circumvent this difficulty, [17] first carried out a nonlinear transformation where 2λ is the length of the branch cut on the Re(z) axis which is here assumed to be centered at ω = 0.In the complex ξ plane the branch cut is mapped onto a circle of radius λ.The interior of the circle contains the unphysical sheet of F (ω) to which any spurious singularities are confined.The exterior of the circle is the physical sheet.Using a nonlinear transformation of Eq. (B7) and applying series expansion of the integrand in Eq. (B3) about ξ → ∞ and integrating over dω ′ , we get where a modified-moment function G m (λ) is introduced as a linear combination of the raw moments ⟨ω i ⟩ up to order m (0 see Tab.III.From the inverse of Eq. (B7) we obtain (using "+" for ω ≥ 0 and "−" for ω < 0) Although the radius of convergence of the above series is increased compared to that of Eq. (B6), one can further improve the behavior of F (ω) for finite number (m ≤ K) of available moments ⟨ω m ⟩ in G m (λ) by employing the Padé approximants to the series in Eq. (B8), which do not have poles on the branch cut in the range −2λ < ω < 2λ.Using a diagonal Padé approximant P/P S(ξ;K,λ) to series S(ξ; K, λ) about the point ξ → ∞, we finally obtain  In practice, the lim ϵ→0 is replaced by a sufficiently small ϵ (e.g., ϵ/λ ∼ 0.001) to provide a smooth F (ω) in the range −2λ < ω < 2λ.

Appendix C: Finite element aggregate model
A polycrystalline aggregate model with 4000 randomly shaped and randomly oriented grains was generated upon Voronoi tessellation [18] with periodic microstructure in all three spatial directions [19].A FE mesh was generated with quadratic tetrahedral elements (∼5 million total elements of type C3D10) to preserve the geometry of the grains.The aggregate model used in this study is shown in Fig. 5.A general uniform loading Σ was applied to the aggregate such that Σ ij = ⟨σ ij ⟩ for the averages taken over the entire volume of the aggregate model.Since grains were assumed ideally elastic, a unit loading was applied using a small strain approximation.
The constitutive equations of the generalized Hooke's law were solved with FE solver Abaqus [20].Numerically calculated stress field σ ij was then used to obtain one σ nn value at each GB-element facet by first averaging the six Cauchy stresses σ ij at nearby Gauss points within the two elements sharing the facet and then projecting the averaged stress σij onto the facet normal n i , σ nn = ij σij n i n j .It was verified that the aggregate size and FE mesh density were sufficiently large to produce quasiisotropic responses and negligible finite size effects in the calculation of INS distributions 17 .
The grains were assumed to be made either of gamma iron (γ-Fe) or calcium sulfate (CaSO 4 ) with material properties listed in Tab.IV.This choice was made arbitrarily for the purpose of demonstrating the feasibility of the method in Sec.III.

FIG. 1 .
FIG.1.Feasibility of the method for the reconstruction of INS distribution PDF(σnn/Σmis) from available first K central moments.The original distribution (labeled FE) and its central moments µ m (m ≤ K) were calculated from the FE simulation of γ-Fe polycrystal under tensile loading.The µ m were then used to reconstruct the PDF using Eq.(18) with different parameter sets: (a) effect of K, (b) effect of P , (c) effect of λ (in units of Σmis).All the lines were calculated with ϵ/Σmis = 0.001 in the range of −2λ ≤ σnn − ⟨σnn⟩ ≤ 2λ.

F 2 .FIG. 5 .
FIG. 5. Periodic Voronoi aggregate model with 4000 grains and 31154 GBs.Different grains are denoted by different colors.FE mesh is shown for one selected grain.

TABLE I .
The first 20 central moments µ m dev of the general INS distribution for the assumed deviatoric stress loading Σ dev , expressed in terms of the two deviatoric stress invariants J2 ≡ tr(Σ 2 dev )/2 and J3 ≡ det(Σ dev ), and material invariants Mi,j with 2i + 3j = m.