© L. Fourel et al., Published by EDP Sciences, 2025

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.

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.

Fig. 1 Analysis domain.

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 σ^M from contact pressure.

3. Computation of heterogeneous stress σ 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.

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$$(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$$(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)={\displaystyle {\int }_{-\infty }^{+\infty }p}\left(x^{\prime }\right)G\left(x-x^{\prime }\right)\text{d}{x}^{\prime }$$(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]: $$ℱ(u):\xi \mapsto \widehat{u}(\xi )={\displaystyle {\int }_{-\infty }^{+\infty }u}(x)e^{-i\xi x} \text{d}x$$(4)

The convolution of functions expressed in the space domain becomes an algebraic product in the frequency domain: $$\begin{array}{c}ℱ(p\ast G)=ℱ(p)ℱ(G)\\ \widehat{u}(\xi )=\widehat{p}(\xi )\widehat{G}(\xi )\end{array}$$(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 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.

Fig. 2. Computation steps.

Fig. 3. Initial gap h0, gap h, displacement u, approach δ.

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₁ and a Poisson’s ratio ν₁ and the second body has a Young’s modulus E₂ and a Poisson’s ratio ν₂. 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}$$(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 δ occurs. The surfaces deform to prevent interpenetration. Figure 3 illustrates the relationship between the gap h, the initial gap h₀, the normal displacement u and the approach δ: $$h(x,y)=h_0(x,y)+u(x,y)-\delta$$(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{array}{l}h(x,y)>0\to p(x,y)=0\hfill \\ h(x,y)=0\to p(x,y)>0\hfill \end{array}$$(8)

According to the force equilibrium, the integral of the contact pressure on the surface equals the load W: $$W={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }p}}(x,y)\text{d}x \text{d}y$$(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={\displaystyle {\int }_{{\text{Ω}}_C}\left(\frac{1}{2}pu+p{h}_0\right)}\text{d}S$$(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)$$(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→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)$$(12) $$G^{p\to u}(x,y)=\frac{1}{\pi E^{*}}\frac{1}{\sqrt{x^2+y^2}}$$(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 σ^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)$$(14)

The relationship between the concentrated contact pressure p located in (x = 0, y = 0) and the homogeneous stress σ^M is given by the Green’s function G^p→σ: $${\sigma }_{ij}^M(x,y,z)=p(x,y)*G_{ij}^{p\to \sigma }(x,y,z)$$(15)

Love [29] and Kalker [27] derived the Green’s function $G_{ij}^{p\to \sigma }$ and the discrete convolution kernel $K_{ij}^{p\to \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: $$\{\begin{array}{l}{\widehat{K}}_{ij}^{p\to \sigma }=\mathrm{FFT}\left(K_{ij}^{p\to \sigma }\right)\hfill \\ \widehat{p}=\mathrm{FFT}(p)\hfill \\ {\widehat{\sigma }}_{ij}^M=\widehat{p}{\widehat{K}}_{ij}^{p\to \sigma }\hfill \\ {\sigma }_{ij}^M=\mathrm{iFFT}\left({\widehat{\sigma }}_{ij}\right)\hfill \end{array}$$(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 ε is considered to be the sum of the homogeneous strain ε^M and the eigenstrain also called fluctuation strain $\tilde{\epsilon }$: $${\epsilon }_{ij}={\epsilon }_{ij}^M+{\tilde{\epsilon }}_{ij}$$(17)

${\epsilon }_{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 ς is defined by: $${\varsigma }_{ij}=\text{Δ}{C}_{ijkl}{\epsilon }_{kl}$$(18)

With: $$\text{Δ}{C}_{ijkl}=C_{ijkl}-C_{ijkl}^0$$(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⁰ that is constant for all points in the domain.

The fluctuation strain $\widehat{\tilde{\epsilon }}$ is computed in the frequency domain with coordinates ξ from the polarization stress $\widehat{Ϛ}$ using: $${\widehat{\tilde{\epsilon }}}_{ij}=\frac{1}{2}\left({\widehat{\text{Γ}}}_{ijkl}^{\varsigma \to \tilde{\epsilon }}+{\widehat{\text{Γ}}}_{jikl}^{\varsigma \to \tilde{\epsilon }}\right){\widehat{\varsigma }}_{kl}$$(20)

With the Green’s operator ${\widehat{\Gamma }}_{ijkl}$: $${\widehat{\text{Γ}}}_{ijkl}={\widehat{G}}_{ik,lj}=-{\widehat{G}}_{ik}{\xi }_j{\xi }_l$$(21)

With the Green’s function ${\widehat{G}}_{ik}$: $${\widehat{G}}_{ik}={\left(C_{ijkl}^0{\xi }_j{\xi }_l\right)}^{-1}$$(22)

The fixed-point computation scheme proposed by Moulinec and Suquet [13,14] is used: $$\begin{array}{l}\text{ Initialization: }{\epsilon }_{ij}^M={\left(C_{ijkl}^0\right)}^{-1}{\sigma }_{kl}^M\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_{ij}^{(0)}={\epsilon }_{ij}^M\\ \text{ Iteration: }{\varsigma }_{ij}^{(n)}=\text{Δ}{C}_{ijkl}{\epsilon }_{kl}^{(n-1)}\\ \begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\widehat{\varsigma }}_{ij}^{(n)}=FFT\left({\varsigma }_{ij}^{(n)}\right)\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\widehat{\hspace{0.17em}\tilde{\epsilon }}}_{ij}^{(n)}=\frac{1}{2}\left({\widehat{\text{Γ}}}_{ijkl}^{\varsigma \to \tilde{\epsilon }}+{\widehat{\text{Γ}}}_{jikl}^{\varsigma \to \tilde{\epsilon }}\right){\widehat{\varsigma }}_{kl}^{(n)}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_{ij}^{(n)}={\epsilon }_{ij}^M+iFFT\left({\widehat{\tilde{\epsilon }}}_{ij}^{(n)}\right)\hfill \end{array}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{Convergence test}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{Result:}{\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}^{(n)}\end{array}$$(23)

With the following convergence test: $$\frac{\langle \Vert {\epsilon }_{ij}^{(n)}-{\epsilon }_{ij}^{(n-1)}\Vert \rangle }{\langle \Vert {\epsilon }_{ij}^M\Vert \rangle }<e$$(24)

With the L2 norm $\Vert \cdot \Vert$, the spatial arithmetic mean $\langle \cdot \rangle$ 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⁻³. 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₀ at the center. The lengths are made dimensionless using the contact half-width a and the stress distributions using the Hertz pressure p₀ . 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 × 256 × 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

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=\left[\begin{array}{cccccc}{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{array}\right]$$(25)

With the cubic elastic properties C₁₁, C₁₂ and C₄₄. 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: $$\begin{array}{l}{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)\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times R_{k_1{k}_2}(x,y,z)R_{l_1{l}_2}(x,y,z)C_{i_2{j}_2{k}_2{l}_2}^c\hfill \end{array}$$(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 σ_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 σ_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 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|$$(27)

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 Lb 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 tha 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.

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.

Fig. 5. Heterogeneous normal stress σzz in surface elements at the center of the point contact.

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 σ_zz = -p is imposed and the other components are set to zero. First-order polynomial elements are used.

The difference ∆σ 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₀ : $$\text{Δ}{\sigma }_{ij}=\left|\frac{{{\sigma }_{ij}|}_{\text{Green-FFT }}-{{\sigma }_{ij}|}_{\text{FEM}}}{p_0}\right|$$(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.

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 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 3s 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.

Fig. 6 (a) Periodization error ezz without buffer zone Lb = 0 for a grain size DG = 0.080 a, (b) maximum periodization error depending on the buffer size Lb relative to the grain size DG.

Fig. 7. Equivalent FEM model.

Fig. 8. Comparison between the Green-FFT model and the equivalent FEM model.

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₀, respectively.

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

$$r(x,y)=\sqrt{x^2+y^2}$$(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.

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

Fig. 9. Execution time of the Green-FFT model and the equivalent FEM model depending on the degrees of freedom (DOF) of the models.

Fig. 10. Maximum shear stress with a subsurface defect, (a) inclusion, (b) porosity.

Fig. 11 Dent profile.

Fig. 12. Maximum shear stress with a dented surface.

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.

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

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\left(x_m\right),\forall x\in \left[x_m-\frac{\text{Δ}x}{2},x_m+\frac{\text{Δ}x}{2}\right]$$(A.1)

With the element size Δx. The discrete convolution is then defined as a sum: $$u\left(x_m\right)={\displaystyle \sum _{n}p}\left(x_n\right)K\left(x_m-x_n\right)$$(A.2)

With the discrete convolution kernel K: $$K\left(x_m-x_n\right)={\displaystyle {\int }_{-\frac{\text{Δ}x}{2}}^{\frac{\text{Δ}x}{2}}G}\left(x_m-x_n-x^{\prime }\right)d{x}^{\prime }$$(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→u defined in equation (13) is integrated over each rectangular element of size Δx in the x direction and size Δy in the y direction to obtain the discrete convolution kernel K^p→u: $$K^{p\to u}\left(x_m-x_n,y_q-y_r\right)=\frac{1}{\pi E^{*}}{\displaystyle {\int }_{-\frac{\text{Δ}x}{2}}^{\frac{\text{Δ}x}{2}}{\displaystyle {\int }_{-\frac{\text{Δ}y}{2}}^{\frac{\text{Δ}y}{2}}\frac{1}{\sqrt{{\left(x_m-x_n-x^{\prime \prime }\right)}^2+{\left(y_q-y_r-y^{\prime \prime }\right)}^2}}}}d{x}^{\prime \prime }d{y}^{\prime \prime }$$(A.4)

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

The analytical expression of the discrete convolution kernel K^p→u is [29]: $$\begin{array}{l}{K}^{p\to u}\left(x_m-x_n,y_q-y_r\right)=\frac{1}{\pi E^{*}}\\ \begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}[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)\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+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)\hfill \end{array}\\ \begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-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)\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-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)]\hfill \end{array}\end{array}$$(A.5)

With: $$\begin{array}{c}{x}_P=x_m-x_n+\text{Δ}x/2\\ y_P=y_q-y_r+\text{Δ}y/2\\ x_M=x_m-x_n-\text{Δ}x/2\\ y_M=y_q-y_r-\text{Δ}y/2\end{array}$$(A.6)

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

The homogeneous stress σ_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)={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }p}}\left(x^{\prime },y^{\prime }\right)G_{ij}^{p\to \sigma }\left(x-x^{\prime },y-y^{\prime },z\right)d{x}^{\prime }d{y}^{\prime }$$(B.1)

Green’s function $G_{ij}^{p\to \sigma }$ is defined by the following expressions: $$\begin{array}{l}\begin{array}{c}{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{z{y}^2}{{\rho }^3}\right)-\frac{3z{x}^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{z{x}^2}{{\rho }^3}\right)-\frac{3z{y}^2}{{\rho }^5}\right)\end{array}\\ \begin{array}{c}{G}_{zz}^{p\to \sigma }(x,y,z)=-\frac{3z^3}{2\pi {\rho }^5}\\ G_{yz}^{p\to \sigma }(x,y,z)=-\frac{3y{z}^2}{2\pi {\rho }^5}\\ G_{xz}^{p\to \sigma }(x,y,z)=-\frac{3x{z}^2}{2\pi {\rho }^5}\end{array}\\ 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{array}$$(B.2)

With Poisson’s ratio ν and: $$\begin{array}{l}r=\sqrt{x^2+y^2}\hfill \\ \rho =\sqrt{x^2+y^2+z^2}\hfill \end{array}$$(B.3)

In a discrete domain with coordinates (x_m, y_q, z_s) where the fields σ_ij and p are constants in each element of size Δx in x direction and size Δy in y direction: $${\sigma }_{ij}\left(x_m,y_q,z_s\right)={\displaystyle \sum _n{\displaystyle \sum _r{K}_{ij}^{p\to \sigma }}}\left(x_m-x_n,y_q-y_r,z_s\right)p\left(x_n,y_r\right)$$(B.4)

With the kernel $K_{ij}^{p\to \sigma }$ defined by: $$K_{ij}^{p\to \sigma }\left(x_m-x_n,y_q-y_r,z_s\right)={\displaystyle {\int }_{-\frac{\text{Δ}x}{2}}^{\frac{\text{Δ}x}{2}}{\displaystyle {\int }_{-\frac{\text{Δ}y}{2}}^{\frac{\text{Δ}y}{2}}{G}_{ij}^{p\to \sigma }}}\left(x_m-x_n-x^{\prime \prime },y_q-y_r-y^{\prime \prime },z\right)d{x}^{\prime \prime }d{y}^{\prime \prime }$$(B.5)

The components of the kernel $K_{ij}^{p\to \sigma }$ can be expressed analytically: $$\begin{array}{l}{K}_{ij}^{p\to \sigma }(x,y,z)=S_{ij}^{p\to \sigma }\left(x+\frac{\text{Δ}x}{2},y+\frac{\text{Δ}y}{2},z\right)+S_{ij}^{p\to \sigma }\left(x-\frac{\text{Δ}x}{2},y-\frac{\text{Δ}y}{2},z\right)\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-S_{ij}^{p\to \sigma }\left(x+\frac{\text{Δ}x}{2},y-\frac{\text{Δ}y}{2},z\right)-S_{ij}^{p\to \sigma }\left(x-\frac{\text{Δ}x}{2},y+\frac{\text{Δ}y}{2},z\right)\hfill \end{array}$$(B.6)

With: $$\begin{array}{l}\begin{array}{c}{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 \left(x^2+z^2\right)}\\ 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 \left(y^2+z^2\right)}\end{array}\\ \begin{array}{l}{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)\hfill \\ \begin{array}{l}{S}_{yz}^{p\to \sigma }(x,y,z)=\frac{x{z}^2}{2\pi \rho \left(y^2+z^2\right)}\\ \begin{array}{l}{S}_{xz}^{p\to \sigma }(x,y,z)=\frac{y{z}^2}{2\pi \rho \left(x^2+z^2\right)}\hfill \\ S_{xy}^{p\to \sigma }(x,y,z)=-\frac{z}{2\pi \rho }-\frac{1-2\nu }{2\pi }\ln (\rho +z)\hfill \end{array}\end{array}\hfill \end{array}\end{array}$$(B.7)

All Tables

All Figures

	Fig. 1 Analysis domain.
In the text

	Fig. 2. Computation steps.
In the text

	Fig. 3. Initial gap h0, gap h, displacement u, approach δ.
In the text

	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.
In the text

	Fig. 5. Heterogeneous normal stress σzz in surface elements at the center of the point contact.
In the text

	Fig. 6 (a) Periodization error ezz without buffer zone Lb = 0 for a grain size DG = 0.080 a, (b) maximum periodization error depending on the buffer size Lb relative to the grain size DG.
In the text

	Fig. 7. Equivalent FEM model.
In the text

	Fig. 8. Comparison between the Green-FFT model and the equivalent FEM model.
In the text

	Fig. 9. Execution time of the Green-FFT model and the equivalent FEM model depending on the degrees of freedom (DOF) of the models.
In the text

	Fig. 10. Maximum shear stress with a subsurface defect, (a) inclusion, (b) porosity.
In the text

	Fig. 11 Dent profile.
In the text

	Fig. 12. Maximum shear stress with a dented surface.
In the text