# Green-FFT model for 3D contacts considering microstructure and defects

Lucas Fourel, Jean-Philippe Noyel, Xavier Kleber, Philippe Sainsot, Fabrice Ville

► To cite this version:

Lucas Fourel, Jean-Philippe Noyel, Xavier Kleber, Philippe Sainsot, Fabrice Ville. Green-FFT model for 3D contacts considering microstructure and defects. Mechanics & Industry, 2025, 26, pp.10.1051/meca/2025004. hal-05015275

HAL Id: hal-05015275

https://hal.science/hal-05015275v1

Submitted on 1 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Green-FFT model for 3D contacts considering microstructure and defects

Lucas Fourel<sup>1,2,3</sup>①, Jean-Philippe Noyel<sup>2,*</sup>①, Xavier Kleber<sup>3</sup>, Philippe Sainsot<sup>1</sup>, and Fabrice Ville<sup>1</sup>,

Received: 24 July 2024 / Accepted: 7 February 2025

**Abstract.** Physical phenomena at various scales are involved in energy dissipation and material degradation induced by mechanical contacts. Stress computation at the microstructural level is particularly valuable for the study of rolling contact fatigue or wear. For this purpose, a new efficient and versatile numerical model based on the Green-FFT method is proposed in this paper. Rough or dented surfaces can be modeled along with arbitrary material heterogeneities such as polycrystalline microstructure, inclusions or porosities. The novelty of this work lies in the efficiency of the Green-FFT method, which requires significantly lower computational cost than FEM for similar results. This enables detailed 3D analyses. The model entails three successive steps: contact pressure, stress distribution in a homogeneous material and stress distribution in a heterogeneous material computations. The periodization effect, which arises from the use of frequency domain, is mitigated by implementing zero-padding and buffer zones.

**Keywords:** FFT / 3D / contact / surface / material / microstructure

## 1 Introduction

The contacts between mechanical parts are essential to the operation of most vehicles and industrial equipment. They enable relative motion between different elements. Examples of contact components are rolling element bearings, gears, railroad tracks and wheels. Such components play a crucial role in the performance of mechanisms. Energy dissipation and material degradation are often concentrated at the contact level. Indeed, contact loading induces stresses that can cause various phenomena including plasticity, abrasive or adhesive wear and rolling contact fatigue.

Several models have been proposed to estimate the contact stress distribution. One of the first modeling approach was an analytical description of the contact using the theory of potential developed by Boussinesq [1] and Cerruti [2]. Hertz [3] established a seminal model that describes the contact between two smooth, homogeneous and elastic solids. Analytical expressions describing the stress distribution under these conditions were derived by Hills et al. [4].

However, surface defects such as roughness or dents, as well as subsurface defects such as inclusions or

porosities, can significantly affect the stress distribution. Additionally, various physical phenomena occur at the microstructural level. Indeed, contact components are typically made of metallic materials with a heterogeneous polycrystalline microstructure such as steel. Therefore, the previously mentioned homogeneous description can be unsatisfactory.

The versatility of Finite Element Method (FEM) allowed for the numerical computation of contact stress distributions in a broader variety of cases. The contact stress distributions in polycrystalline materials have mostly been computed in 2D due to the high computational cost of this method with fine meshes [5,6]. Other alternatives have been developed, such as the multigrid finite difference method [7,8].

Recently, a numerical method has gained attention due to its computational efficiency for determining stress distribution in heterogeneous solids [9]. This method is based on the concept of eigenstrain, which Eshelby [10,11] used in the case of elliptical inclusions and was later generalized by Mura [12]. Moulinec and Suquet [13,14] proposed a numerical method to compute stress fluctuations caused by heterogeneities in the frequency domain using the Fast Fourier Transform (FFT) algorithm. This method relies on superposing a fundamental solution called Green's function to compute the stress distribution.

<footnote n="1">¹ Univ Lyon, INSA Lyon, CNRS, LaMCoS, UMR5259, 69621 Villeurbanne, France</footnote>

<footnote n="2">² Univ Lyon, ECAM LaSalle, LabECAM, 69005 Lyon, France</footnote>

<footnote n="3">³ Univ Lyon, INSA Lyon, CNRS, MATEIS, UMR5510, 69621 Villeurbanne, France</footnote>

<footnote>* e-mail: jean-philippe.noyel@ecam.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Analysis domain.</figcaption></figure>

The method is often referred to as the “FFT method” and is closely related to certain implementations of the boundary element method. This paper uses the term “Green-FFT” to refer to the numerical method in order to differentiate it from the actual FFT algorithm.

In this paper, a numerical model based on the Green-FFT method is proposed to compute the stress distribution in a material with heterogeneities under contact conditions. An approach, combining multiple previous works, is used to compute the contact pressure, the stress distribution in a homogeneous material called homogeneous stress, and finally, the stress distribution in a heterogeneous material called heterogeneous stress. The originality of this work lies in the computational efficiency and versatility of the proposed model. The computational efficiency enables the study of physical phenomena at scales much finer than the contact scale by considering detailed microstructures and defects in 3D. The versatility of the approach allows the modeling of rough, dented or scratched surfaces and arbitrary shaped heterogeneities in the material.

First, the modeling strategy is introduced, followed by the presentation of the Green-FFT implementation. The accuracy and efficiency of the Green-FFT method are evaluated by comparison with FEM on equivalent models. Finally, various stress distributions with different surface and subsurface defects are presented to illustrate the capabilities of the proposed model.

## 2 Modeling strategy

The analysis of contact stress distributions is restricted to a small domain in one of the two contacting solids to avoid unnecessarily high computational cost. As illustrated in Figure 1, the second solid is implicitly modeled by the contact pressure distribution at the surface of the first solid. This approach is frequently used to calculate stress distribution at the contact scale [5,6,15].

The stress distribution in the analysis domain is obtained from three successive steps:

1. Computation of contact pressure $p$ from contact load and surface geometries.

2. Computation of homogeneous stress $\sigma^M$ from contact pressure.

3. Computation of heterogeneous stress $\sigma$ from homogeneous stress and microstructure.

Figure 2 illustrates the three computation steps and their main input parameters. The arrows indicate that the communication occurs exclusively in one direction. The addition of an iterative feedback loop, which could be a potential avenue for future development, is not considered in this study. The Green-FFT method can be employed at the different steps, leading to a significant reduction in computational cost compared to direct approaches. The numerical implementations of this method are introduced in the next section.

## 3 Numerical methods

The Green-FFT method is used to compute contact pressure, homogeneous stress and heterogeneous stress. First, the general principles of the method are presented. Then, the peculiarities of each step are detailed.

As previously stated, these methods are not novel, and the originality of this paper lies in their combination. This section provides readers with the principles for better understanding. For more details, refer to Wang et al.'s review on FFT-based contact solvers [16] and Lucarini et al.'s review on FFT-based heterogeneous solvers [17].

### 3.1 Green-FFT method

The Green-FFT method is a numerical method that solves differential equations. This method is based on the concepts of Green's function, convolution and frequency domain.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Computation steps.</figcaption></figure>

Green's functions are mathematical tools used to solve differential equations. A linear differential equation can be expressed using a linear differential operator L which gives the relationship between the response u and the source p:

$$Lu = p \quad (1)$$

George Green introduced the Green's functions in 1828 [18]. These functions are related to the concepts of influence coefficients, propagators, fundamental solutions and elementary solutions. Depending on the scientific field and the authors, these terms may refer to the functions themselves or to derived operators. A Green's function G, is the solution of a given differential equation under a unit impulse δ, also referred as Dirac Delta distribution, located at a specific point in the domain with given boundary conditions, satisfying the following relationship:

$$LG = \delta \quad (2)$$

The linearity of the differential operator L allows the superposition of the solutions of equation (1). Since any source distribution p is a superposition of concentrated sources at different points, the response u can be obtained by the same superposition of the Green's function G. This superposition is achieved using the convolution operation:

$$u(x) = (p * G)(x) = \int_{-\infty}^{+\infty} p(x') G(x - x') dx' \quad (3)$$

Computing the discrete convolution directly requires n² operations, where n represents the number of elements in the domain. The number of operations can be reduced by

performing computation in the frequency domain instead of the space domain.

A function u defined in the space domain with coordinates x can be expressed in the frequency domain with coordinates ξ using the Fourier transform F [19]:

$$\mathcal{F}(u) : \xi \mapsto \hat{u}(\xi) = \int_{-\infty}^{+\infty} u(x) e^{-i\xi x} dx \quad (4)$$

The convolution of functions expressed in the space domain becomes an algebraic product in the frequency domain:

$$\begin{gathered} \mathcal{F}(p * G) = \mathcal{F}(p) \mathcal{F}(G) \\ \hat{u}(\xi) = \hat{p}(\xi) \hat{G}(\xi) \end{gathered} \quad (5)$$

The computation of u from p can be performed using the frequency domain, with a complexity of O(n log n). The computational efficiency of the Green-FFT method arises from the lower complexity of the frequency method for the computation of discrete convolutions. As a result, this method is highly scalable as model size increases.

However, using the frequency domain induces a periodization of the space domain that is undesirable in some cases. This periodization error can be minimized using zero-padding, which consists of extending the space domain and imposing a null source on the extension.

The Green-FFT method is limited to structured meshes with constant element size. It is therefore unsuitable for problems with complex geometries or those requiring adaptive discretization. On the other hand, the numerical efficiency of the Green-FFT method makes it particularly relevant for problems where a structured mesh is
<figure><img src="image_4.png" /><figcaption>Fig. 3. Initial gap $h_0$, gap $h$, displacement $u$, approach $\delta$.</figcaption></figure>

an option, such as contacts and heterogeneous materials. The grid mesh may not conform to the geometry of the heterogeneities, however, this limitation can be mitigated by increasing the mesh resolution.

## 3.2 Contact pressure

The method for computing the contact pressure used in the numerical model is presented in this section. This method has been proposed by Polonsky and Keer [20], Sainsot and Lubrecht [21], Liu et al. [22] and Frérot et al. [23], among others. The method and its variations are described in detail in the papers cited above. The paper of the contact mechanics challenge compares several alternative methods [24].

The solver uses the Green-FFT method, combined with the conjugate gradient method [25].

The contact occurs between two deformable bodies. The first body has a Young's modulus $E_1$ and a Poisson's ratio $\nu_1$ and the second body has a Young's modulus $E_2$ and a Poisson's ratio $\nu_2$. The equivalent contact representation is used in this analysis to reduce the number of parameters [26]. In the equivalent representation, the first body is semi-infinite with equivalent elastic modulus $E^*$:

$$ \frac{1}{E^*} = \frac{1 - \nu_1^2}{E_1} + \frac{1 - \nu_2^2}{E_2} \qquad (6) $$

The second equivalent body is non-deformable (disk in Fig. 3).

The main modeling assumptions include the homogeneity and isotropy of the linear elastic semi-infinite solid, alongside frictionless and dry contact.

When a load $W$ is transmitted through the contact, an approach $\delta$ occurs. The surfaces deform to prevent inter-penetration. Figure 3 illustrates the relationship between the gap $h$, the initial gap $h_0$, the normal displacement $u$ and the approach $\delta$:

$$ h(x, y) = h_0(x, y) + u(x, y) - \delta \qquad (7) $$

The contact pressure $p$ is positive where the gap $h$ is null and the contact pressure is null where the gap is positive:

$$ \begin{cases} h(x, y) > 0 \rightarrow p(x, y) = 0 \\ h(x, y) = 0 \rightarrow p(x, y) > 0 \end{cases} \qquad (8) $$

According to the force equilibrium, the integral of the contact pressure on the surface equals the load $W$:

$$ W = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} p(x, y) \, dx \, dy \qquad (9) $$

The complementary elastic energy is derived from the geometric compatibility (Eq. (7)), the contact conditions (Eq. (8)) and the force equilibrium (Eq. (9)) [27]:

$$ U_c = \int_{\Omega_c} \left( \frac{1}{2} p u + p h_0 \right) dS \qquad (10) $$

Determining the contact pressure involves finding the distribution $p$ that minimizes the complementary elastic energy which has a unique solution $U_c$ [27,28].

The Green-FFT method described in the previous section is used to calculate the displacement $u(x, y)$ resulting from a pressure distribution $p(x, y)$. The expression of the general equation (1) applied in this step is:

$$ Lu(x, y) = p(x, y) \qquad (11) $$

The relationship between the concentrated contact pressure $p$ located in $(x = 0, y = 0)$ and the normal displacement $u$ is given by the Green's function $G^{p \to u}$. An analytical expression has been derived by Boussinesq [1] and Cerruti [2]:

$$ u(x, y) = p(x, y) * G^{p \to u}(x, y) \qquad (12) $$

$$ G^{p \to u}(x, y) = \frac{1}{\pi E^*} \frac{1}{\sqrt{x^2 + y^2}} \qquad (13) $$
With the equivalent elastic modulus $E^*$ defined in equation (6).

This Green's function is then integrated to obtain the discrete convolution kernel $K^{p→u}$ defined in Appendix A.

The solution is computed iteratively using a conjugate gradient method proposed by Polonsky and Keer [20]. The Green-FFT method is used at each iteration step to obtain the displacement $u$ induced by the contact pressure $p$ using the discrete convolution kernel $K^{p→u}$. In case of a non-periodic contact, zero-padding is performed on an extended domain. The extended domain is twice the size of the analysis domain in each non-periodic dimension.

### 3.3 Homogeneous stress

The homogeneous stress $\sigma^M$ is computed with the Green-FFT method using the contact pressure $p$ as the source. The expression of the general equation (1) applied in this step is:

$$L\sigma^M(x, y, z) = p(x, y) \quad (14)$$

The relationship between the concentrated contact pressure $p$ located in $(x = 0, y = 0)$ and the homogeneous stress $\sigma^M$ is given by the Green's function $G^{p→\sigma}$:

$$\sigma_{ij}^M(x, y, z) = p(x, y) * G_{ij}^{p→\sigma}(x, y, z) \quad (15)$$

Love [29] and Kalker [27] derived the Green's function $G_{ij}^{p→\sigma}$ and the discrete convolution kernel $K_{ij}^{p→\sigma}$. Their analytical expressions are defined in Appendix B. Other authors also used this method to compute homogeneous stress distributions under contact conditions [30–35].

The computation scheme is:

$$\left\{
\begin{aligned}
\hat{K}_{ij}^{p→\sigma} &= \text{FFT}(K_{ij}^{p→\sigma}) \\
\hat{p} &= \text{FFT}(p) \\
\hat{\sigma}_{ij}^{M} &= \hat{p} \hat{K}_{ij}^{p→\sigma} \\
\sigma_{ij}^{M} &= \text{iFFT}(\hat{\sigma}_{ij})
\end{aligned}
\right.
\quad (16)$$

### 3.4 Heterogeneous stress

Local stress variations are caused by heterogeneities such as microstructure, phases, inclusions and fibers. The Green-FFT method is used in this model to calculate heterogeneous stress based on homogeneous stress. This method, applied to heterogeneous materials, was introduced by Moulinec and Suquet [13,14] and is based on the concepts of eigenstrain and polarization stress.

Eigenstrain refers to nonelastic strains such as thermal expansion strain, phase transformation strains, plastic strains, initial strains or strains due to heterogeneities. The total heterogeneous strain $\varepsilon$ is considered to be the sum of the homogeneous strain $\varepsilon^M$ and the eigenstrain also called fluctuation strain $\tilde{\varepsilon}$:

$$\varepsilon_{ij} = \varepsilon_{ij}^{M} + \tilde{\varepsilon}_{ij} \quad (17)$$

$\varepsilon_{ij}^M$ represents the strain computed from the homogeneous stress $\sigma_{ij}^M$.

Polarization stress in the context of a heterogeneous material was introduced by Kroner [36] to refer to stresses induced by heterogeneities. In this model the polarization stress $\zeta$ is defined by:

$$\zeta_{ij} = \Delta C_{ijkl} \varepsilon_{kl} \quad (18)$$

With:

$$\Delta C_{ijkl} = C_{ijkl} - C_{ijkl}^0 \quad (19)$$

With the heterogeneous material elasticity tensor $C$ varying at each point in the domain depending on the local material properties and a homogeneous elasticity tensor $C^0$ that is constant for all points in the domain.

The fluctuation strain $\hat{\varepsilon}$ is computed in the frequency domain with coordinates $\xi$ from the polarization stress $\hat{\zeta}$ using:

$$\hat{\varepsilon}_{ij} = \frac{1}{2} (\hat{\Gamma}_{ijkl}^{\zeta \to \tilde{\varepsilon}} + \hat{\Gamma}_{jikl}^{\zeta \to \tilde{\varepsilon}}) \hat{\zeta}_{kl} \quad (20)$$

With the Green's operator $\hat{\Gamma}_{ijkl}$:

$$\hat{\Gamma}_{ijkl} = \hat{G}_{ik,lj} = -\hat{G}_{ik} \xi_j \xi_l \quad (21)$$

With the Green's function $\hat{G}_{ik}$:

$$\hat{G}_{ik} = (C_{ijkl}^0 \xi_j \xi_l)^{-1} \quad (22)$$

The fixed-point computation scheme proposed by Moulinec and Suquet [13,14] is used:

$$\text{Initialization: } \varepsilon_{ij}^{M} = (C_{ijkl}^{0})^{-1} \sigma_{kl}^{M} \\
\varepsilon_{ij}^{(0)} = \varepsilon_{ij}^{M}$$

$$\text{Iteration: } \zeta_{ij}^{(n)} = \Delta C_{ijkl} \varepsilon_{kl}^{(n-1)} \\
\hat{\zeta}_{ij}^{(n)} = \text{FFT}(\zeta_{ij}^{(n)}) \quad (23) \\
\hat{\varepsilon}_{ij}^{(n)} = \frac{1}{2} (\hat{\Gamma}_{ijkl}^{\zeta \to \tilde{\varepsilon}} + \hat{\Gamma}_{jikl}^{\zeta \to \tilde{\varepsilon}}) \hat{\zeta}_{kl}^{(n)} \\
\varepsilon_{ij}^{(n)} = \varepsilon_{ij}^{M} + \text{iFFT}(\hat{\varepsilon}_{ij}^{(n)})$$

Convergence test

$$\text{Result: } \sigma_{ij} = C_{ijkl} \varepsilon_{kl}^{(n)}$$

With the following convergence test:

$$\frac{\langle \|\varepsilon_{ij}^{(n)} - \varepsilon_{ij}^{(n-1)}\| \rangle}{\langle \|\varepsilon_{ij}^{M}\| \rangle} < e \quad (24)$$
<table><caption>**Table 1.** Homogeneous properties of the solid.</caption><tr><td>Young's modulus <i>E</i></td><td>212 GPa</td></tr><tr><td>Poisson's ratio <i>ν</i></td><td>0.296</td></tr></table>

With the L2 norm $||·||$, the spatial arithmetic mean $⟨·⟩$ and the tolerance threshold $e$.

The discretization scheme used in this study does not strictly converge to numerical precision for infinite stiffness contrast. However, sufficiently accurate solutions can be obtained using a tolerance threshold of $10^{-3}$. This value is used in the following computations. More advanced discretization schemes such as the one proposed by Willot [37] can be used to circumvent this potential limitation.

This method is applied to the numerical simulation of stress in polycrystalline materials by several authors [38–41].

For non-periodic problems, the analysis domain is extended by adding a buffer zone. This solution is employed by several researchers to address non-periodic problems [9,38,42–44]. The effect of the buffer properties and the influence of the buffer size on the periodicity error are analyzed in Section 4.

## 4 Buffer zone

As previously stated, the use of the frequency domain in the Green-FFT method induces a periodization of the analysis domain. Although this can be convenient for periodic problems, it introduces a periodization error for non-periodic problems.

The Green's functions are expressed in the space domain for the first two steps, which consist of the computation of contact pressure and homogeneous stress. This allows the use of zero-padding over a doubled analysis domain to eliminate the periodization error. However, for the final step consisting of the computation of the heterogeneous stress, the Green's function is expressed in the frequency domain which implies a non-zero periodization error. Nevertheless, by extending the analysis domain with a buffer zone, the periodization error can be negligible [38].

In order to determine the size and properties of the buffer zone, a point contact loading under Hertz conditions is modeled. Therefore, the contact surfaces are perfectly smooth and have constant curvatures, resulting in a circular contact area of radius $a$ and a maximum contact pressure $p_0$ at the center. The lengths are made dimensionless using the contact half-width $a$ and the stress distributions using the Hertz pressure $p_0$. The modeling strategy described in Section 2 is applied. First, the contact pressure is calculated, then homogeneous stress is calculated and finally, heterogeneous stress is computed. The analysis domain is a cube of edge length $L = 3a$ discretized into $256 \times 256 \times 256$ elements.

The isotropic elastic properties used to calculate homogeneous stress are presented in Table 1.

The solid used to compute heterogeneous stress is a polycrystal with a random microstructure generated by

<table><caption>**Table 2.** Cubic elastic properties of α-iron [47].</caption><tr><td><i>C</i><sub>11</sub></td><td>237 GPa</td></tr><tr><td><i>C</i><sub>12</sub></td><td>141 GPa</td></tr><tr><td><i>C</i><sub>44</sub></td><td>116 GPa</td></tr></table>

Voronoi tessellation using Neper [45]. The grains have a cubic elastic behavior described by the elasticity tensor $C^c$ which can be defined by a matrix using the Voigt notation [46]:

$$C^c = \begin{bmatrix} C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{12} & C_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{44} \end{bmatrix} \quad (25)$$

With the cubic elastic properties $C_{11}$, $C_{12}$ and $C_{44}$. Table 2 presents the values provided by Courtney [47] for α-iron which are used in this study.

An individual crystal orientation is assigned to each grain, allowing the elasticity tensor $C$ to be oriented as a function of position in the solid:

$$C_{i_1 j_1 k_1 l_1}(x, y, z) = R_{i_1 i_2}(x, y, z) R_{j_1 j_2}(x, y, z) \times R_{k_1 k_2}(x, y, z) R_{l_1 l_2}(x, y, z) C^c_{i_2 j_2 k_2 l_2} \quad (26)$$

With the rotation matrix $R$, which depends on the grains and therefore on the location in the solid.

Below the surface, the same periodic microstructure is extended to ensure a continuity of the grains between the analysis domain and the buffer zone. The same behavior and properties as the solid are imposed in the buffer zone. Regarding the properties of the buffer zone above the surface, two options are available:

* A perfectly compliant behavior is prescribed above the surface as illustrated in Figure 4b.

* The solid behavior and periodic microstructure are identical above the contact area, while above the free surface the behavior is perfectly compliant as illustrated in Figure 4c.

As shown in Figures 4b and 5, the fully compliant buffer zone of a size $L_b$ above the surface causes the normal surface stress $\sigma_{zz}$ to converge to the contact pressure $p$. Figures 4c and 5 indicate that in contrast, if the solid behavior is applied above the contact area, normal stress fluctuations can build up at the surface and $\sigma_{zz}$ converges to a contact pressure, modified by the properties of the microstructure. Figures 4a and 5 show that without a buffer zone, the periodicity induced by the frequency domain causes the bottom of the domain to interact directly with the top, allowing normal stress fluctuations in the contact area, but also at the free surface. The choice among the two buffer types depends on the modeling approach employed: a fully compliant buffer zone above the surface to enforce the contact pressure or a solid buffer
<figure><img src="image_5.png" /><figcaption>Fig. 4. Heterogeneous stress computation using the Green-FFT method (a) without buffer zone, (b) with a compliant buffer zone above the surface, (c) with a solid buffer zone above the contact area.</figcaption></figure>

zone above the contact area to account for the fluctuations in contact pressure due to material heterogeneities.

To estimate the periodization error $e$, the stress distributions $\sigma|_{L_b}$ obtained with various compliant buffer sizes $L_b$ are compared to the reference stress distribution $\sigma|_{ref}$ obtained with a large buffer size of $L_b = 5 D_G$, where

$D_G$ represents the grain size. The periodization error $e$ is defined by the following relationship:

$$e_{ij} = \left| \frac{\sigma_{ij}|_{L_b} - \sigma_{ij}|_{ref}}{p_0} \right| \quad (27)$$
<figure><img src="image_6.png" /><figcaption>Fig. 5. Heterogeneous normal stress $\sigma_{zz}$ in surface elements at the center of the point contact.</figcaption></figure>

Figure 6a illustrates the periodization error without buffer zone ($L_b = 0$) for a grain size $D_G = 0.08$ a. The periodization error is concentrated at the domain boundaries due to the field incompatibilities. Figure 6b shows a decrease of the maximum periodization error as the buffer size $L_b$ increases compared to the grain size $D_G$. The analysis of three different grain sizes $D_G$ indicates that the maximum periodization error does not directly depend on the size of the analysis domain, but on the characteristic size of the heterogeneities. The stress distributions converge as the buffer zone reaches the size of several grains. In this case, the maximum periodization error is lower than 1% when $L_b = 2 D_G$. To minimize the effect of periodization when using the Green-FFT method, this value is used throughout the study.

## 5 Method evaluation

The proposed Green-FFT model computes the heterogeneous stress distribution under contact conditions. This information is difficult to acquire experimentally, therefore the accuracy of the Green-FFT method is evaluated by comparison with the Finite Element Method (FEM), which is extensively used in solid mechanics computations. An equivalent FEM model has been developed for this purpose, using the same microstructure, elastic properties, mesh and contact conditions as the Green-FFT model to ease comparisons. The efficiency of the Green-FFT method is then evaluated for different model size.

Conducting a comparison between Green-FFT and FEM in 3D while preserving the geometry of the

microstructure with a structured mesh is too computationally expensive with FEM. Therefore, the comparison is conducted on equivalent Green-FFT and FEM models in 2D with plane strain hypothesis. The 2D FEM model is similar to other existing models [6,48–51]. As shown in Figure 7, the domain of the FEM model is extended by height adjacent domains with the same periodic microstructure. This allows to minimize the influence of boundary conditions by moving the boundaries of the extended domain away from the analysis domain. Zero normal displacements are imposed at the boundaries of the extended domain, except at the surface ($z = 0$) where the contact stress $\sigma_{zz} = -p$ is imposed and the other components are set to zero. First-order polynomial elements are used.

The difference $\Delta\sigma$ is defined as the difference between the stress distribution $\sigma|_{\text{Green-FFT}}$ obtained with Green-FFT and the stress distribution $\sigma|_{\text{FEM}}$ obtained with FEM, normalized by the maximum contact pressure $p_0$:

$$ \Delta\sigma_{ij} = \left| \frac{\sigma_{ij}|_{\text{Green-FFT}} - \sigma_{ij}|_{\text{FEM}}}{p_0} \right| \quad (28) $$

Since the contact pressure is imposed in the FEM model, a fully compliant buffer zone above the surface is used in the Green-FFT model. The contact pressure is determined with Hertz analytical model in the FEM model while it is computed numerically in the Green-FFT model. This allows to reduce the uncertainty in the FEM model and to evaluate all three steps of the Green-FFT model. The analysis domain is a square of size 5 a. The element size is set to 0.01 a.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** (a) Periodization error $e_{zz}$ without buffer zone $L_b = 0$ for a grain size $D_G = 0.080\ a$, (b) maximum periodization error depending on the buffer size $L_b$ relative to the grain size $D_G$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Equivalent FEM model.</figcaption></figure>

Both methods result in similar heterogeneous stress distributions as shown in Figure 8. The maximum difference is lower than 5%. The difference is highest at the grain boundaries, where the interfaces between the solid properties cause slight oscillations. Although solutions to reduce

this effect have been proposed in the literature [37], they have not been implemented in this model, as the results are considered satisfactory.

Different mesh resolutions are used to compare the execution time of Green-FFT and FEM as a function of
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison between the Green-FFT model and the equivalent FEM model.</figcaption></figure>

degrees of freedom (DOF). Although the execution time may vary depending on the hardware and software implementation, the simulations were performed on the same hardware and using optimized libraries for core functions such as the tensor product, solving the linear system of equations, and FFT.

Figure 9 shows the efficiency of the Green-FFT method compared to the finite element method. Within this model size range, the Green-FFT method provides roughly 30 times higher efficiency than the FEM with the same computational resources. The Green-FFT method successfully solves 3 million DOF in roughly 3 s on a laptop. This efficiency gain allows consideration of much larger 3D models, which would have been prohibitively expensive using FEM. The following section presents multiple 3D case studies with several hundred millions DOF.

## 6 Surface and subsurface defect simulations

In this section, multiple stress distributions are computed in order to illustrate the capabilities of the proposed Green-FFT model. The point contact described in Section 4 is used as the contact loading. Therefore, the lengths and stress distributions are made dimensionless using the contact half-width $a$ and the Hertz pressure $p_0$, respectively.

The first example incorporates a spherical inclusion of diameter $D_i = 0.2a$ into the polycrystalline solid at the point $(x = 0, y = 0, z = -0.48a)$ which corresponds to the location of maximum shear stress for Hertz point contacts. An isotropic elastic behavior is assigned to the inclusion,

<table><caption>Table 3. Properties of the inclusion.</caption><tr><td>Young's modulus <i>E</i></td><td>424 GPa</td></tr><tr><td>Poisson's ratio <i>ν</i></td><td>0.3</td></tr></table>

<table><caption>Table 4. Simulation parameters.</caption><tr><td>Size of analysis domain</td><td>3a × 3a × 2a</td></tr><tr><td>Number of elements</td><td>480 × 480 × 320</td></tr><tr><td>Element size</td><td>6.25 × 10<sup>-3</sup> a</td></tr><tr><td>Number of grains</td><td>100 000</td></tr><tr><td>Grain size</td><td>7.01 × 10<sup>-2</sup> a</td></tr></table>

which is two times stiffer than that of the solid. The elastic properties of the inclusion are presented in Table 3.

A variation of this case is also simulated. This second example is identical to the first, except that the Young's modulus of the inclusion is set to zero to model a spherical porosity.

The simulation parameters common to all cases are presented in Table 4.

The resulting maximum shear stress distributions are shown in Figure 10a for the inclusion and Figure 10b for the porosity.

The third example involves contact between a smooth and a dented surface. The same load and equivalent curvature radius of the surfaces as in the previous examples are used, however, the contacting surface is modified with a circular dent profile. The dent profile is modeled by the analytical expression proposed by Coulon et al. [52]:

$$z(r) = Z_d \exp\left(-K_d \frac{r^2}{4D_d^2}\right) \cos\left(\frac{\pi r}{D_d}\right) \quad (29)$$
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Execution time of the Green-FFT model and the equivalent FEM model depending on the degrees of freedom (DOF) of the models.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Maximum shear stress with a subsurface defect, (a) inclusion, (b) porosity.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Dent profile.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Maximum shear stress with a dented surface.</figcaption></figure>

<table><caption>**Table 5.** Dent parameters.</caption><tr><td>Dent maximum depth Z<sub>d</sub></td><td>15 × 10<sup>-6</sup> a</td></tr><tr><td>Dent diameter D<sub>d</sub></td><td>0.5 a</td></tr><tr><td>Dent shoulder parameter K<sub>d</sub></td><td>12</td></tr></table>

The resulting maximum shear stress distribution is presented in Figure 12.

## 7 Conclusion

A new numerical model has been proposed to compute the stress distribution under contact conditions, taking into account surface defects as well as material heterogeneities such as microstructure and subsurface defects. The model combines multiple existing methods using various Green's functions and the FFT algorithm. This enables the resulting Green-FFT model to handle detailed 3D microstructures that were not possible with conventional methods such as FEM.

$$r(x, y) = \sqrt{x^2 + y^2} \qquad (30)$$

With dent maximum depth $Z_d$, dent diameter $D_d$ and dent shoulder parameter $K_d$. The dent profile and parameters are illustrated in Figure 11.

The values of the dent parameters that are used in the example are presented in Table 5.
Zero-padding and buffer zones are employed to minimize the periodization effect caused by the use of the frequency domain.

A comparison was made between equivalent Green-FFT and FEM models to assess the accuracy and efficiency of the proposed model. Both models produced comparable stress distributions. Nevertheless, the Green-FFT method demonstrated a significantly better efficiency than FEM.

Examples were presented to demonstrate how the Green-FFT model effectively handles a wide variety of cases. Large models consisting of several hundred million DOF can be computed on a personal computer in less than an hour. However, these cases require significant memory resources.

The proposed model is applicable to the study of microscopic physical phenomena such as rolling contact fatigue [53]. The finer resolution made possible by the efficiency of the method allows the development of new microscopic fatigue criteria. In future works, implementation of other material behaviors, such as plasticity and damage, could extend the applicability of this Green-FFT model to the study of a wider range of phenomena related to contact mechanics.

## Funding

This research received no external funding.

## Conflicts of interest

The authors declare that they have no conflict of interest.

## Data availability statement

This article has no associated data generated.

## Author contribution statement

All authors participated in the study and agreed with this version of the manuscript. Conceptualization, all authors; Methodology, all authors; Software, Lucas Fourel and Philippe Sainsot; Writing – Original Draft Preparation, Lucas Fourel, Jean-Philippe Noyel and Fabrice Ville.

## References

[1] J. Boussinesq, Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques, Gauthier-Villars, Paris, 1885

[2] V. Cerruti, Ricerche Intorno All'equilibrio de'corpi Elastici Isotropi. Reale Accaemia dei Lincei, Roma, 1882

[3] H. Hertz, über die Berührung fester elastischer Körper, Zeitung Reine Angew. Math. **92**, 156–171 (1882)

[4] D.A. Hills, D. Nowell, A. Sackfield, Mechanics of Elastic Contacts. Butterworth-Heinemann, Oxford [England]; Boston, 1993

[5] N.R. Paulson, J.A.R. Bomidi, F. Sadeghi, R.D. Evans, Effects of crystal elasticity on rolling contact fatigue, Int. J. Fatigue **61**, 67–75 (2014)

[6] J.-P. Noyel, F. Ville, P. Jacquet, A. Gravouil, C. Changenet, Development of a granular cohesive model for rolling contact fatigue analysis: crystal anisotropy modeling, Tribol. Trans. **59**, 469–479 (2016)

[7] H. Boffy, C.H. Venner, Multigrid solution of the 3D stress field in strongly heterogeneous materials, Tribol. Int. **74**, 121–129 (2014)

[8] B. Zhang, H. Boffy, C.H. Venner, Multigrid solution of 2D and 3D stress fields in contact mechanics of anisotropic inhomogeneous materials, Tribol. Int. **149**, 105636 (2020)

[9] R.A. Lebensohn, A.D. Rollett, Spectral methods for full-field micromechanical modelling of polycrystalline materials, Computat. Mater. Sci. **173**, 109336 (2020)

[10] J.D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. Roy. Soc. Lond. Ser. A Math. Phys. Sci. **241**, 376–396 (1957)

[11] J.D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proc. Roy. Soc. Lond. Ser. A Math. Phys. Sci. **252**, 561–569 (1959)

[12] T. Mura, Micromechanics of Defects in Solids. No. 3 in Mechanics of Elastic and Inelastic Solids, 2 edn., M. Nijhoff; Distributors for the U.S. and Canada, Kluwer Academic Publishers, Dordrecht, Netherlands ; Boston : Hingham, MA, USA, 1987

[13] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, Comptes Rendus Acad. Sci. Sér II Méc. Phys. Chim. Astron. **318**, 1417–1423 (1994).

[14] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Eng. **157**, 69–94 (1998)

[15] F. Sadeghi, B. Jalalahmadi, T.S. Slack, N. Raje, N.K. Arakere, A review of rolling contact fatigue, J. Tribol. **131**, 041403 (2009)

[16] Q.J. Wang, L. Sun, X. Zhang, S. Liu, D. Zhu, FFT-based methods for computational contact mechanics, Front. Mech. Eng. **6**, 61 (2020).

[17] S. Lucarini, M.V. Upadhyay, J. Segurado, FFT based approaches in micromechanics: fundamentals, methods and applications, Model. Simul. Mater. Sci. Eng. **30**, 023002 (2022).

[18] G. Green, An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism, T. Wheelhouse, Nottingham, 1828

[19] J. Fourier, Théorie Analytique de la Chaleur, Firmin Didot, Paris, 1822

[20] I. Polonsky, L. Keer, A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques, Wear **231**, 206–219 (1999)

[21] P. Sainsot, A.A. Lubrecht, Efficient solution of the dry contact of rough surfaces: a comparison of fast Fourier transform and multigrid methods, Proc. Inst. Mech. Eng. J **225**, 441–448 (2011).

[22] S. Liu, Q. Wang, G. Liu, A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses, Wear **243**, 101–111 (2000)

[23] L. Frérot, G. Anciaux, V. Rey, S. Pham-Ba, J.-F. Molinari, Tamaas: a library for elastic-plastic contact of periodic rough surfaces, J. Open Source Softw. **5**, 2121 (2020)

[24] M.H. Müser, W.B. Dapp, R. Bugnicourt, P. Sainsot, N. Lesaffre, T.A. Lubrecht, B.N.J. Persson, K. Harris, A. Bennett, K. Schulze, S. Rohde, P. Ifju, W.G. Sawyer, T. Angelini, H. Ashtari Esfahani, M. Kadkhodaei, S. Akbarzadeh, J.-J. Wu, G. Vorlaufer, A. Vernes,
S. Solhjoo, A.I. Vakis, R.L. Jackson, Y. Xu, J. Streator, A. Rostami, D. Dini, S. Medina, G. Carbone, F. Bottiglione, L. Afferrante, J. Monti, L. Pastewka, M.O. Robbins, J.A. Greenwood, Meeting the contact-mechanics challenge, Tribol. Lett. **65**, 118 (2017)

[25] M. Hestenes, E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bureau Standards **49**, 409 (1952)

[26] K.L. Johnson, Contact Mechanics. Cambridge University Press, Cambridge, 1987

[27] J.J. Kalker, Three-Dimensional Elastic Bodies in Rolling Contact, vol. 2 of Solid Mechanics and Its Applications. Springer Netherlands, Dordrecht, 1990

[28] G. Duvaut, JL. Lions, Les inéquations en mécanique et en Physique, Dunod, Paris, 1972

[29] A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, 4th edn., Cambridge University Press, 1927

[30] C. Mayeur, P. Sainsot, L. Flamand, A numerical elasto-plastic model for rough contact, J. Tribol. **117**, 422–429 (1995)

[31] P. Sainsot, C. Jacq, D. Nelias, A numerical model for elasto-plastic rough contact, Comput. Model. Eng. Sci. **3**, 497–506 (2002)

[32] F. Wang, L.M. Keer, Numerical simulation for three dimensional elastic-plastic contact with hardening behavior, J. Tribol. **127**, 494–502 (2005)

[33] K. Zhou, W.W. Chen, L.M. Keer, Q.J. Wang, A fast method for solving three-dimensional arbitrarily shaped inclusions in a half space, Comput. Methods Appl. Mech. Eng. **198**, 885–892 (2009)

[34] J. Leroux, B. Fulleringer, D. Nelias, Contact analysis in presence of spherical inhomogeneities within a half-space, Int. J. Solids Struct. **47**, 3034–3049 (2010)

[35] K. Zhou, W. Wayne Chen, L.M. Keer, X. Ai, K. Sawamiphakdi, P. Glaws, Q. Jane Wang, Multiple 3D inhomogeneous inclusions in a half space under contact loading, Mech. Mater. **43**, 444–457 (2011)

[36] E. Kröner, Berechnung der elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls, Z. Phys. **151**, 504–518 (1958)

[37] F. Willot, Fourier-based schemes for computing the mechanical response of composites with accurate local fields, Comptes Rendus Méc. **343**, 232–245 (2015)

[38] R.A. Lebensohn, R. Brenner, O. Castelnau, A.D. Rollett, Orientation image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper, Acta Materialia **56**, 3914–3926 (2008)

[39] P. Suquet, H. Moulinec, O. Castelnau, M. Montagnat, N. Lahellec, F. Grennerat, P. Duval, R. Brenner, Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep, Procedia IUTAM **3**, 76–90 (2012)

[40] J.-B. Gasnier, F. Willot, H. Trumel, B. Figliuzzi, D. Jeulin, M. Biessy, A Fourier-based numerical homogenization

tool for an explosive material, Matér. Tech. **103**, 308 (2015)

[41] A. Schmidt, C. Gierden, J. Waimann, B. Svendsen, S. Reese, Two-scale FE-FFT-based thermo-mechanically coupled modeling of elasto-viscoplastic polycrystalline materials at finite strains, PAMM **22**, e202200172 (2023)

[42] N.B. Nkoumbou Kaptchouang, L. Gelebart, Extension des solveurs FFT aux conditions aux limites non périodiques et aux méthodes multi-grilles locales, in: 15ème Colloque National En Calcul Des Structures, (83400 Hyères-Palmiers, France), Université Polytechnique Hauts-de-France [UPHF], May 2022.

[43] C. Cocke, H. Mirmohammad, M. Zecevic, B. Phung, R. Lebensohn, O. Kingstedt, A. Spear, Implementation and experimental validation of nonlocal damage in a large-strain elasto-viscoplastic FFT-based framework for predicting ductile fracture in 3D polycrystalline materials, Int. J. Plast. **162**, 103508 (2023)

[44] M. Zecevic, R.A. Lebensohn, L. Capolungo, Non-local large-strain FFT-based formulation and its application to interface-dominated plasticity of nano-metallic laminates, J. Mech. Phys. Solids **173**, 105187 (2023)

[45] R. Quey, P. Dawson, F. Barbe, Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing, Comput. Methods Appl. Mech. Eng. **200**, 1729–1745 (2011)

[46] W. Voigt, Lehrbuch Der Kristallphysik (Mit Ausschluss Der Kristalloptik), Teubner, Leipzig, Berlin, 1910

[47] T.H. Courtney, Mechanical Behavior of Materials, 2nd edn. Waveland Press, Long Grove, Illinois, 2005

[48] E. Bossy, J. Noyel, X. Kleber, F. Ville, C. Sidoroff, S. Thibault, Competition between surface and subsurface rolling contact fatigue failures of nitrided parts: a Dang Van approach, Tribol. Int. **140**, 105888 (2019)

[49] G. Vouaillat, J.-P. Noyel, F. Ville, X. Kleber, S. Rathery, From Hertzian contact to spur gears: Analyses of stresses and rolling contact fatigue, Mech. Ind. **20**, 626 (2019)

[50] L. Fourel, J.-P. Noyel, E. Bossy, X. Kleber, P. Sainsot, F. Ville, Towards a grain-scale modeling of crack initiation in rolling contact fatigue – Part 1: shear stress considerations, Tribol. Int. **164**, 107224 (2021)

[51] L. Fourel, J.-P. Noyel, E. Bossy, X. Kleber, P. Sainsot, F. Ville, Towards a grain-scale modeling of crack initiation in rolling contact fatigue – Part 2: persistent slip band modeling, Tribol. Int. **163**, 107173 (2021)

[52] S. Coulon, F. Ville, and, A.A. Lubrecht, Effect of a dent on the pressure distribution in dry point contacts, J. Tribol. **124**, 220–223 (2002)

[53] L. Fourel, J.-P. Noyel, X. Kleber, P. Sainsot, F. Ville, Numerical analysis of 3D crack initiation under rolling contact fatigue in polycrystals with surface dents and subsurface inclusions, in preparation, 2023
## Appendix A: Normal displacement in an elastic half-space subjected to a surface pressure

To compute the response $u$ numerically, the distribution $p$ and the Green's function $G$ must be defined in a discrete finite domain. The distributions are assumed constant in each element and computed at their centers $x_m$:

$$u(x) = u(x_m), \forall x \in \left[x_m - \frac{\Delta x}{2}, x_m + \frac{\Delta x}{2}\right] \qquad (A.1)$$

With the element size $\Delta x$. The discrete convolution is then defined as a sum:

$$u(x_m) = \sum_n p(x_n) K(x_m - x_n) \qquad (A.2)$$

With the discrete convolution kernel $K$:

$$K(x_m - x_n) = \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} G(x_m - x_n - x') dx' \qquad (A.3)$$

The components of the discrete convolution kernel $K$ correspond to the influence of a uniform source $p$ on an element. They are obtained analytically or numerically by integrating Green's function $G$.

In a discrete domain, the Green's function $G^{p \to u}$ defined in equation (13) is integrated over each rectangular element of size $\Delta x$ in the $x$ direction and size $\Delta y$ in the $y$ direction to obtain the discrete convolution kernel $K^{p \to u}$:

$$K^{p \to u}(x_m - x_n, y_q - y_r) = \frac{1}{\pi E^*} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \int_{-\frac{\Delta y}{2}}^{\frac{\Delta y}{2}} \frac{1}{\sqrt{(x_m - x_n - x'')^2 + (y_q - y_r - y'')^2}} dx'' dy'' \qquad (A.4)$$

With the equivalent elastic modulus $E^*$ defined in equation (6).

The analytical expression of the discrete convolution kernel $K^{p \to u}$ is [29]:

$$K^{p \to u}(x_m - x_n, y_q - y_r) = \frac{1}{\pi E^*} \left[ \begin{aligned} & x_P \ln \left( y_P + \sqrt{x_P^2 + y_P^2} \right) + y_P \ln \left( x_P + \sqrt{x_P^2 + y_P^2} \right) \\ & + x_M \ln \left( y_M + \sqrt{x_M^2 + y_M^2} \right) + y_M \ln \left( x_M + \sqrt{x_M^2 + y_M^2} \right) \\ & - x_P \ln \left( y_M + \sqrt{x_P^2 + y_M^2} \right) - y_M \ln \left( x_P + \sqrt{x_P^2 + y_M^2} \right) \\ & - x_M \ln \left( y_P + \sqrt{x_M^2 + y_P^2} \right) - y_P \ln \left( x_M + \sqrt{x_M^2 + y_P^2} \right) \end{aligned} \right] \qquad (A.5)$$

With:

$$\begin{aligned} x_P &= x_m - x_n + \Delta x/2 \\ y_P &= y_q - y_r + \Delta y/2 \\ x_M &= x_m - x_n - \Delta x/2 \\ y_M &= y_q - y_r - \Delta y/2 \end{aligned} \qquad (A.6)$$

## Appendix B: Stress in an elastic half-space subjected to a surface pressure

The homogeneous stress $\sigma_{ij}$ is obtained from the contact pressure $p$ by convolution with Green's functions $G_{ij}^{p \to \sigma}$ [27,29]:

$$\sigma_{ij}(x, y, z) = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} p(x', y') G_{ij}^{p \to \sigma}(x - x', y - y', z) dx' dy' \qquad (B.1)$$
Green's function $G_{ij}^{p \to \sigma}$ is defined by the following expressions:

$$
\begin{align*}
G_{xx}^{p \to \sigma}(x, y, z) &= \frac{1}{2\pi} \left( \frac{1-2\nu}{r^2} \left( \left(1-\frac{z}{\rho}\right) \frac{x^2-y^2}{r^2} + \frac{zy^2}{\rho^3} \right) - \frac{3zx^2}{\rho^5} \right) \\
G_{yy}^{p \to \sigma}(x, y, z) &= \frac{1}{2\pi} \left( \frac{1-2\nu}{r^2} \left( \left(1-\frac{z}{\rho}\right) \frac{y^2-x^2}{r^2} + \frac{zx^2}{\rho^3} \right) - \frac{3zy^2}{\rho^5} \right) \\
G_{zz}^{p \to \sigma}(x, y, z) &= -\frac{3z^3}{2\pi\rho^5} \\
G_{yz}^{p \to \sigma}(x, y, z) &= -\frac{3yz^2}{2\pi\rho^5} \\
G_{xz}^{p \to \sigma}(x, y, z) &= -\frac{3xz^2}{2\pi\rho^5} \\
G_{xy}^{p \to \sigma}(x, y, z) &= \frac{1}{2\pi} \left( \frac{1-2\nu}{r^2} \left( \left(1-\frac{z}{\rho}\right) \frac{xy}{r^2} - \frac{xyz}{\rho^3} \right) - \frac{3xyz}{\rho^5} \right)
\end{align*}
\tag{B.2}
$$

With Poisson's ratio $\nu$ and:

$$
\begin{align}
r &= \sqrt{x^2 + y^2} \notag \\
\rho &= \sqrt{x^2 + y^2 + z^2} \tag{B.3}
\end{align}
$$

In a discrete domain with coordinates $(x_m, y_q, z_s)$ where the fields $\sigma_{ij}$ and $p$ are constants in each element of size $\Delta x$ in $x$ direction and size $\Delta y$ in $y$ direction:

$$
\sigma_{ij}(x_m, y_q, z_s) = \sum_n \sum_r K_{ij}^{p \to \sigma}(x_m - x_n, y_q - y_r, z_s) p(x_n, y_r) \quad (B.4)
$$

With the kernel $K_{ij}^{p \to \sigma}$ defined by:

$$
K_{ij}^{p \to \sigma}(x_m - x_n, y_q - y_r, z_s) = \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \int_{-\frac{\Delta y}{2}}^{\frac{\Delta y}{2}} G_{ij}^{p \to \sigma}(x_m - x_n - x'', y_q - y_r - y'', z) dx'' dy'' \quad (B.5)
$$

The components of the kernel $K_{ij}^{p \to \sigma}$ can be expressed analytically:

$$
\begin{equation}
\begin{split}
K_{ij}^{p \to \sigma}(x, y, z) ={}& S_{ij}^{p \to \sigma}\left(x + \frac{\Delta x}{2}, y + \frac{\Delta y}{2}, z\right) + S_{ij}^{p \to \sigma}\left(x - \frac{\Delta x}{2}, y - \frac{\Delta y}{2}, z\right) \\
& - S_{ij}^{p \to \sigma}\left(x + \frac{\Delta x}{2}, y - \frac{\Delta y}{2}, z\right) - S_{ij}^{p \to \sigma}\left(x - \frac{\Delta x}{2}, y + \frac{\Delta y}{2}, z\right)
\end{split}
\tag{B.6}
\end{equation}
$$

With:

$$
\begin{align*}
S_{xx}^{p \to \sigma}(x, y, z) &= \frac{\nu}{\pi} \arctan \left( \frac{z^2 + y^2 - y\rho}{zx} \right) + \frac{1-2\nu}{\pi} \arctan \left( \frac{\rho - y + z}{x} \right) + \frac{xyz}{2\pi\rho(x^2 + z^2)} \\
S_{yy}^{p \to \sigma}(x, y, z) &= \frac{\nu}{\pi} \arctan \left( \frac{z^2 + y^2 - y\rho}{zx} \right) + \frac{1-2\nu}{\pi} \arctan \left( \frac{\rho - x + z}{y} \right) + \frac{xyz}{2\pi\rho(y^2 + z^2)} \\
S_{zz}^{p \to \sigma}(x, y, z) &= \frac{1}{2\pi} \arctan \left( \frac{z^2 + y^2 - y\rho}{zx} \right) - \frac{xyz}{2\pi\rho} \left( \frac{1}{x^2 + z^2} + \frac{1}{y^2 + z^2} \right) \tag{B.7} \\
S_{yz}^{p \to \sigma}(x, y, z) &= \frac{xz^2}{2\pi\rho(y^2 + z^2)} \\
S_{xz}^{p \to \sigma}(x, y, z) &= \frac{yz^2}{2\pi\rho(x^2 + z^2)} \\
S_{xy}^{p \to \sigma}(x, y, z) &= -\frac{z}{2\pi\rho} - \frac{1-2\nu}{2\pi} \ln(\rho + z)
\end{align*}
$$