FFT-Based Methods for Computational Contact Mechanics

Computational contact mechanics seeks for numerical solutions to contact area, pressure, deformation, and stresses, as well as flash temperature, in response to the interaction of two bodies. The materials of the bodies may be homogeneous or inhomogeneous, isotropic or anisotropic, layered or functionally graded, elastic, elastoplastic, or viscoelastic, and the physical interactions may be subjected to a single field or multiple fields. The contact geometry can be cylindrical, point (circular or elliptical), or nominally flat-to-flat. With reasonable simplifications, the mathematical nature of the relationship between a surface excitation and a body response for an elastic contact problem is either in the form of a convolution or correlation, making it possible to formulate and solve the contact problem by means of an efficient Fourier-transform algorithm. Green's function inside such a convolution or correlation form is the fundamental solution to an elementary problem, and if explicitly available, it can be integrated over a region, or an element, to obtain influence coefficients (ICs). Either the problem itself or Green's functions/ICs can be transformed into a space-related frequency domain, via a Fourier transform algorithm, to formulate a frequency-domain solution for contact problems. This approach converts the original tedious integration operation into multiplication accompanied by Fourier and inverse Fourier transforms, and thus a great computational efficiency is achieved. The conversion between ICs and frequency-response functions facilitate the solutions to problems with no explicit space-domain Green's function. This paper summarizes different algorithms involving the fast Fourier transform (FFT), developed for different contact problems, error control, as well as solutions to the problems involving different contact geometries, different types of materials, and different physical issues. The related works suggest that (i) a proper FFT algorithm should be used for each of the cylindrical, point, and nominally flat-flat contact problems, and then (ii) the FFT-based algorithms are accurate and efficient. In most cases, the ICs from the 0-order shape function can be applied to achieve satisfactory accuracy and efficiency if (i) is guaranteed.


INTRODUCTION
Contact of materials is a common engineering phenomenon, and the solution to a contact problem, in terms of pressure, deformations, and stresses, as well as flash temperature, is usually among first steps in the design and analysis of an engineering system or a functional device. A contact problem is solved first for the information of the contacting interface, such as contact pressure, surface interaction, contact area, and interfacial friction, followed by a boundary-value solution process for the stresses in each contacting body. At least two sets of convolutions, each for a group of excitations and Green's functions, are involved, i.e., displacements and stresses, in response to surface tractions, and the former are calculated first (Conry and Seireg, 1971;Kalker and van Randen, 1972;Kalker, 1986;Polonsky and Keer, 1999). When solving the contact of bodies involving an inclusioncontaining material, the mathematical correlation between an eigenstrain and a Green's function appears . In most cases, the convolution nature makes it possible to formulate and solve a contact problem by means of an efficient fast Fourier transform (FFT) algorithm (Ju and Farris, 1996;Stanley and Kato, 1997;Ai and Sawamiphakdi, 1999;Hu et al., 1999;Keer, 1999, 2000;Liu et al., 2000;. We will discuss and summarize the theories, algorithms, and numerical methods of FFT-based contact modeling approaches. If the contact body can be treated as a mathematically "semiinfinite" medium, or a half-space material, analytical influence coefficients (ICs), or frequency-response functions (FRFs), may be derived for certain problems. Normally such a simplification can be made if deformations are small and the radii of curvature of contact bodies are sufficiently larger than the effective contact size of these two bodies. In the case of the macroscopic contact of two spheres, or two cylinders, of the same size and material, if the contact radius, or half width, a, is 10% of the radius of the contact body, the maximum error from using the z-direction deformation to replace the radial deformation at the contact edge is within 0.5%. However, stresses need more attention. Two contact problems, (1) two steel spheres of the same radius and (2) their equivalent sphere and a half space, are analyzed using the finite-element method (FEM) without considering plasticity. Figure 1 shows maximum relative errors of the von Mises stress along the central axis, due to ignoring curvature, as a function of a/R e . The error is about 7% when a/R e reaching 10%, and this situation is beyond the "small elastic deformation" assumption. In the elastic range, using steels as an example, the maximum error of the von Mises stress is <3% if the half-space solution is pursued. Similar small errors due to the use of the halfspace assumption in the elastic range have also been reported by Londhe et al. (2018), in the comparison of the results from the FEM and Hertz formulas for different types of contacts.
In a contact simulation, the computational complexity of evaluating a convolution via direct summation of the products of ICs and surface traction is on the order of O(N 2 ), where N is the mesh number. If N is large and the convolution has to be repeatedly calculated in an iteration process, the computation burden is very heavy. The works by Ju and Farris (1996), Nogi and Kato (1997), Stanley and Kato (1997), Ai and Sawamiphakdi (1999), Hu et al. (1999), Polonsky and Keer (2000), and Liu et al. (2000) are in a chain of studies to apply the FFT to evaluate the convolutions for elastic deformation and stresses efficiently in the field of contact mechanics and tribology. Several papers have reviewed such efforts of solving tribology problems via the FFT and analyzed the sources of errors (Liu et al., 2000;Wang et al., 2003;Wang and Zhu, 2019). Although most of the effort is on solving non-conformal contact problems, the FIGURE 1 | Maximum errors of the von Mises stress, σ VM , along the z axis, or the central axis in the depth direction, by using the half-space approach to solve the problem of contact of two equal spheres, without considering plasticity, calculated with the FEM; a is the half contact width, R e is the equivalent radii of the contact bodies, and 1/R e = 1/R 1 + 1/R 2 .
FFT methods suit for certain conformal-contact problems as well if they involve a convolution and have ICs obtained analytically or numerically. Actually, the circular nature of a cylindrical structure fits the circular convolution theorem perfectly. Liu and Chen (2012) and Liu (2013) reported an FFT-based conformalcontact model for two-dimensional (2D) problems. Wang and Jin (2004) conducted the fluid-film lubrication analysis for artificial joints, which requires the determination of the elastic deformation of the bearing surface of both the acetabular cup and the femoral head. They used FFT along with the spherical distance and numerical ICs from an FEM calculation.
The FFT methods greatly help reduce the computation burden. For example, for a three-dimensional (3D) point-contact problems, the FFT operation is on the order of O(12Nlog 2 [4N]), to be discussed in detail later. Its ratio to the operation needed for calculating the convolution is 12Nlog 2 (4N)/N 2 = 12log 2 (4N)/N = 0.00024, if N is 1024 * 1024. This is a significant saving of computational time. The key issues to be addressed are (1) how should the Fourier transform (FT) method be properly used to solve contact problems? (2) Can one solution algorithm be used for all problems?
Theoretically, the Fourier transform can be used for infinitedomain and the Fourier series for periodic problems, but most contact problems do not satisfy these conditions. For example, a point-contact problem has its pressure only on a small region of contacting surfaces. If the FFT is directly used to solve such a problem, the results near the borders have notable errors. In order to reduce the periodicity error, Ju and Farris (1996) substantially extended the domain, Ai and Sawamiphakdi (1999) decomposed the total pressure into a smooth portion and a zero-mean fluctuating portion, and Polonsky and Keer (2000) developed a hybrid algorithm by adding a special correction procedure. Each of these brings in a certain accuracy improvement while introducing new complications. We have investigated the theories of contact mechanics and signal analyses, and realized that proper convolution theorems should be considered in solving different contact problems of nonconformal and conformal configurations (Figure 2) (Liu et al., 2000;Liu and Hua, 2009;Liu and Chen, 2012), as summarized by Wang and Zhu (2019, Chapter 4).
In this paper, we will discuss the basic issues of the FFT methods for contact analyses from the convolution theorems and the tree of the Fourier-transform algorithms for solving different contact problems, such as (1) the algorithm of discrete-convolution and fast-Fourier-transform (DC-FFT), with double domain extension in each dimension, for non-periodic problems, and the discrete-convolution and fast-Fourier-transform algorithm (DC-FFT) without domain extension for journal bearing problems, (2) the algorithm of continuous-convolution and Fourier-transform (CC-FT) for periodic (or infinite) contact problems, (3) the algorithm of discrete convolution with duplicated padding and FFT (DCD-FFT), that of discrete-continuous convolutions and FFT (DC-CC-FFT), and that of the discrete convolution with IC summation and FFT (DCS-FFT) for 3D line-contact problems that are periodic (or infinite) in one direction but non-periodic in the other direction, (4) the algorithm of discrete-correlation and FFT (DCR-FFT) for inclusion problems, (5) the FRF-IC conversion method, as well as the applications of them to solve the contact problems involving layered materials, anisotropic elastic materials, and viscoelastic polymers, or those subjected to multifields, and (6) a non-uniform DCS-FFT method, recently developed by Sun et al. (2020), for solve large-scale contact problems. Most of the contents in this paper are based on the works by Liu et al. (2000, , Boucly et al. (2005), Chen et al. (2008), Liu and Hua (2009), Yu et al. (2014, Zhang X. et al. (2017Zhang X. et al. ( , 2018, , Sun (2020), Sun et al. (2020), Zhang et al. (2020a,b), and Chapter 4 in the book by Wang and Zhu (2019). Table 1 summarizes the FFF-based approaches and their applications.

Name and references
Algorithm and method Problem to solve DC-FFT (Liu et al., 2000;Liu, 2001) (Liu and Chen, 2012;Liu, 2013) Discrete-convolution and FFT Point-contact problems Cylindrical contact problems, counterformal and conformal CC-FT (Ju and Farris, 1996;Liu, 2001; DCSS-FFT (Sun, 2020) Continuous-convolution and FFT Nominally flat-flat contact problems DCD-FFT (Chen et al., 2008) DC-CC-FFT (Liu et al., 2006) DCS-FFT (Liu and Hua, 2009;Sun et al., 2020) Discrete convolution with duplicated padding and FFT Discrete-convolution, continuous-convolution and FFT Discrete convolution with IC summation and FFT 3D line-contact problems DCR-FFT  Discrete-convolution-correlation and fast-Fourier-transform Materials with residual strains, inclusions and/or inhomogeneities, contact plasticity problems IC conversion (Liu et al., 2000;Liu, 2001; FRFs known, which can be transformed to ICs, followed by the DC-FFT algorithm or other proper ones Layered, viscoelastic, transversely isotropic materials, coupled-stress problems, multifield contact problems It should be mentioned that the pressure and contact area that satisfy the Kuhn-Tucker type complementary conditions can be solved with different methods, and the conjugate gradient method (CGM) (Polonsky and Keer, 1999;Jin et al., 2013) is currently the widely accepted one. In the following, the elastic field means the distributions of stresses and displacements, and the target domain means the physical domain, on which a physical contact problem is defined. The algorithms and methods will be explained mainly through deformation calculations; details of the FFT-based computations of stresses, flash temperature, and other physical fields can be found in the reports by Liu (2001), , Chen et al. (2008), Zhang X. et al. (2018), and Wang and Zhu (2019), as well as those mentioned in the previous paragraph.

Convolution and ICs
Let's use the pressure-displacement relationships, such as the Flamant and Boussinesq equations (Johnson, 1987), as examples. Here, an excitation at ξ , or (ξ , η), and a response at x, or (x, y), are related to each other through a Green's function defined with the distance between the two, which is either |x − ξ | or R I = The surface normal displacement of a cylinder in a line contact, u z (x), due to pressure p(x) on surface region S x is (1) where C = −4/(πE ′ ), E' is the effective Young's modulus and the corresponding Green's function is C ln |x|. The integral kernel, G, is defined as In numerical modeling, the equation above can be discretized and re-written as the summation of the products of influence coefficients D(k, i), or D i,j , and nodal pressures p i .
where N p , N d are the total numbers of nodes for pressure and deformation, respectively. A shape function, Ys, may be used to distribute pressure, or other excitations, around a nodal point, and the commonly used shape functions are 0-order, 1st-order, and 2nd-order polynomials. Detailed expressions and use of these shape functions can be found in the book by Wang and Zhu (2019).
The ICs are from the elementary integration of Green's function and shape function Ys, implying unit nodal pressure, or from G and shape function Ys without including coefficient C. The latter is used here. In general, where f is the integration result, and x u and x l are the upper and lower boundaries of the element integration at x i , with 1 and If the zero-order shape function is used, the influence coefficient expression becomes Note that here the influence coefficients, D, depend only on the geometric factors of the grid. When a uniform grid is used with the constant mesh spacing x i+1 -x i = x i -x i−1 = 2 , the above becomes For point-contact problems, surface normal displacement u z (x, y) due to pressure p(x, y) on surface area is, x 2 +y 2 . Likewise, the equation above can be re-written as the summation of the products of influence coefficients D k,l i,j and nodal pressures p i,j . The ICs are from the elementary integration of Green's function and a shape function, Ys, or G and shape function Ys without including the material properties, and the latter is used here.
where N dx , N dy are the total numbers of nodes for deformation in the x, y directions, and N px , N py are the total numbers of nodes for pressure in the x, y directions, respectively.
where f is the integration result, and (x u , y u ) (x l , y l ) are the upper and lower boundaries of the element integration.
If the zero-order shape function is used, the influence coefficient expression is (Love, 1929), Replacing x k−i = x k − x i and y l−j = y l − y j leads to the following.
where a and b are the half length of the rectangular integration element.

Frequency Response Functions and ICs
The Fourier transform can be applied to the IC matrix, D, Equation (4), to obtain,D

=G ·Ỹs
andG =D/Ỹs whereG is the frequency response function excluding the elastic parameter, C. These two equations show how to obtain one from the other. If discrete Fourier transformD has already been obtained from a set of known ICs via the FFT,G can be solved fromG = D/Ỹs once Fourier series coefficientsD can be obtained from D. Based on the sampling theorem, the one-dimensional (1D) relationship between the FT (∼) and discrete Fourier transform (DFT) ( ∧ ) of the ICs, sampled with mesh interval 2 , can be obtained aŝ Frontiers in Mechanical Engineering | www.frontiersin.org An aliasing control parameter, AL, can be introduced instead of summation for r from -∞ to ∞, in order to satisfy a required accuracy while saving computation time.
The above can be simplified if the sampling frequency is sufficiently high, or the datum interval is sufficiently small.
The above becomes exact if and only if m s > 2m max (Morrison, 1994) with ±m max as the band limit (or m max as the highest frequency component), beyond which there is no Fourier transform results (Walker, 1996). Then, only the term at r = 0 is needed. Likewise, the 2D relationship between the FT (∼) and DFT ( ∧ ) of the ICs iŝ where x and y are the mesh intervals in the x and y directions, and x = y = 2 if the meshes are uniform in both directions. This means two operations.
(2) Upon knowingG, we can getD fromD using Equations (13) and (17), or (18). In many cases where the solutions to frequency response functions are more convenient than those to Green's functions, and this operation can be utilized to convert the FRFs to the discrete Fourier transformed ICs.

CONVOLUTION THEOREMS
Equations (1) and (7) are both convolutions, which can be solved efficiently via Fourier transform followed by inverse Fourier transform. However, because the pressure may be in a discrete form, i.e., rough-surface contact pressure, and its application domain may be in different sizes, i.e., finite or infinite, accurate solutions to these equations, and others in the same nature, require the use of different convolution theorems. The following explains these theorems with 1D datum series for convenience.

Continuous Linear Convolution
If a set of continuous functions of t, f and G, follows the convolution in Equation (20), resulting in O, then the Fourier transform of the convolution results,Õ(ω), satisfies Equation (21), which convert the integration in Equation (20) to multiplication of continuous Fourier transforms of f and G.
where O is accounted as the response of the continuous linear convolution, and symbol " * " is the convolution operator;Õ(ω),

Periodic Convolution
If D p (i) and f p (i) are periodic functions with period N, the product of their discrete Fourier series (DFS) coefficients, and ⌣ f p (m), with m for the frequency, is ⌣ O p that can be expressed as (Oppenheim et al., 1999), The periodic sequence, O p , with the same period, N, is the periodic convolution analyzed in one period, N, shown below.
For contact problems, D p (i) and f p (i) are D p (x i ) and f p (x i ). Here, the subscript can be negative.

Cyclic (Circular) Convolution
Based on Oppenheim et al. (1999), if D c (i) and f c (i) are sequences of finite length N, and their DFT results areD c (m) andf c (m), with m for the frequency, then their term-by-term product isÔ c , expressed asÔ The finite sequence O c is actually the circular convolution analyzed in length N, shown below.
where H(r-j) is the Heaviside unit step function, which is 1 when r-j is positive, or 0 otherwise. The term NH(r-j) contributes to the subscript numbering only when the step function is not zero to avoid negative subscript. Here, subscript -r reverses the sequence of the D series and subscript j shifts it in a circular fashion. More details has been given by Liu et al. (2000). This means that the Fourier transform operation in Equation (24) is valid only when the convolution of D c (i) and f c (i) is in the form given in Equations (25) or (26). Although Equations (23) and (25) are for different events, they are expressed in the same form and lead to the same results in one period (Oppenheim et al., 1999).
The cyclic (circular) convolution is for series of finite lengths, and the circular fashion of the D series make it suitable for problems with either periodic features (e.g., nominally flat contact with a special IC treatment) or a circular configuration (e.g., cylindrical and journal bearing problems). Liu and Chen (2012) and Liu (2013) reported an FFT-based conformalcontact model for 2D problems with two concentric cylindrical interfaces, Figure 1F, for which we can directly apply the cyclic convolution theorem and 1D FFT operations to obtain the shaft deformation. However, for problems without the periodic features, such as counterformal line/point contact problems, special measures are needed to make the D series in such a needed circular fashion so that the cyclic convolution can be properly performed. The ICs and the pressure series can be properly handled based on the characteristics of the problems, e.g., Liu S. et al. (2000, , and Chen et al. (2008), so that they fit the need for the circular-convolution analyses. Because the FFT is a collection of algorithms for fast execution of DFT, the cyclic convolution of the two datum series (Equation 26) should also satisfy Equation (24) when the DFT is replaced by the FFT.

FFT ALGORITHMS FOR CONTACT MECHANICS Cyclic Convolution and the DC-FFT Algorithm for Non-periodic Contact Problems
Consider the general line-contact displacement problem shown in Equation (3), subjected to a Hertzian pressure, or any localized pressure in a certain distribution, where influence coefficient D(k, i) means u z /C at node k caused by a unit pressure at node i, on a uniform grid of mesh spacing 2 . This is a problem of the convolution of two series of finite lengths; it is not infinite, nor periodic. Therefore, the cyclic convolution theorem, Equation (26), should be applied in order to solve it with the Fourier transform method, for which the IC matrix has to be a cyclic matrix. This section explains how such a matrix is constructed from the original IC matrix via wrap-around order, and how this problem is solved properly and efficiently via the FFT. The wrap-around order requires one-to-one extension of the target domain on which the physical problem is defined, and the pressure on the extended domain should be set to zero (zero padding).
It should be mentioned that if such a problem were solved with the continuous convolution theorem, Equation (22), via FRF and Fourier transform of pressure in the finite target domain, a noticeable error would appear at the borders because it periodizes the problem mathematically. Ju and Farris (1996) depressed this error by five times domain extension. The analysis by Liu et al. (2000) indicated that a complete error removal would require 16 times domain extension, as if the problem were infinite. Error analyses will be discussed in a later section.
Because influence coefficient D(k, i) only depends on the grid geometry, or, more specifically, the distance between points k and i, for a given uniform grid, it relies solely on k − i no matter what k or i is. We can define X k−i as the non-dimensional distance (normalized by a characteristic length, a) from X k to X i , or X k−i = X k to X i , then the IC component can be expressed as D k−i . Subscript k-i marks each element in the IC matrix. Obviously here for Equation (3) (3) can be expressed as follows in a non-dimensional form, For example, if a and p h are the Hertzian contact half width and maximum pressure, P(X) = √ 1 − X 2 , and the Hertzian pressure distribution is then p = P p h , and if the problem is digitized with The solution requires 5 × 5, or NxN, multiplication operations. The IC matrix above is a Toeplitz matrix, or a diagonalconstant matrix. This matrix has to be converted to a cyclic one in order to utilize the cyclic convolution theorem . This can be done with the operation of the wrap-around order (Bracewell, 1978;Brigham, 1988;Press et al., 1992) by adding the reversed first column without the first element, which is [D 4 , D 3 , D 2 , D 1 ], to the end of the first row, . Then, the extended first row becomes The total number of the extended series is N c = 2N−1 = 9. By using the Heaviside step function notation, j-r+N c H(r-j), with H(r−j) =1 for r-j >0, Equation (29) becomes Equation (30). Series D c in Equation (26) can be written as [< D >] N c =2N−1 to show the cyclic nature by < > and the size information by the subscript. We can also use D 2N−1 to express this extended wraparound matrix. A short vertical bar is used to separate the original and extended terms for clarity.
where EX means the extended components. The nodal pressure vector is also extended by zero padding as follows where "•" means the operation of term-by-term type complex multiplication. Note that the equal sign with arrows indicates that the vector at the right-hand side contains exactly that at the lefthand side, but the former has extra useless terms in the extension. This is the expression for the discrete convolution (cyclic convolution) and fast Fourier transform (DC-FFT) algorithm, named by Liu et al. (2000), for the deformation calculation. Similar expressions can be obtained for stress calculations (Liu, 2001;. Because the FFT operation of an N-number series is in the order of Nlog 2 N, the operation of Equation (32) is in the order of 3N c log 2 N c , much smaller than that of the direct summation (DS) operation, which is N×N, especially when N is large, as shown earlier in Introduction. Wang and Zhu (2019) offer a detailed numerical example, which shows that the direct summation, Equation (28), and the DFT-IDFT operation of the DC-FFT algorithm, Equation (32), lead to the same results in the analyzed accuracy. In order to show the cyclic nature of this operation, the fully extended matrix, or the cyclic matrix, is completely constructed from [< D >] N c = 2N−1 by circulating the last element in one row to the first position in the next row, given below.
Only the top portion of the matrix is written out because they are related to the physical target domain of the original problem in the matrix operation. When using the DC-FFT algorithm, only the first row is needed, and only the first N rows of the IFFT results are the needed solutions while the others should be discarded. IC wrap around order and pressure zero padding, shown in Figure 3, are two important operations for the DC-FFT algorithm to utilize the cyclic convolution and solve contact mechanics problems. Different implementation variations can be made; however, these two operations are necessary. The wraparound order for deformation calculation can be done by shifting the negative side or flip those between 1∼N-1. Caution should be paid for stress analyses because the shear-stress ICs are antisymmetric. Sometimes, ICs also need zero padding, which should be done at node N where the IC is the smallest. Analogous to the line-contact problem discussed above, the IC matrix for a 3D point-contact problem can be constructed by extending the original physical domain in the two lateral directions, as shown in Figure 4A. Equation (9) becomes the following for a 3D contact problem, where a is a reference length, which could be the Hertzian radius of a spherical contact, and p h a reference pressure, which could be the maximum Hertzian pressure.
where N px = N dx = N x and N py = N dy = N y are the numbers of nodes in the x and y directions, respectively. The extended circular convolution IC matrix is where D on the right-hand side means the ICs corresponding to the physical domain. It should be emphasized that two domains are involved, which are (1) the target domain, or the physical domain where the contact problem is defined, and (2) the extended domain, or the computation domain. At the end of the calculation, only the data for the target domain should be retained for the results. It is worth mentioning that the total lengths of P and D here are the minimum requirements; longer series than these should work, simply at a higher cost of computational efficiency. With these in mind, we can construct certain variations to implement this extension in order to apply the circular convolution for different situations, which are not discussed here.
The solution process involves IC matrix calculation, IC wraparound order, pressure zero padding, and the FFT-IFFT. Its procedure is detailed as follow (Liu et al., 2000).
-N x to N x − 1 in the x direction and -N y to N y − 1 in the y direction; 2) Apply the wrap-around order and zero padding to convert the IC matrix into a cyclic matrix, [< D >] 2N x ×2N y in the calculation domain, marked as 1 : 2Nx, 1 : 2Ny (or 0 : 2Nx − 1, 0 : 2Ny − 1 ), in the x and y directions, as shown in Figure 4B, and then apply the two-dimensional FFT to obtain the Fourier transformed IC matrix, < D > ; 3) Input pressure, conduct zero padding to convert pressure [P] N x ×N y into[P] 2N x ×2N y and then apply the two-dimensional FFT to get < P > ; 4) Obtain a temporary frequency series using element-byelement product of the two, i.e., < D > • < P > , where "•" means the operation of element complex multiplication; 5) Conduct two-dimensional IFFT < D > • < P > to obtain the surface deformation and keep the result data within the original physical domain.
The error analyses by Liu et al. (2000) and Wang et al. (2003) have convinced that (1) the DC-FFT algorithm generates no additional inaccuracy beyond the discretization error, (2) its accuracy for solving elastic contact problems is nearly independent of the computation domain size accept for the necessary extension, and (3) the DC-FFT algorithm is the fastest among the commonly used contact analysis methods. Figure 4B presents a series of calculation results for the contacts of two honed rough surfaces subjected to several normal loads. The composite root-mean square (RMS) roughness is R q = 0.5 micron. The contact ellipticity is K = 2.0, radii of curvature of the contact bodies are R x = 19.05 mm, R y = 54.165 mm, the equivalent elastic modulus is E ′ = 226.4 GPa, and the maximum Hertzian pressure is 2.72 GPa at the load of 7,680 N. No plasticity is considered in this set of analyses. The complementary conditions for contact modeling are given in the Appendix.

Continuous Convolution and Fourier Transform (CC-FT) for Nominally Flat-Flat Contact Problems
The contact of nominally flat but rough surfaces is a problem with an infinite domain, and it can also be considered as a periodic contact problem, i.e., the contact characteristics in a representative finite region repeat periodically in lateral directions. Since the FRF of Green's function exists and the periodic pressure distribution can be made into Fourier series, the continuous convolution theorem is applicable.
The continuous convolution and Fourier transform (CC-FT) algorithm has been suggested by Liu et al. (2000) and  for this type of problems, which is so named because in the theoretical nature, the continuous convolution theorem is applied. If the FT in Equation (21) and the final IFT are replaced by the discrete Fourier transform and the inverse discrete Fourier transform (IDFT), which is actually that the FT divided by mesh interval is replaced by the DFT and the IFT multiplied by the mesh interval is replaced by the IDFT, the equation below can be executed directly. Here, the FFT and IFFT are applicable to execute the DFT and IDFT efficiently. The CC-FT algorithm is built upon the frequency response functions (FRFs)G and FFT of pressure,p, The FRFs are singular at the coordinate origin, which can be processed with the Gauss quadrature integration method. This CC-FT method is actually what Ju and Farris (1996) and others used before 2000 for non-periodic contact problems where it involved periodic errors. With the CC-FT algorithm, the solution can be computed only in a representative target domain as if this were one period of the rough surface laterally. Therefore, this method is highly efficient, as well as accurate, in computation.

FFT Algorithms for 3D Line-Contact Problems
The 3D line-contact problems involve a limited domain size in one of the lateral directions but a significantly long length in the other. Such a problem can be simplified as a 2D plane-strain issue if the contacting surfaces are ideally smooth and the materials homogeneous. To explore more detail, we are facing 3D problem with mixed issues. Three FFT-based approaches have been developed, to be discussed below, which are the algorithms of mixed discrete-continuous convolution with duplicated padding and FFT (DCD-FFT) (Chen et al., 2008), the hybrid discrete convolution, continuous convolution and FFT (DC-CC-FFT) (Liu and Hua, 2009), and the discrete-periodic convolution with IC summation and FFT (DCS-FFT) (Liu and Hua, 2009;Sun et al., 2020) to consider the effect of roughness and material inhomogeneity. They are mathematically and numerically the same in the finite direction (x), but different in the other dimension (y).

DCD-FFT Algorithm
The DCD-FFT algorithm ( Figure 5A) modifies the DC-FFT algorithm with a treatment in the y direction, simply by duplicating the pressures in the target domain to the extended region in the y direction, called duplicated padding (Chen et al., 2008). The length of the extended region should be at least the same as that of the target domain. This method is not as accurate as the DC-FFT one for non-periodic problems, especially when the extended domain in y is not sufficiently long, because it ignores the deformation influences from the region not included in the calculation, and the IC truncation error plays a role. Therefore, this method is mentioned here for a reference use only. However, the DCD-FFT method can effectively solve the normal deformation of finite-cylinder surfaces (a quarterspace problem), if the cylinders are not too short, with pressure duplication in the extended domain .

DC-CC-FFT Algorithm
The infinite extension of the problem in the y direction qualifies the direct use of the CC-FT method. Thus, the DC-CC-FFT algorithm combines the features of the discrete convolution theorem in the x direction, and the continuous convolution theorem in the y direction, shown in Figure 5B, involving hybrid FRFs and ICs (named ICs-FRF). Because there are two ways to obtain the FRFs for the y-direction solution, two variations of the DC-CC-FFT algorithm can be constructed.

DC-CC-FFT with IC-conversion.
A simple way, with known ICs, described by Sun et al. (2020), to build the DC-CC-FFT algorithm is to process the 2D Fourier transform of the ICs sequentially in the y and the x directions to obtain ICs-FRF, where the FRF is from the IC conversion. After the y-direction FFT, the FRFs are calculated from Equations (14) and (17). Then the wraparound order in the x direction is conducted, followed by the 1D Fourier transform in x. The calculation can be performed in the following steps with a slightly different procedure for the wrap around order.  (14) and (17) in the y direction; 7) Treat D 2N x ×N y with the wrap-around order only in the x direction and get <D > 2N x ×N y which is only term flip because the domain has been enlarged in 4). Here, <D > 2N x ×N y is the ICs after the wrap around order in the x direction but FRF in the y direction, and the latter is converted from the ICs done in the previous step; 8) Apply 1D FFT to <D > 2N x ×N y in the x direction, and the resultant series is denoted as <D > 2N x × N y ; 9) Obtain a temporary frequency series by element-by-element production between < P > 2N x ×N y and <D > 2N x ×N y , which can be expressed as <D > • < P > , where "•" means complex multiplication in an element-byelement manner; 10) Apply a 2D IFFT to <D > • < P > ; 11) Obtain the surface displacement by keeping values within the original target domain, and the resultant series are shown as IFFT <D > • < P > .
The advantage of the algorithm of DC-CC-FFT with ICconversion is that it does not require known FRFs if ICs can be calculated more easily. It is more accurate than the DCD-FFT method but still involves some IC truncation and conversion errors. The IC truncation error can be greatly reduced by using the IC summation method, to be discussed in section DCS-FFT Algorithm. Precise DC-CC-FFT algorithm, presented by Liu and Hua (2009), is to conduct the analytical Fourier transform of Green's function in the length (y) direction, and calculate the ICs with respect to x, to obtain the hybrid ICs-FRFs, and then conduct the wrap-around order in x and the onedimensional FFT of the new ICs-FRFs. This approach is theoretically accurate.
Both DC-CC-FFT algorithms only require extending the domain in the x direction twice the size of the target domain, while the domain size in y is unchanged, which can be as small as possible, as long as it is sufficient represent the necessary features of surfaces and materials.

DCS-FFT Algorithm
Liu and Hua (2009) suggested another approach to deal with 3D line-contact problems, which is to consider the elasticity effect of the entire domain, or a sufficiently large domain, by IC summation, and Sun et al. (2020) have extended this method to solve the contact problems involving inhomogeneous materials. Figure 6A illustrates this idea. We can include M segments of equal length on each side of the target domain in the length direction, while making the target domain size as small as possible, as long as the main contact features are included. N x in the x direction and N y in the y direction may be different.
The ICs in each segment are the same, periodically shifted from those in the target domain. The influences of the pressures in the other segments on the deformation within the target domain can be included via the IC summation in the target domain, given below, with the superscripts for the segment numbers and the subscripts for the y direction coordinates in the Because the ICs in each segment are the same, only one copy of the IC, marked with superscript (0), is necessary, shown below with proper shifts for other copies.  Figure 6B presents the 3D cylindrical contact of a rigid infinitelength cylinder and an elastic half space material (E = 200 GPa, υ = 0.3) with a sinusoidal rough surface, analyzed with the DCS-FFT algorithm. The 3D roughness is periodic in both the x and y directions, and the load is treated periodically in the y direction. With the DCS-FFT method, no edge effect appears at the borders of the target domain, which means that the elasticity effect of neighboring duplicated domains has been properly taken into account. Compared with the DCD-FFT algorithm, the DCS-FFT method does not require a large target domain because of the IC summation. As mentioned before, N y can be very small, e.g., as small as one period of the pressure variation in the y direction, or as short as the length of a representative rough surface area and/or material inhomogeneity region. Among all the three algorithms for 3D line-contact problems, the DCS-FFT and the accurate DC-CC-FFT algorithms are recommended, and the former may be more preferred because it uses the same DC-FFT solution framework, convenient for programming, especially when the DC-FFT software is available as a set of the open-source codes (http://othello.mech.northwestern.edu/ qwang/OpenSourceCodeDCFFT/DC-FFTWeb.htm). It should be mentioned that the DCS-FFT algorithm can also be used to construct a mechanism to solve the nominally flat-flat contact problems, named DCSS-FFT, which further modifies the ICs with 2D IC summations in both x and y directions (Sun, 2020).

DCR-FFT Algorithm
In comparison to the contact issues involving an excitation and response defined in a convolution, the relationship between eigenstrains and the components of the elastic field in a material involves both the convolution and the correlation , and the correlation theorem should be implemented as well. Here, eigenstrain is a generic term for various non-elastic strains, defined by Mura (1993), including thermal strain, plastic strain, fit-induced strain, phase transformation induced strain, and residual strain in a general sense. Many inhomogeneity problems can be solved via the Equivalent-eigenstrain method (EIM) (Eshelby, 1957), and correlations of functions are involved in the mathematical expressions of the eigenstrain-induced field.

Correlation Theorem
The Fourier transform of the correlation of two real datum series is the multiplication of the complex conjugate of Fourier transform of one function and the Fourier transform of the other.
with R Gf (τ ) = R fG (−τ ) , and the Fourier transform is, wheref is the complex conjugate off .

DCR-FFT Algorithm
The DCR-FFT algorithm is an analogy to the DC-FFT algorithm, to be proven below, and it can be combined with the DC-FFT algorithm for a hybrid convolution-correlation operation to solve the inhomogeneity problem illustrated in Figure 1D. Details are given by Liu and Wang (2005), and applications can be found in the works by Liu et al.
Similar to the expression of convolution related to potential function R I = (x − ξ) 2 + y − η 2 for the solution to contact elasticity, Equation (41) involves Here, R has a plus sign inside the third term, which means a 1D cross correlation with [e] on z for the solution.
The discrete form of the 1D cross correlation of influence coefficient D and [e] in Equation (41) is given below where N e is the total number of nodes, and IC components D k+i involves R. More generally, for two series of complex numbers, f and g, where f (t) is the complex conjugate of f (t). Equation (42) is equivalent to the following correlation, This means that the convolution theorems should also be true to Equation (44) or (44'), no matter the series are real or complex. The numerical calculation of correlation can be done either by direct use of the correlation theorem (Equations 39-40) or the convolution given in Equation (44) (or Equation 44') with proper treatment of the two series, and the DC-FFT algorithm directly works for the latter. Caution should be given in conducting the Fourier transform. By analogy, the DC-CC-FFT and CC-FT algorithms can all be extended to include the correlation operation for the analyses of the field due to eigenstrains. Figure 7 shows a case of a cylinder in contact with the rough surface of an inhomogeneous half-space material , solved with the DCS-FFT algorithm. The representative piece of a ground surface is given in Figure 7A) in a mesh of 128 × 128, and R q = 1.18µm. A virtual ground rough surface subjected to the cylindrical contact is formed through periodically extending this patch along the y direction. The target domain has the dimensions of l x ×l y ×l z , and the length are all set to be l x = l y = 8a 3 , l z = 4a 3 . All the cuboidal inhomogeneities are identical, with a x = a y = a z = 0.3a, and they are distributed in the subsurface in the x = 0 plane in the depth of z d = 0.45a. Figure 7B shows the 3D pressure distribution mapped on the contact surface, where sporadic pressure peaks can be found. The contours of the pressure and the von Mises stress distributions in the XOZ and YOZYOZ cross sections are given in Figures 7C,D. The 3D features of the pressure and stress are captured while their length-direction periodicities are also retained. Figure 1E shows the contact involving a multi-layered material. The core analytical solutions to this type of problems are obtained, ignoring the body forces, by solving the governing differential equations in the frequency domain through the Fourier transform. The Fourier-domain solutions become the frequency response functions (FRFs) if surface tractions are unit valued. The conversion through Equations (13)-(18) leads to influence coefficients for the DC-FFT algorithm Yu et al., 2014Yu et al., , 2016. Note that here, the direct inverse Fourier transform of the frequency-domain solution may result in inaccuracy when solving a non-periodic non-infinite contact problem (Liu et al., 2000).

IC Conversion From FRFs
The IC conversion method can be used to other contact cases as well. Our recent studies on the contacts involving magnetoelectroelastic or viscoelastic materials are all through the path of FRFs-ICs conversion and then the DC-FFT algorithm.
Two such examples are given in Figures 8, 9. Figure 8A shows the contact of two magnetoelectroelastic materials subjected to interfacial contact pressure, surface electric and magnetic charges. For this type of multifield material systems, the Fourier-domain solutions can be obtained by solving the coupled mechanical governing equations and the Maxwell equations in the frequency domain through the Fourier transform. The Fourier-domain solutions become the FRFs if the mechanical tractions, surface electric and magnetic charges are unit valued. The inverse Fourier transform through Equations (13) and (18) leads to influence coefficients to construct the DC-FFT algorithm. Figure 8B shows a case of an ellipsoid in contact with the surface of a magnetoelectroelastic material. The major and minor radii of the ellipsoid are R x = 300 mm and R y = 200 mm, and the normal force is P z = 500 N. More details can be found in the work by Zhang X. et al. (2017Zhang X. et al. ( , 2018. Viscoelastic contact problems have drawn a great deal of attention (Goryacheva and Sadeghi, 1995;Chen et al., 2011;Putignano et al., 2015;Stepanov and Torskaya, 2018). It is also convenient to solve certain viscoelastic contact problems in the frequency domain. Here, two types of frequencies are involved, one with respect to time and the other to space. Figure 9A shows the sliding contact of a layer-substrate system, which actually represents four cases for each contacting bodies, (1) a viscoelastic layer on an elastic substrate, as shown, (2) a viscoelastic half space with the layer thickness = ∞, (3) an elastic half space with the layer thickness = 0, and 4) an elastic layer on a viscoelastic substrate by exchanging material properties under certain conditions, all solvable with the same model (Zhang et al., Figure 9B shows a case of a sphere sliding on the surface of a viscoelastic layer-elastic substrate system at a contact speed V = 50 mm/s. The layer thickness is 1 mm, the sphere radius is 10 mm, and the normal load is 1.48 N. Note that in Figure 9A, the interface between the layer and the substrate can be imperfect, and the spatial domain ICs have been solved by Wang et al. (2017a,b) and . More studies related to viscoelastic contact solutions of layered materials with perfect or imperfect interfaces subjected to steady-state and transient conditions can be found in the work by Zhang et al. (2020a,b).

FFT With Non-uniform Mesh
In engineering systems, such as gears, journal bearings, and manufacturing tools, the contact regions can be large but the computation scale in a numerical simulation is limited. It is always a challenge to balance efficiency and accuracy. Take the line contact in a gear or a roller bearing as an example; the middle portion of the contact zone can be analyzed accurately with a coarse mesh, but the edges have to be modeled with a fine mesh in order to describe the drastic pressure variations there. The FEM deals with this type of issues with nonuniform meshes; however, the FFT-based methods built so far, although efficient in a single mesh, are largely confined by the requirement of uniform grids. Sun (2020) has developed a method to extend the FFT-based algorithms to meshes of different densities, in which a finer mesh system is used in specific regions involving high pressure peaks while a coarser one is set in other regions under relatively smooth pressures. Figure 10A illustrates this non-uniform-mesh idea for a roller contact problem, where one side of the regions under the edge effects is meshed denser, while the other side is meshed with a coarse grid. In this particular example, the density of the fine grid is three times that of the coarse grid. The solution on the coarse-grid domain may be run for the entire region, depending on the size of the problem, but that on the finemesh domain is pursued just in a region somewhat larger than what it is designated. Extra data are discarded, and the joint deformations from the two meshes in both designated regions are used to evaluate the gap. This process involves overlapped calculations; however, for problems like the roller contact shown here, the larger the physical domain is, the more the saving of the computation time. Figure 10B shows the comparison of the pressure calculated with a uniform mesh system and non-uniform meshes. The density of the uniform mesh is the same as that of the fine mesh of the non-uniform meshes. The pressure obtained on the effective zone of the fine mesh from the non-uniform mesh solution well matches that on the uniform fine mesh system. When the contact region is large and the pressure distribution is strongly non-uniform, the FFT method with non-uniform meshes offers a more efficient and flexible way to detail the  (Liu et al., 2000), Reprinted by permission from Elsevier, Wear, γ means the discretization density in the frequency domain and χ the ratio for spatial calculation domain extension, (B) CC-FT and DCD-FFT algorithms in solving the contact of a 3D sinusoidal wavy surface and a flat , Reprinted by permission from Elsevier, Tribology International; this is an infinite-domain, or periodic, problem. p means the average pressure.
regions of special concerns, such as edges, interfaces, dents, and defects. Figure 11A is for the results of a finite-domain problem with pressure acting on a region of 2a in length (Liu et al., 2000); it compares the errors from the DC-FFT algorithm structured in different routes with respect to the solution from the direct summation (DS) method. The continuous convolution and Fourier transform (CC-FT) algorithm is used to solve the same finite-domain problem, and its result is also compared. In this figure, the solution methods are named routes, which are Route 2: DC-FFT with ICs from Green's function, Route 3: DC-FFT with ICs from conversion of frequency response functions, Route 4: CC-FT with frequency response functions, Route 5: the classic DS method with ICs from Green's function. Route 1 is not included here, which uses the FEM to calculate the ICs. The DC-FFT solution method (Route 2) appears to be accurate; its results overlap with those by the DS method. This means that the solution is not affected by the calculation domain size. Route 3, however, is different, depending on the discretization density in the frequency domain. The discrete influence coefficient from the frequency response function is the key for error control, the frequency domain sampling intervals, m and n , should satisfy N m ≥2N x and N n ≥2N y , and Equations (16), or (17), or (18) should be implemented, followed by a proper wrap-around order .  Figure 11B shows the behavior of the CC-FT algorithm in solving the periodic problem of a 3D sinusoidal wavy surface in contact with a flat; it compares the numerical solution with the corresponding analytical results results . In addition, the DCD-FFT algorithm is also evaluated. It is evident that the CC-FT algorithm yields the solutions of high accuracy that is within machine precision. The CC-FT algorithm with either the zero-order or the first-order shape function yields nearly the same high accuracy. This indicates that a higher order shape function may be unnecessary if the algorithm is formulated properly.

DISCUSSION AND CONCLUSIONS Accuracy Comparisons
The CC-FT algorithm is superb in analyzing the contact of nominally flat but actually rough surfaces. Figure 12 compares the 3D cylindrical-contact algorithms of DCS-FFT, DCD-FFT, and DC-CC-FFT with IC-sum-conversion in calculating the maximum shear stress τ 1 , given below: The identical target domain (2a × 2a × 2a) is used in all the three FFT-based methods. A large gap is shown between the result curves for the analytical and the DCD-FFT results, caused by the IC truncation error. The utilization of the DC-CC-FFT algorithm (IC-summation-conversion) greatly reduces the truncation error; however, some difference is still visible. The advantage of this DC-CC-FFT algorithm with the summarized ICs is the reduced computational burden because no duplicated-padding and wraparound order are needed in the length direction. The DCS-FFT algorithm shows its accuracy and efficiency in dealing with this type of contact problems without knowing analytical FRFs. Theoretically the precise DC-CC-FT algorithm is more accurate , which is similar to the CC-FT shown in Figure 11B.

Range of Applications of the FFT Approaches
The FFT-based approaches are efficient. Many publications have shown the applications of these FFT algorithms on the contact analyses of elastic fields, plastic transition and yield, flash temperature, thermal stress, partial slip, and contact electrical and magnetic fields, dealing with science and mechanics issues for various systems from traditional mechanical components to emerging sensors and batteries Zhang et al., 2020c). The FFT algorithms discussed above are advantageous in solving the problems that are mathematically described in convolutions, correlations, combined convolutions and correlations, and in frequency response functions as well.
The last one makes the FFT-based approaches more widely applicable and powerful because many governing differential equations can be more easily solved in the frequency domain through the Fourier transform, such as those for functionally graded materials (Ke and Wang, 2006), thermoelastically graded materials (Choi and Paulino, 2008), thermally graded materials (Zhang H. B. et al., 2018), materials with coupled stresses , and magnetoelectroelastic materials (Zhou and Lee, 2013). In addition, because the surface deformation analysis is an intermittent process of the elastohydrodynamic lubrication (EHL) calculation and related modeling, the FFT-based solutions are also building blocks in the models of mixed EHL in general (Liu et al., 2006, 3D line-contact EHL (Ren et al., 2009), finite-roller EHL (Zhang H. B. et al., 2017 with IC-overlapping DC-FFT, Liu et al., 2020 with DCD-FFT), coating EHL (Liu Y. et al., 2007 for single coatings, Wang et al. (2015) for multilayered coatings), plasto-EHL (Ren et al., 2010), EHL of inhomogeneous materials , wear in EHL , EHL of transversely isotropic materials , EHL of artificial joints (Wang and Jin, 2004), and EHL of 2D bearings (Liu and Chen, 2012). Moreover, adhesive contact problems can be solved with the FFT-based methods as well (Pohrt and Popov, 2015;Popov et al., 2017;Rey et al., 2017). It is important to note that different contact types require different convolution theorems and thus different FFT algorithms. As summarized by Wang and Zhu (2019), contact with nominally flat surfaces should be tackled with the continuous convolution and Fourier transform (CC-FT) algorithm. Point-contact problems, either circular or elliptical, are non-periodic, and they should be formulated with the cyclic convolution and solved with the discrete convolution and FFT (DC-FFT) method utilizing zero padding and wraparound order. Three-dimensional cylindrical (line) contact problems involve infinite domain extension in one direction but a finite domain width in the other, for which the discretecontinuous convolutions and FFT (DC-CC-FFT) and the discrete convolution with IC summation and FFT (DCS-FFT) methods should be used. The DCS-FFT algorithm is recommended if ICs are already known. Both the CC-FT and DCS-FFT algorithms are capable of solving periodic problems, and if the latter is used, the IC summation in the other direction is also needed, which makes it the DCSS-FFT algorithm mentioned in section DCS-FFT Algorithm. Figure 13 summarizes the FFT algorithm tree and the application field of each. The roots of this tree are the essential analytical solutions to mechanics and physical problems in the forms of Green's functions, ICs, or FRFs.
Recently, the fast Fourier transform method has been incorporated with the boundary-element method for solving problems with an arbitrary columnar geometry, not confined by the half-space assumption (Benad, 2018(Benad, , 2019. This further extends the breadth of FFT applications. For the plane-strain problems associated with a columnar geometry formed by extruding a 2D shape in the length direction, e.g., the one solved by Benad (2018Benad ( , 2019, the FFT solutions have a lower computational complexity, O(n 3 log n 1.5 ), than the inversion of a standard BEM matrix, O(n 4 ), where n x n is the total number of surface nodes. These problems are, mathematically, in the same nature as that in Figure 2F, automatically meeting the DC-FFT requirement with no need of the domain extension, as indicated in the first row of Table 1 and Figure 13. Likewise, the FFT method should also be directly applicable to plane-stress problems of a disk-like geometry of any shape.

Limitations and Future Developments
The FFT-based methods mentioned here well fit the solutions of many engineering problems, not limited to mechanical contacts, as long as their model formulations contain convolution and/or correlation, or are solvable in the frequency domain, subjected to the assumptions of small deformation (linear or piecewise linear). So far, for counterformal contacts, the characteristic body dimension, such as radius, should be much larger than that of the contact area; for conformal contacts, bodies involved should have axisymmetry or columnar geometries. Generally, these FFTbased solution approaches are confined by uniform meshes. Although section FFT With Non-uniform Mesh has briefly discussed the use of non-uniform meshes, more work is needed to make the non-uniform-mesh FFT algorithms more flexible and more efficient for contact analyses. Large deformation and the effect of body forces are also among the challenges to further developments of FFT-based methods for contact mechanics.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, open-source codes are available online in http://othello.mech.northwestern.edu/qwang/Open SourceCodeDCFFT/DC-FFTWeb.htm; http://othello.mech. northwestern.edu/qwang/OpenSourceCodeEigenstrainFiledHalf Space/2012Web.htm. Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
QW leads the overall writing, LS is on some of the methods and results in Introduction, FFT Algorithms for 3D Line-Contact Problems, DCR-FFT Algorithm, FFT With Nonuniform Mesh, XZ on IC Conversion From FRFs, SL on technical and algorithm review, and DZ on a part of Frequency Response Functions and ICs and rough surface contact results. All authors contributed to the article and approved the submitted version.
If G 3 (x 1 , x 2 ) < 0, make G 3 (x 1 , x 2 ) = 0 and add point (x 1 , x 2 ) to A c (A2) The overall the load balance is whereP M is the normal load. The contact area and pressure distributions can be solved by using the conjugate gradient method (CGM) (Polonsky and Keer, 1999).
In addition, the contact should also satisfy the interfacial conditions for fluxes.