Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 16, N o 2, 2018, pp. 127 - 138 https://doi.org/10.22190/FUME180526023B © 2018 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper FAST NUMERICAL IMPLEMENTATION OF THE MDR TRANSFORMATIONS UDC 539.3 Justus Benad Berlin University of Technology, Berlin, Germany Abstract. In the present paper a numerical implementation technique for the transformations of the Method of Dimensionality Reduction (MDR) is described. The MDR has become, in the past few years, a standard tool in contact mechanics for solving axially-symmetric contacts. The numerical implementation of the integral transformations of the MDR can be performed in several different ways. In this study, the focus is on a simple and robust algorithm on the uniform grid using integration by parts, a central difference scheme to obtain the derivatives, and a trapezoidal rule to perform the summation. The results are compared to the analytical solutions for the contact of a cone and the Hertzian contact. For the tested examples, the proposed method gives more accurate results with the same number of discretization points than other tested numerical techniques. The implementation method is further tested in a wear simulation of a heterogeneous cylinder composed of rings of different material having the same elastic properties but different wear coefficients. These discontinuous transitions in the material properties are handled well with the proposed method. Key Words: Normal Contact, Method of Dimensionality Reduction, Stress, Wear 1. INTRODUCTION The Method of Dimensionality Reduction (MDR) is a simple and convenient tool for the calculation of contact forces between elastic and viscoelastic bodies. It is particularly easy to use for the simulation of axially-symmetric contacts. Since it was first proposed in 2007 [1] the MDR has been applied to a wide range of problems [2]. The method maps a given three-dimensional contact problem to an equivalent contact problem of a transformed indentation profile with a one-dimensional elastic or viscoelastic foundation of independent elements. From a numerical perspective, the solution of the contact problem in the transformed MDR domain is then trivial due to the decoupled degrees of freedom. A variety Received May 26, 2018 / Accepted June 30, 2018 Corresponding author: Justus Benad TU Berlin, Institut für Mechanik, FG Systemdynamik und Reibungsphysik, Str. d.17. Juni 135, 10623 Berlin, DE E-mail: mail@jbenad.com 128 J. BENAD of problems can be solved directly in this domain after the initial transformation to the equivalent problem was performed (see for example [3]). However, there are other problems which require multiple transformations to the MDR domain and back, for example due to a continuously changing indentation profile as it appears in wear simulations, see [4-6]. With such kinds of problems, the main difficulty in achieving an accurate and efficient numerical simulation is the implementation of the MDR transformations. These are given by Abel-like integral equations and it is well known that their numerical treatment is challenging [7, 8]. This work is dedicated to providing a simple and fast numerical method for the implementation of the MDR transformations for axially-symmetric contact problems. The transformations have an integrable singularity which is handled well with the proposed method. The parts of this work are organized as follows: In Section 2, the MDR transformations are rewritten using integration by parts to a form which is well suited for numerical implementation. In Section 3, this numerical implementation technique is explained in detail. Section 4 gives some advice on optimizing the implementation for maximum speed. Section 5 shows exemplary results of the newly proposed technique and highlights its advantages and weaknesses. In Section 6, a small addition to the method is presented to further improve it. In Section 7, the accuracy of the introduced numerical method is compared to other known implementation techniques which rely on the original form of the transformations. In Section 8, an exemplary wear simulation is conducted with the newly introduced technique and compared to the results of other numerical implementation methods. A conclusion is presented in Section 9. This work can be regarded as a small addition to the paper “Method of Dimensionality reduction in contact mechanics and friction: A users handbook” [9]. In the following, only homogeneous elastic material is considered. However, the MDR is applicable also to gradient media [10] and to viscoelastic media [11], which can be treated in a similar manner as described in the present paper. 2. FORMULATION OF THE MDR TRANSFORMATIONS FOR SIMPLE NUMERICAL IMPLEMENTATION The general MDR procedure is fully described in [9]. Three main transformations occur in the method: The transformation of three-dimensional profile f(r) to a one-dimensional profile is 2 2 0 ( ) ( ) x f r g x x dr x r     , (1) the transformation of one-dimensional foundation displacement w1D(x) to three-dimensional normal surface displacement w(r) is 1D 2 2 0 ( )2 ( ) r w x w r dx r x    , (2) and the transformation of one-dimensional force density q(x) to three-dimensional pressure distribution p(r) is Fast Numerical Implementation of the MDR Transformations 129 2 2 1 ( ) ( ) r q x p r dx x r       . (3) The singularity arising at x = r in the numerical summation can be avoided when rewriting Eqs. (1), (2) and (3) to 2 2 0 ( ) ( ) atan ( ) 2 x r g x x f x x f r dr x r             , (4) 1D 1D 2 2 0 2 ( ) ( ) atan ( ) r x w r w r w x dx r x           , (5) and 2 21 1 ( ) log( ) ( ) log( ) ( ) r p r r q r x r x q x dx         (6) using integration by parts. The following example shall illustrate a possible numerical implementation of the three transformations (4), (5) and (6). 3. EXEMPLARY NUMERICAL PROCEDURE Consider a uniform discretization of r∈ [0,L] and x∈ [0,L] with N points each and the same step size 1 L h N   , (7) so that ( 1), ( 1), , {1, 2,..., } n k r h n x h k n k N     . (8) The first and second derivatives of a discretized indentation profile fn = f(rn) can be obtained via central differences: 1 1 2 n n n f f f h      , (9) 1 1 2 2 n n n n f f f f h       . (10) Some care must be taken at the borders. To obtain f1 ′ and f1 ′′ recall that in the present framework of the MDR profile f is axially-symmetric and f(0) = 0. Thus, it is 1 0f  , (11) and 2 1 2 2 f f h   . (12) 130 J. BENAD At the other border the values for fN ′ and fN ′′ can remain undetermined. One-dimensional profile gk can now be obtained with Eq. (4). It is , for 2, 3,..., 1 2 k k k k k g x f x t k N      , (13) where tk is the result of the integral in (4). Using the trapezoidal rule, it is 1 2 2 2 , for 2 4 atan , for 3, 4,..., 1 4 k k k n n k n k n f h k t r f h f h k N x r                            (14) Again, some care must be taken at the borders. To obtain g1, recall that in the framework of the MDR it is g(0) = 0. Thus, it is 1 0g  . (15) At the other border the value for gN can remain undetermined. In a quite similar fashion, normal surface displacement wn can be obtained: The first derivative in Eq. (5) can be obtained as in Eq. (9) using central differences, and the integral can be calculated as in Eq. (14) using the trapezoidal rule. Subsequent smoothing of wn with wn := (wn–1 + wn + wn+1) / 3 increases the accuracy of wn. The third transformation to obtain pn is again similar to the first and second transformation. The derivatives in Eq. (6) can once more be obtained as in Eqs. (9) and (10) using central differences. Then it is 1 1 log( ) , for 2, 3,..., 2 n n n n p r q s n N       (16) where sn is the result of the integral in (6). Using the trapezoidal rule, it is 2 2 2 1 log( ) log( ) , for 2, 3,..., 3 2 log( ) , for 2 2 N n n k n k k k n n n n x q h x r x q h n N s x q h n N                  (17) Note that at kinks of q the term qk ′′ h = (qk+1 – 2qk + qk–1) / h converges to finite values for decreasing step-sizes h. Note also that the summation in Eq. (17) stops at N – 2 because qN ′′ and qN–1 ′′ are undetermined. This is not problematic because in the framework of the MDR it is q ′′ = 0 for sufficiently large x in any way (x > a, where a is the contact radius). Once again, some care must be taken at the borders. One way to approximate p1 is via Taylor series. A first order expansion yields simply 3 2 1 2 2 p p p p    . (18) Fast Numerical Implementation of the MDR Transformations 131 At the other border the values for pN and pN–1 remain undetermined. Again, this is not problematic as long as it is ensured that these last points lie outside the contact area. Then they can simply be set to 1 0 N N p p    . (19) 4. PERFORMANCE Often, the MDR transformations need to be performed repeatedly. One example is that of wear simulations where the transformations (1) and (3) need to be performed many times after each other for a changing indentation profile. In such cases, consider optimizing the implementation of the MDR transformations for maximum speed. For example, when using the transformation technique presented in the example above, note that the summation in Eqs. (13) and (16) can be regarded as a matrix vector product in which the matrix is a kernel which is independent of the indentation profile and can be predefined. This enables full vectorization of the transformations when they are used repeatedly for changing indentation profiles. Also consider the possibility of calculating the derivatives in the transformations such as (9) and (10) via matrix vector multiplication using predefined sparse matrices. 5. EXEMPLARY RESULTS Fig. 1 shows the results of the previously described implementation technique for a conic and parabolic indenter at an exemplary indentation depth d. It becomes apparent that already for as few as N = 51 discretization points a fairly good approximation of the analytical solutions can be achieved. The maximum error of gk, wn, and pn with respect to the analytical solutions for g, w, and p decreases when the number of discretization points N is increased as can be seen in Fig. 2. For most N, the maximum error of pressure distribution pn (a thin grey oscillating line in Fig. 2) is given by the error at the very last discretization point lying within the contact area (highlighted with a star in Fig. 1). Index n of this particular point shall be denoted with n = s. At all other points a much better accuracy is achieved: If point s would be disregarded in the assessment of the maximum error, the upper limit of the grey oscillating line would move down from the dotted line to the dashed one. In the exemplary case of N = 51 which is displayed in Fig. 1 this relatively high error of pn at the point n = s does not immediately become apparent to the viewer due to the large slope of p at the end of the contact area which puts the numerical value close to the analytical curve even if there is a relatively high error. 132 J. BENAD a) b) Fig. 1 Results of the MDR transformations carried out with the numerical procedure described above for N = 51 discretization points, exemplary input parameters of L = 1, E * = 1, d = 0.3 and an exemplary conic indenter (left) given with f(r) = r tan(π/8) and an exemplary parabolic indenter (right) given with f(r) = r 2 /2. The pressure which is obtained at last discretization point within the contact area in this example is highlighted with a star a) b) Fig. 2 Maximum error of gk, wn , and pn for a discretization of N = 51, 52, 53 … 5000, shown for the exemplary inputs L = 1, E * = 1, d = 0.3 and the exemplary conic indenter (left) given with f(r) = r tan(π/8) and the exemplary parabolic indenter (right) given with f(r) = r 2 /2. The oscillating thin grey line shows the maximum error of pn. Its upper limit is marked with the dotted black line. Neglecting the error of pn at the point n = s in the assessment of the maximum error would cause a much lower upper limit which is marked with the dashed black line. Note also that the maximum error of gk for the exemplary conic indenter lies at around 10 -15 and is thus outside the chosen region displayed in the figure Fast Numerical Implementation of the MDR Transformations 133 6. ADJUSTMENTS If needed, one possible way of obtaining better results for ps is by inserting one or more discretization points on a finer grid after point s so that it is no longer the last point in the contact area. For example, one can use a technique such as the one illustrated in Fig. 3. Here, a new discretization point is added right at contact radius a, which is approximated from gk with a simple linear interpolation 1 1 ( )s s s s s s r r a d g r g g        , (20) another discretization point is added in between rs and a at rs + (a – rs) / 2, and at both ends two more points are added, one at rs – (a – rs) / 2 and one at a + (a – rs) / 2. The values for the one-dimensional profile are interpolated linearly from gk to these points. The resulting five equally spaced points are marked with crosses in Fig. 3. Fig. 3 Detailed view of the graph in Fig. 1b, here with additional discretization points at the end of the contact area which are marked with black crosses The desired value for the pressure at point s can now be calculated as in Eqs. (16) and (17) using a new refined grid. The three inner points are used for the summation while the two additional outer points are only there to obtain the derivatives with the central difference scheme. Note that the above method of obtaining a more accurate value for ps does not practically increase the computational time. It is a simple addition of three values, and the three linear interpolations which are needed are also given with small algebraic equations. Compared to the time for the main transformation steps, the time for these additional steps is negligible. However, the small correction reduces the maximum error norm (see Fig. 4). At the end of this section it shall be noted that higher order methods for the calculation of the derivatives and for the numerical integration do not necessarily lead to more accurate results. It is observed that the use of more neighboring points than in Eqs. (9) and (10) for calculating the derivatives tends to decrease the accuracy of the transformations for the contact of the cone and the Hertzian contact. It was also observed that using Simpson’s rule to perform the summation in Eqs. (13) and (16) instead of the trapezoidal rule decreases the accuracy of the transformations. 134 J. BENAD a) b) Fig. 4 Maximum error of gk, wn , and pn as displayed in Fig. 2 here however ps is corrected with the technique explained above. The old upper limit of the maximum error of pn from Fig. 2 is shown with the grey dotted line. The new upper limit of the maximum error of pn after the correction of ps is shown with the black dotted line. In the case of the cone (a) this upper limit falls onto the dashed line marking the maximum error of pn when the point s is disregarded 7. COMPARISON WITH OTHER NUMERICAL TECHNIQUES Recall that at the beginning of the numerical scheme presented above the MDR transformations were rewritten using integration by parts. However, the MDR transformations can also be implemented numerically using their original form of Eqs. (1), (2) and (3) without rewriting them to Eqs. (4), (5) and (6). Here, two such methods which will be called “Method I” and “Method II” shall briefly be discussed. Their accuracy will be compared to the partial integration methods introduced above, which are referred to as “Method III” and “Method IV” in the following. Method I – insertion of h at singularity: One technique for the implementation of the transformations using their original form of Eqs. (1), (2) and (3) is to insert a single increment h at the singularity where x = r, as in 1 2 2 1 k n k k k n k n f f g x h hx r                 , (21) where the first derivative is computed as in Eq. (9) using central differences. This method, however, delivers only very poor results. As can be seen in Fig. 5, the technique requires a number of discretization points which is several orders of magnitude higher in order to reach the accuracy which is achieved by the other implementation techniques. This method is not recommended. Fast Numerical Implementation of the MDR Transformations 135 Method II – implementation of the kernel with its antiderivative: A far better technique for the implementation of the transformations using their original form is to implement the kernel of the transformation using its antiderivative. For the transformation to gk, this translates to 1 2 2 2 2 1 1 atan atann nk k k k n k n n n r r g x h x r x r f h                (22) and for the transformation to pn, one can use 2 2 2 2 1 1 1 log( ) log( ) N k n k k n k n k k n x r x x r x p h q h                  . (23) The first derivatives can once more be obtained via central differences. As can be seen in Fig. 5, the Method II provides a much better accuracy than Method I. Method III – partial integration method: This technique was described in great detail in the first sections of this work. Here the singularity at x = r is avoided through partial integration of the transformations. This leads to alternative formulations of the transformations in which the second derivative of the three-dimensional indentation profile and the deformed elastic foundation occur. Thus, singularities now occur at kinks of these profiles; however, they disappear in the numerical integration, similarly to Method II where small increment h cancels out in Eqs. (22) and (23). Recall, however, that the singularity which is overcome in Method II occurs in the kernel. Method III, however, overcomes singularities which may occur through the shape of the indentation profile or the deformed one-dimensional foundation. Also, the singularity in Method II always influences the transformation values at all discretization points whereas in Method III the singularities through kinks may leave transformation values at some discretization points uninfluenced. In Fig. 5 it can be seen that with Method III the number of discretization points can substantially be reduced to achieve the same accuracy as in Method II. However, it stands out that the maximum error in Method III is still fairly close to the maximum error in Method II. This relatively high maximum error of Method III is generally attained at the end of the contact area. Method IV – partial integration method with small adjustment: The previously described relatively high maximum error of Method III is reduced in Method IV. The simple adjustment through the insertion of an additional discretization point at the end of the contact area is described in Section 6 above. In Fig. 5 it can be seen that with the Method IV the number of discretization points can further be reduced to achieve a certain desired accuracy. 136 J. BENAD a) b) Fig. 5 Upper limits of the maximum absolute error of pn (left graph a) and the mean absolute error of pn (right graph b) compared for the different numerical methods: Method I – insertion of h at singularity, Method II – implementation of the kernel with its antiderivative, Method III – partial integration method (technique from this paper), Method IV – partial integration method with small adjustment (refined technique from this paper). As before, the curves are displayed for the exemplary inputs of L = 1, E * = 1, d = 0.3 and here only for the exemplary parabolic indenter given with f(r) = r 2 /2 8. EXEMPLARY WEAR SIMULATION Apart from a high accuracy, Method III and Method IV may also show an advantage when they are used multiple times on a changing indentation profile, such as in wear simulations. As an example, consider a heterogeneous cylinder which is pressed onto an elastic half-space with normal force FN and moves tangentially with velocity v0. The cylinder shall be composed of rings of different material having the same elastic properties but different wear coefficients k1 and k2 (see Fig. 6a). This setup has recently been studied with the MDR by Li et al. [6] using Archard’s law [12] N wear 0 ( ) ( ) F s V k r r   (24) to model the change of the indentation profile due to wear. Therein, kwear(r) and σ0(r) are wear coefficient, that is, hardness, and with k(r) = kwear(r) / σ0(r) the linear wear is 0 ( ) ( ) ( )f r k r p r v t   . (25) In the following, the same procedure is adopted. It shows that the numerical method which is used for the MDR transformations has a significant impact on the quality of the simulation results. The limiting profile and pressure reached after a long enough wear process are both displayed in Fig. 6b. Profile f is normalized with initial indentation depth d0 = FN / (2aE * ) and the pressure distribution is normalized with p0 = FN / (2πa 2 ). As can be seen in Fig. 6b, the use of Method II to perform the transformations leads to an oscillating error in the results for both the profile and the pressure (a thin grey jagged line). This error does not occur when Method III or Method IV are used (a smooth bold line). For an increasing Fast Numerical Implementation of the MDR Transformations 137 number of discretization points or smaller time steps in the wear simulation the oscillating error which occurs with Method II does not vanish although it can be smoothed out in the post-processing. Method III and IV, however, deliver the undistorted results straight away without the requirement for subsequent corrections. Note that these raw results of the exemplary simulation obtained using Method III and Method IV also reproduce results obtained for validation purposes in [6] with the Boundary Element Method (BEM) [13]. a) b) Fig. 6 Left graph a): A heterogeneous cylinder composed of rings of different material having the same elastic properties but different wear coefficients k1 and k2 is pressed onto an elastic half-space with the normal force FN and moves tangentially with velocity v0. Right graph b): Simulation results for the limiting profile and pressure after a long enough running-in process as obtained with Method II (a thin grey jagged line), and the techniques from this paper – Method III and Method IV (a smooth bold line) with N = 201 discretization points and k2/k1 = 10 9. CONCLUSION A simple implementation technique for the MDR transformations is presented in this work. It relies on integration by parts of the transformations, a central difference scheme to obtain the derivatives, and the trapezoidal rule to perform the summation. It is shown that the results of the method for the contact of a cone and the Hertzian contact converge to the corresponding analytical solutions for an increasing number of discretization points. Therein, the highest error occurs at the border of the contact area. A small refinement to the numerical method has been presented to reduce this error. The introduced method and its refinement are then compared to other numerical techniques which rely on the original form of the transformations. For the tested examples, the newly introduced method and its refinement deliver more accurate results at the same number of discretization points (see Fig. 5). Furthermore, it is shown that apart from a higher accuracy when used once, the presented method and its refinement may have another benefit when used multiple times in 138 J. BENAD wear simulations. In an exemplary simulation, the wear of a heterogeneous cylinder composed of rings of different material having the same elastic properties but different wear coefficients is modeled. These discontinuous transitions in the material properties are handled well by the newly introduced methods, whereas the tested numerical techniques which rely on the original form of the transformations deliver results with a high oscillating error (see Fig. 6). Acknowledgments: The author would like to thank V. L. Popov, M. Heß, Q. Li, E. Willert and P. Diercks for many valuable discussions on the topic and critical comments. REFERENCES 1. Popov, V.L., Psakhie, S., 2007, Numerical simulation methods in tribology, Tribology International, 40(6), pp. 916-923. 2. Popov, V.L., Heß, M., 2016, Method of dimensionality reduction in contact mechanics and friction, Springer, Berlin. 3. Popov, M., Benad, J., Popov, V.L., Heß, M., 2015, Acoustic Emission in Rolling Contacts, Method of Dimensionality Reduction in Contact Mechanics and Friction, Springer, Berlin, pp. 207-214. 4. Dimaki, A., Dmitriev, A., Chai, Y., 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. 5. Dimaki, A., Dmitriev, A., Menga, N., Papangelo, A., Ciavarella, M., Popov, V.L., 2016, Fast high-resolution simulation of the gross slip wear of axially symmetric contacts, Tribology Transactions, 59(1), pp. 189-194. 6. Li, Q., Forsbach, F., Schuster, M., Pielsticker, D., Popov, V.L., 2018, Wear Analysis of a Heterogeneous Annular Cylinder, Lubricants, 6(1), 28. 7. Murio, D., Hinestroza, D., Mejía, C., 1992, New stable numerical inversion of Abel's integral equation, Computers & Mathematics with Applications, 23(11), pp. 3-11. 8. Hansen, E., Law, P., 1985, Recursive methods for computing the Abel transform and its inverse, Journal of the Optical Society of America A, 2(4), pp. 510-520. 9. Popov, V.L., Heß, M., 2014, Method of dimensionality reduction in contact mechanics and friction: A users handbook. I. Axially-symmetric contacts, Facta Universitatis-Series Mechanical Engineering, 12(1), pp. 1-14. 10. 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. 11. Popov, V.L., Heß, M., Willert, E., 2017, Handbuch der Kontaktmechanik, Springer, Berlin. 12. Archard, J., Hirst, W., 1956, The wear of metals under unlubricated conditions, Proc. R. Soc. Lond. A, 236(1206), pp. 397-410. 13. Pohrt, R., Li, Q., 2014, Complete Boundary Element Formulation for Normal and Tangential Contact Problems, Physical Mesomechanics, 17(4), pp. 334-340.