FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 19, No 4, 2021, pp. 805 - 816 https://doi.org/10.22190/FUME210415057W © 2021 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper FFT-BASED IMPLEMENTATION OF THE MDR TRANSFORMATIONS FOR HOMOGENEOUS AND POWER-LAW GRADED MATERIALS Emanuel Willert Technische Universität Berlin, Institute of Mechanics, Germany Abstract. It is shown how the Abel transform solution to the general axisymmetric normal contact problem for homogeneous and power-law graded elastic materials, which is paramount for the solution of different classes of tribological problems with the help of the method of dimensionality reduction (MDR), can be written in terms of explicit convolutions. These can be very efficiently evaluated with the 1D Fast Fourier Transform (FFT), which reduces the numerical complexity of the transformations from the order of N2 for the standard matrix-vector-multiplication (MVM) algorithm to the order of N. Convergence and performance of the proposed method are studied in detail. As an illustrative example a fretting wear simulation based on the new implementation is shown, the results of which are compared to the standard MVM implementation. Key words: Method of dimensionality reduction, Fast Fourier transform, Boussinesq problem in contact mechanics, Power-law graded materials 1. INTRODUCTION The Boussinesq problem, i.e., the frictionless normal contact problem of elastic bodies, is one of the basic fundaments of contact mechanics. In the case of axisymmetric, convex contacting bodies within the half-space approximation, several qualitatively different classes of contact problems can be reduced to the frictionless, non-adhesive normal contact [1], e.g. the JKR-adhesive normal contact [2] the tangential contact in Cattaneo-Mindlin approximation ([3, 4]) or – based on the correspondence principle between boundary value problems of linear elasticity and linear viscoelasticity – the normal contact of viscoelastic bodies ([5]). The solution to the general axisymmetric Boussinesq problem for homogeneous materials was first discovered by Schubert [6]. A very clear and convenient interpretation of the solution is given by the method of dimensionality reduction (MDR, [7]); the solution Received April 15, 2021 / Accepted August 02, 2021 Corresponding author: Willert Emanuel Technische Universität Berlin, Institut für Mechanik, Sekr. C8-4, Straße des 17. Juni 135, 10623 Berlin, Germany E-mail: e.willert@tu-berlin.de 806 E. WILLERT consists of a series of Abel-like integral transformations, which will be given in detail below and whose numerical implementation in terms of matrix-vector-multiplications (MVM) is straight-forward and usually does not pose severe problems [8]. However, the computational complexity of the MVM is of the order of N2, where N denotes the number of elements/nodes used for the numerical transformations. If the spatial resolution needs to be very high, or the transformations need to be executed very often, e.g., for the simulation of fretting wear [9] or other highly transient phenomena, this complexity might still lead to undesirably high computational costs. On the other hand, as will be shown below, the transformations can be written in terms of simple convolutions, whose numerical implementation in one dimension can be accelerated with the (discrete) Fast Fourier Transform (FFT). The 1D-FFT reduces the numerical complexity to the order of N (omitting logarithmic factors). Therefore, in the present manuscript it will be shown, how the general solution to the axisymmetric Boussinesq problem (in its MDR interpretation) can be obtained extremely fast based on the 1D-FFT. For homogeneous materials this is done in Section 2. Nevertheless, the MDR has also been extended to explicitly incorporate contacts of functionally graded materials (FGM). Within an FGM, the mechanical properties vary continuously over the volume, the material is inhomogeneous. Examples are hardened surfaces as well as diverse biological and biotechnological systems such as bones, joints or their artificial variants, adhesive devices (e.g., on the feet of geckos) and cell membranes [10]. One of the most important classes of FGM are materials in which the elastic modulus E varies with the depth z. Over time, different concrete functions E = E(z) have been considered (see e.g. the treatise by Selvadurai [11]); the only case that allows a largely analytical treatment seems to be the power law E ~ zk. The solution to the axisymmetric Boussinesq problem for power-law graded materials stems from Giannakopoulos & Suresh [12] and Jin et al. [13]. In the latter work as well as by Chen et al. [14] also the respective JKR-adhesive problem was solved. The implementation of these solutions into the framework of the MDR has been demonstrated by Heß [15]. Moreover, tangential contact problems of power-law graded materials can, within the Cattaneo-Mindlin approximation, be analysed with the MDR [16]. As the correspondence principle applies to viscoelastic graded materials, if the spatial and temporal variations of the moduli are separable [17], also contacts of viscoelastic power-law graded materials can be studied within the framework of MDR [18]. This is particularly interesting for biological graded materials, since these usually have visco- or poroelastic properties [10]. The axisymmetric normal contact solution for power-law graded materials is most conveniently formulated in terms of generalized Abel transforms – the homogeneous case is, of course, always included as a limiting case, if the exponent k of the power-law equals zero – which can also be written as explicit convolutions and thus very efficiently evaluated by the 1D-FFT, as will be shown in Section 3. The application of this method for the simulation of fretting wear can be very fruitful in several applications, as it has been demonstrated that graded materials can offer a solution to the wear-fatigue dilemma in fretting [19]. Finally, Section 4 is devoted to discussions and conclusive remarks. FFT-based Implementation of the MDR transformations for homogeneous and power-law graded materials807 2. TRANSFORMATIONS FOR HOMOGENEOUS MATERIALS 2.1. Profile Transformation Consider a rotationally symmetric rigid indenter with the convex profile f(r), that is pressed into an homogeneous, isotropic, linearly elastic half-space with the effective Young’s modulus E* = E/(1–υ2), with Young’s modulus E and the Poisson ratio υ (the determination of f and E* for the contact of two elastic bodies is elementary and can be found in any textbook on contact mechanics [1]). Then, the contact solution is most conveniently formulated in terms of the equivalent profile [1] 2 2 0 d ( ) ( ) , ( ) : d . d x G rf r g x G x r x x r = = −  (1) For example, the function g will give the relations between the macroscopic contact quantities – contact radius a, indentation depth δ and normal force FN [1] * * 0 ( ), 2 2 ( ) d . a N g a F E a E g x x= = −   (2) To write the transform (1) as an explicit convolution, following Bracewell [20] we introduce the substitutions 2 2 ˆ ˆ: , : , ( ) ( ), ( ) ( ).x r f r f G x G= = = =    (3) Note that the function g will be given by the scaled derivative ˆ ˆd d d ˆ( ) 2 ( ). d d d G G g x g x = = =      (4) After substitution and integration by parts one obtains 0 0 ˆ1 ( ) ˆˆ ( ) d ( ) d , 2 f G f = = − −              (5) where the prime denotes the first derivative with respect to the given argument. The last integral shall be evaluated numerically on equidistant one-dimensional arrays, 2 : ( 1) , : ( 1) , : , , 1: , 1 i j L i h j h h i j N N = − = − = = −   (6) where N is the number of elements and L is the length of the simulation area for the spatial coordinates x and r. Note that the discretization is equidistant in the squares of the spatial coordinates. Applying the trapezoidal rule, it is (by convention the sum over an empty set shall be zero instead of ill-defined) 1 3 1 2 1 ˆ ˆˆ 1 . 2 i i j j G h f i f i j − =    = − + −     (7) 808 E. WILLERT To demonstrate the performance advantage of the 1D-FFT, the sum in the last equation has been evaluated both with MVM and as an FFT-accelerated convolution. The calculation time for the pure transformation step (without build-up of a priori known matrices or vectors) as a function of the number of array elements N is shown in Fig. 1 in double- logarithmic scale. While the absolute value for the calculation time is not too relevant, as it depends on the computational setup, it can be easily seen, that the numerical complexity of the MVM is of the order of N2, while the FFT complexity only increases with N. The scatter in the calculation time for the FFT is probably due to performance variance of the FFT depending on N. Note that the numerical error mainly originates from the discretization and the trapezoidal rule. Therefore, when comparing the MVM and FFT solutions with analytical solutions, both variants exhibit the same convergence behaviour. Fig. 1 Comparison of the calculation time for the transformation step (7) with matrix- vector-multiplication (MVM) and as a FFT-accelerated convolution, depending on the number of array elements used for the transformations. 2.2. Transformation of Local Displacements Several classes of tribological problems can already be solved with the MDR only based on the one transformation given in the previous subsection, as the relations between the macroscopic contact quantities (contact radii, forces, and macroscopic displacements) are the same in the MDR model as in the axisymmetric original system [7]. This, e.g., already allows the extensive study of different classes of dynamic contact problems, like impacts, with the MDR [21]. However, local quantities, i.e., local displacements and contact stresses, do not exist explicitly within the MDR model but must obtained from the respective MDR quantities by yet different integral transforms. The local displacements (both normal and tangential) in the axisymmetric original system, u(r), and its MDR image, U(x), are connected via the transform [1] 2 2 0 2 ( ) ( ) d . r U x u r x r x = −   (8) FFT-based Implementation of the MDR transformations for homogeneous and power-law graded materials809 Note that formally this is also the inverse transform to Eq. (1), to obtain the axisymmetric profile f from the equivalent profile g. Introducing ( ) : ( ) (0),U x U x U= − (9) and the substitutions ˆˆ( ) ( ), ( ) ( ),u u r U U x= =  (10) one obtains 0 ˆ2 ( ) ˆ( ) (0) d . 2 U u U= + −          (11) Now, there are (among others) two ways to integrate this expression by parts, 0 2 ˆˆ( ) (0) ( ) arccos d ,u U U   = +              (12) or, 0 0 ˆ2 ( ) d ( ) ˆ( ) (0) lim ( ) d , ( ) : . dx U x U u U x→     = + − + − =                          (13) Obviously, the second variant has several disadvantages: the occurring limit value must exist and β(0) should be finite for the trapezoidal rule to be applicable without problems (that is e.g. not the case if U(x) is parabolic); however, the remaining integral in Eq. (13) is a convolution which can be evaluated once again very efficiently with the 1D-FFT. On the other hand, if the given restrictions are not met, a matrix-vector-multiplication based on the less restrictive Eq. (12) will be preferable. 2.3. Transformation between Line Loads and Stresses The line loads in the MDR model, q(x), and the contact stresses in the axisymmetric original, σ(r), are connected via the Abel inverse transform [1] 2 2 1 ( ) ( ) d . a r q x r x x r  = − −   (14) Note that, once again, this transform is valid for both the normal and the tangential stresses or line loads. Introducing the usual substitutions, ˆ ˆ( ) ( ), ( ) ( ),r q x q= =    (15) integrating by parts, 2 2 2 2ˆ1 ( ) 2ˆ ˆ ˆ( ) d ( ) ( ) d , a a q a q a q    = − = − − + −  −                    (16) and discretizing based on Eq. (6) and the trapezoidal rule, one obtains 810 E. WILLERT 2 1 3 1 2 ˆ ˆ( ) , 1: 1, ( 1) , 1 ˆ ˆ . 2 j j j i i j j h q t j a h t h j q q i j − = +  = − − + = − = −      = − + −              (17) The remaining sum will again be evaluated with the 1D-FFT. In Fig. 2 the results of a convergence study for that algorithm with an analytical solution, 2 2 24ˆ ˆ( ) , ( ) (2 ) 3 q a a= = − − +       (18) are given. The mean relative error (the star denotes the analytical solution) * 1 * 1 ˆ ˆ1 : ˆ1 j j j j      − = −  = −  (19) is shown as a function of the number of elements used for the transformation in double- logarithmic scale. All derivatives were determined by second-order finite differences. The convergence is obviously fast and stable. Fig. 2 Convergence study for the mean relative error depending on the number of array elements for the FFT-based implementation of the discrete transform (17). 3. TRANSFORMATIONS FOR POWER-LAW GRADED MATERIALS In this section the convolution formulation of the relevant transformations is shown for power-law graded materials. Naturally, all previous results for homogeneous media are recovered, if the exponent k of the power-law of the elastic grading equals zero. FFT-based Implementation of the MDR transformations for homogeneous and power-law graded materials811 3.1. Profile Transformation For inhomogeneous materials with an elastic grading in the form of a power-law the profile transformation between the axisymmetric profile f and the profile g reads [1] 1 2 2 1 2 2 1 0 0 ( )d d ( ) , ( ) : ( ) ( ) d . 1 d( ) x xk k k k f r r x G g x x G x f r x r r k xx r − − + −  = = = − +−   (20) Introducing the same substitutions as in Eq. (3) we obtain 1 0 ˆˆ ( ) ( ) ( ) d , k G f += −       (21) and 1 ˆ2 d ˆ( ) ( ). 1 d kG g x g k − = = +    (22) With the same discretization as in Eq. (6) and the trapezoidal rule we finally arrive at 1 3 1 1 1 2 1 ˆ ˆˆ ( 1) ( ) , 2 i k k k i j j G h f i f i j − + + + =    = − + −     (23) the sum contribution of which is a discrete convolution, which once again can be efficiently evaluated with the 1D-FFT. 3.2. Transformation of Local Displacements The transformation rule for the local displacements in the MDR model and the original axisymmetric system in the case of power-law graded materials is given by [1] 2 2 1 0 2 ( ) ( ) cos d . 2 ( ) r k k k x U x u r x r x +   =     −    (24) Note that, once again, this formally is also the inverse transform of Eq. (20) to obtain the axisymmetric profile f from the profile g. With the same procedure as in subsection 2.2 we obtain two different possible formulations after integrating by parts. The generalization of Eq. (12) is given by ( ) 0 1 2 1 2ˆ ˆˆ( ) (0) ( ) cos ; d , 2 1 1 1 3 ( ; ) : , ; ; , 1 2 2 2 k k u U U U H k k k k H y k y F y k +    = + −        + + +  =   +             (25) with the hypergeometric function 2F1. The alternative, more restrictive formulation, which on the other hand includes an explicit convolution, i.e., the generalization of Eq. (13), reads 812 E. WILLERT 1 1 10 0 1 2 ( ) ˆ( ) (0) cos lim ( ) ( ) d , 2 ˆd ( ) ( ) : . d k k kx k k U x u U x U − − −→ −     = + − + −            =                    (26) The discretization procedure for Eq. (26) works completely analogous to the previously showed one and shall therefore be omitted here for brevity. In Fig. 3 the results of a convergence study for that algorithm with an analytical solution, 4 4(1 )(3 ) ( ) , ( ) 8 k k g x x f r r + + = = (27) are shown. The exponent of the power-law was chosen to be k = -0.5 (note that the value of k has practically no influence on the convergence behaviour). The mean relative error (defined as before) is shown as a function of the number of elements used for the transformation in double-logarithmic scale. All derivatives were determined by second-order finite differences. The convergence is stable but slower than the one for the pressure transform in Fig. 2, the gradient (in double-logarithmic terms) of the curve in Fig. 3 is close to minus one. Fig. 3 Convergence study for the mean relative error depending on the number of array elements for the FFT-based implementation of the transform (26) of the local displacements for power-law graded materials. 3.3. Transformation to Obtain Stresses For power-law graded materials it is unfortunately not that easy to directly transform between the line loads in the MDR model and the contact stresses in the axisymmetric system. However, in the elastic case, transformations exist between the local displacements in the MDR system and the contact stresses. For example, the pressure distribution is [1] 2 2 1 ( )d ( ) , ( ) a N k r c g x x p r x r −  = −   (28) FFT-based Implementation of the MDR transformations for homogeneous and power-law graded materials813 with a lengthy constant cN, which can be found in the respective literature [1] and which reduces to E* in the homogeneous limit, and a similar transform exists for the determination of the frictional shear stresses in the elastic tangential contact [1]. After completely the same procedure as in subsection 2.3, the discretized pressure distribution, i.e., the generalization of Eqs. (17), is given by 1 1 2 1 3 1 1 1 2 ˆ ˆ( ) , 1: 1, ( 1) , (1 ) 1 ˆ ˆ( ) ( ) . 2 k kN j j k k k j i i j c p j h g t j a h k t h j g g i j + + − + + + = +  = − − = − = −   +    = − + −             (29) 4. EXEMPLARY FRETTING WEAR SIMULATION Let us demonstrate the application of the proposed method to a simple fretting wear simulation. Fretting is a phenomenon that occurs in dry, partial slip, oscillating contacts; it causes both wear and fatigue and, despite its high localization, is highly relevant for the lifetime of various tribological systems in technology and biology [19]. Consider the following contact situation: an initially parabolic indenter with the radius of curvature R is pressed into a half-space (which of the bodies is elastically deformable, does not actually matter) by a constant indentation depth d. Moreover, the contact experiences small tangential oscillations with a displacement amplitude uA. These oscillations shall be too small to cause gross sliding of the contacting bodies; however, due to the partial slip in the contact there will be wear, which can be modelled by the Archard law; as the wear will happen on a much larger timescale than the period of one fretting oscillation, we can treat one fretting cycle as the time step for the wear simulation. Hence, the local change of profile during one fretting cycle is given by [9] wear 0 ( ; ) ( ; ) ( ; ), k f r t p r t u r t =   (30) with the wear coefficient kwear, the hardness of the worn body σ0 and the slip displacement distribution Δu, which can be calculated analytically, as is shown in another manuscript [22] by the author. Therefore, the wear simulation only requires the two-step transformation f → g – see Eq. (1) – and g → p – see Eq. (28) with k = 0 – in every fretting cycle. In Fig. 4 the results for the time-evolution of the worn profile for different numbers of normalized fretting cycles are shown in normalized variables. The characteristic number of fretting cycles is given by [9] 2 0 0 0 wear 0 , A a d N k F u   = (31) with the initial normal force F0. The number of elements was chosen to be N = 2001. The results on the left have been obtained with the standard matrix-vector multiplication procedure; the results on the right were determined with the FFT-based implementation proposed in the present manuscript. There are no notable differences between both variants. However, the FFT-based implementation operates faster by a factor of five. 814 E. WILLERT Moreover, this speed advantage of the new implementation obviously increases, if more elements are used for the transformation to obtain a finer resolution. (a) (b) Fig. 4 Time-evolution of the worn profile f, normalized for the fixed indentation depth d, as a function of the radial coordinate r, normalized for the initial contact radius a0, for tangential fretting wear of an initially parabolic indenter and different numbers of normalized cycles, N/N0. The dashed black line denotes the limiting no-wear profile [23]. (a) Results for the semi-analytic method with matrix-vector-multiplication (b) Results for the semi-analytic method accelerated by the FFT 5. CONCLUSIONS It has been shown, how the Abel-like integral transform solution to the general axisymmetric Boussinesq problem can be written in terms of real convolutions to enable the evaluation of the integral transforms based on the 1D-FFT. Both elastically homogeneous materials and those with an elastic grading in the form of a power-law were considered. While the convolution formulation works perfectly for the profile transformation and the determination of the contact stress distribution, it is significantly more restrictive than the standard formulation in case of the determination of the local displacement distribution. Moreover, it was demonstrated that the FFT-based solution operates faster than the standard formulation, if the number of array elements exceeds 103 (obviously, the advantage of the FFT will grow, if the number of array elements further increases). While modern computers are, of course, able to quickly perform the transformations in the traditional form, there might be applications for which the reduction in calculation time and memory space (for a very fine discretization) can come in handy, if the contact simulation is part of a more global (hybrid) numerical routine, e.g., for the coupled analysis of wear and fatigue in fretting. As the main numerical issues in the transformation result from the finite difference derivatives and the trapezoidal rule, the FFT-based solution exhibits practically the same convergence and stability as the standard procedure. The FFT-formulation can be used for highly transient contact simulations; in this case, one, however, must keep in mind that the one-dimensional arrays used for the FFT are equidistant in the squares of the spatial coordinates; this also leads to a finer resolution at the edge of contact, which might be beneficial in several circumstances. FFT-based Implementation of the MDR transformations for homogeneous and power-law graded materials815 There exists a wide literature on the numerical implementation of the Abel transform and its inversion (see the review paper by Hickstein et al. [24] and the references therein), mainly in the framework of image processing. However, the algorithms used in that context are optimized to handle very noisy data and therefore prioritize stability over efficiency. This is usually not necessary in the framework of axisymmetric contact mechanics. Acknowledgement: This work has been supported by the German Research Foundation under the project number PO 810/66-1. Moreover, the author is very grateful to Mr. Justus Benad for valuable discussions on the topic. REFERENCES 1. Popov, V.L., Heß, M., Willert, E., 2019, Handbook of contact mechanics: exact solutions of axisymmetric contact problems, Springer-Verlag, Berlin Heidelberg, 347 p. 2. Barquins, M., Maugis, D., 1982, Adhesive contact of axisymmetric punches on an elastic half-space: the modified Hertz-Huber’s stress tensor for contacting spheres, Journal de Mécanique Théorique et Appliquée, 1(2), pp. 331-357. 3. Jäger, J., 1998, A New Principle in Contact Mechanics, Journal of Tribology, 120(4), pp. 677-684. 4. Ciavarella, M., 1998, Tangential Loading of General Three-Dimensional Contacts, Journal of Applied Mechanics, 65(4), pp. 998-1003. 5. Lee, E.H., Radok, J.R.M., 1960, The Contact Problem for Viscoelastic Bodies, Journal of Applied Mechanics, 27(3), pp. 438-444. 6. Schubert, G., 1942, Zur Frage der Druckverteilung unter elastisch gelagerten Tragwerken, Ingenieur-Archiv, 13(3), pp. 132-147. 7. Popov, V.L., Heß, M., 2015, Method of Dimensionality Reduction in Contact Mechanics and Friction, Springer- Verlag, Berlin Heidelberg, 265 p. 8. Benad, J., 2018, Fast numerical implementation of the MDR transformations, Facta Universitatis-Series Mechanical Engineering, 16(2), pp. 127-138. 9. Dimaki, A.V., Dmitriev, A.I., Chai, Y.S., Popov, V.L., 2014, Rapid simulation procedure for fretting wear on the basis of the method of dimensionality reduction, International Journal of Solids and Structures, 51(25-26), pp. 4215-4220. 10. Liu, Z., Meyers, M.A., Zhang, Z., Ritchie, R.O., 2017, Functional gradients and heterogeneities in biological materials: Design principles, functions, and bioinspired applications, Progress in Materials Science, 88, pp. 467-498. 11. Selvadurai, A.P., 2007, The Analytical Method in Geomechanics, Applied Mechanics Reviews, 60, 3, pp. 87-106. 12. Giannakopoulos, A.E., Suresh, S., 1997, Indentation of solids with gradients in elastic properties: part II. axisymmetric indentors, International Journal of Solids and Structures, 34(19), pp. 2393-2428. 13. Jin, F., Guo, X., Zhang, W., 2013, A unified treatment of axisymmetric adhesive contact on a power-law graded elastic half-space, Journal of Applied Mechanics, 80(6), 061024. 14. Chen, S., Yan, C., Zhang, P., Gao, H., 2009, Mechanics of adhesive contact on a power-law graded elastic half- space, Journal of the Mechanics and Physics of Solids, 57(9), pp. 1437-1448. 15. Heß, M., 2016, A simple method for solving adhesive and non-adhesive axisymmetric contact problems of elastically graded materials, International Journal of Engineering Science, 104, pp. 20-33. 16. Heß, M., 2016, Normal, tangential and adhesive contacts between power-law graded materials, Presentation at the Workshop on Tribology and Contact Mechanics in Biological and Medical Applications, Technische Universität Berlin, November 14-17, Berlin. 17. Paulino, G.H., Jin, Z.H., 2001, Correspondence principle in viscoelastic functionally graded materials, Journal of Applied Mechanics, 68(1), pp. 129-132. 18. Willert, E., 2020, Ratio of loss and storage moduli determines restitution coefficient in low-velocity viscoelastic impacts, Frontiers in Mechanical Engineering, 6, 3. 19. Willert, E., Dmitriev, A.I., Psakhie, S.G., Popov, V.L., 2019, Effect of elastic grading on fretting wear, Scientific Reports, 9, 7791. 20. Bracewell, R.N., 2000, The Fourier Transform and Its Applications, 3rd Edition, New York: McGraw-Hill Book Company, 636 p. 816 E. WILLERT 21. Willert, E., 2020, Stoßprobleme in Physik, Technik und Medizin: Grundlagen und Anwendungen, Springer Vieweg, Berlin, 241 p. 22. Willert, E., 2021, A simple semi-analytic contact mechanical model for tangential and torsional fretting wear of axisymmetric contacts, Symmetry, submitted manuscript. 23. Popov, V.L., 2014, Analytic solution for the limiting shape of profiles due to fretting wear, Scientific Reports, 4, 3749. 24. Hickstein, D.D., Gibson, S.T., Yurchak, R., Das, D.D., Ryazanov, M., 2019, A direct comparison of high-speed methods for the numerical Abel transform, Review of Scientific Instruments, 90, 065115.