Artificial Neural Network Architecture for Prediction of Contact Mechanical Response

Predicting the contact mechanical response for various types of surfaces is and has long been a subject, where many researchers have made valuable contributions. This is because the surface topography has a tremendous impact on the tribological performance of many applications. The contact mechanics problem can be solved in many ways, with less accurate but fast asperity-based models on one end to highly accurate but not as fast rigorous numerical methods on the other. A mathematical model as fast as an asperity-based, yet as accurate as a rigorous numerical method is, of course, preferred. Artificial neural network (ANN)–based models are fast and can be trained to interpret how in- and output of processes are correlated. Herein, 1,536 surface topographies are generated with different properties, corresponding to three height probability density and two power spectrum functions, for which, the areal roughness parameters are calculated. A numerical contact mechanics approach was employed to obtain the response for each of the 1,536 surface topographies, and this was done using four different values of the hardness per surface and for a range of loads. From the results, 14 in situ areal roughness parameters and six contact mechanical parameters were calculated. The load, the hardness, and the areal roughness parameters for the original surfaces were assembled as input to a training set, and the in situ areal roughness parameters and the contact mechanical parameters were used as output. A suitable architecture for the ANN was developed and the training set was used to optimize its parameters. The prediction accuracy of the ANN was validated on a test set containing specimens not seen during training. The result is a quickly executing ANN, that given a surface topography represented by areal roughness parameters, can predict the contact mechanical response with reasonable accuracy. The most important contact mechanical parameters, that is, the real area of contact, the average interfacial separation, and the contact stiffness can in fact be predicted with high accuracy. As the model is only trained on six different combinations of height probability density and power spectrum functions, one can say that an output should only be trusted if the input surface can be represented with one of these.


INTRODUCTION
Surface topography plays an extremely important role in processes such as wear, friction, lubrication, sealing, contact resistance, and heat conduction. This is due to that the roughness causes local contacts between the surfaces, as in mixed lubrication, thus governing how the surface deforms and behaves in contact and defining the boundary friction and the real area of contact. The surface topography may be characterized by a number of areal roughness parameters defined in ISO 25178, see ISO Central Secretary (2012). These parameters have, however, limited correlation to the real area of contact, as well as to friction and wear processes, especially if only a few out of the complete set of field parameters are considered.
By using computational contact mechanics, we can estimate the real area of contact for a surface with given topography and also show how the areal roughness parameters change inside the contact. Notice that, an accurate and reliable result requires a highly resolved surface topography measurement. Thereby, the mesh considered in the numerical solution procedure will have to be of equal resolution, which, in turn, increases the computational time significantly. From an engineering point of view, a noniterative model that swiftly yields relatively exact predictions of contact mechanics parameters, such as the real area of contact, would therefore be highly desired.
The Greenwood and Williamson (GW) theory (Greenwood and Williamson, 1966;Greenwood and Tripp, 1970) has been and still is very frequently used. Note that it was Archard (1957) who laid the foundation for most of the (multi)asperity-based type of models known of today (Nayak, 1971;Onions and Archard, 1973;Bush et al., 1975;Bush et al., 1979;Carbone, 2009;Greenwood et al., 2011). The Persson contact mechanics theory (Persson, 2006;Yang and Persson, 2008) is also a frequently used tool. Although being highly useful models that provide insight and yield rapid predictions, they are based on assumptions, making them not always very accurate. The GW theory assumes that the asperities at the surfaces exhibit Gaussian probability distributions. The asperities are also assumed to deform independently of each other which leads to that GW theory is applicable only when the contact area is small (compared with the nominal contact area). Persson's theory assumes that the surfaces exhibit Gaussian height probability distributions and it considers interasperity coupling. Although Persson's theory might not be very accurate for small real area of contact, it applies well to study the complete contact (see e.g. ). The study by Müser et al. (2017) summarizes findings obtained with various kinds of models, including asperity-based ones and Persson's theory. Moreover, results from numerical brute force methods, all-atoms-based models, and experiments were presented as well. It was concluded that 1) rigorous numerical brute force approaches yield almost identical results on all properties, 2) Persson's theory, all-atom simulations, and experiments could be used to identify the correct trends, and almost exact numbers for some properties, and 3) asperity models predicted the real area of contact rather well and provided alternative interpretations for other properties. It would be very useful if it was possible to obtain a mathematical model for fast calculation, which is as accurate as the rigorous models are, when predicting contact mechanics parameters such as real area of contact.
The ideal situation would be to describe surface topography by its height probability distribution and its power spectrum, which constitute the complete description. However, this complicates the analysis, and if a subset of the areal roughness parameters ISO Central Secretary (2012) would be sufficient, it would facilitate the analysis tremendously. In this study, we will present an artificial neural network (ANN)-based model. This model acts as a transfer function, taking areal roughness parameters as input and predicts the real area of contact and other contact mechanics parameters. A similar ANN-based approach has been used in contact mechanics before (see (Rapetto et al., 2009)). Other examples where ANN-based approaches have been used in tribology are Nasir et al. (2010), Nirmal (2010),Ćirović et al. (2012, and Moder et al. (2018). If an ANN, which executes much faster than a computational contact mechanics approach, is well designed, trained, and tested, it can thus provide reasonably accurate predictions of tribological performance parameters very rapidly.
The idea with the present work is to generate thousands of surfaces by means of the method developed by Pérez-Ràfols and , and to compute parameters, such as the real area of contact and areal roughness parameters when these surfaces are pressed into contact with a flat rigid counter surface. To this end, we will use the computational contact mechanics approach presented by Almqvist et al.( 2007), which was further developed by Sahlin et al. (2010).
The ANN is trained to find the relationship between the surfaces' original, the in-contact, that is, in situ areal roughness parameters and the contact mechanics parameters, for a range of loads, spanning from no load at all to a load that causes nearly as much as 50% real area of contact.

METHODS
This section presents, in a workflow order, the implementation of the ANN. It starts with describing surface topography generation, followed by preprocessing and a brief description of the contact mechanics approach adopted, and it ends with a presentation of the architecture of artificial neural network that was developed herein.

Surface Topography Generation
Training neural networks requires large data sets. Therefore, it is necessary to generate a wide range of different surface topographies. The surface randomization algorithm developed by Pérez-Ràfols and Almqvist (2019) was employed to randomly generate 2,022 surfaces topographies with given height probability distribution (HPD) and a power spectrum (PS). The HPD and PS can be mathematically modeled by classical distribution and spectrum functions, but they may also be obtained (and adapted) from measured surface topographies. In this work, mathematical models for Gaussian, bi-Gaussian and Weibull functions, and self-affine and exponential PS functions were used. The reader is referred to Pérez-Ràfols and Almqvist (2019) for a precise description of these HPD and PS. A surface topography generation selection scheme is depicted in Figure 1, where one type of HPD and PS is selected together with the corresponding shape-defining parameters, that is, c, k, H, β, α, q 0 , and q 1 . Remark that the HPDs are defined with zero mean value and unit standard deviation. With these constrains, the Gaussian HPD requires no input, while the bi-Gaussian and the Weibull distributions may be defined using one parameter, that is, c and k, respectively. Specifying the PS requires four parameters, that is, the Hurst exponent H for the self-affine and the parameter β, which defines the autocorrelation length 1/β, for the exponential function, plus α for the anisotropy and the wave numbers (q 0 , q 1 ) that specify the frequency bandwidth. A 256 × 256 mesh was considered affordable for the grand total of 7,602 elastoplastic contact mechanics simulations performed. The mesh limits the choice of the high frequency cutoff, and in order to resolve the shortest wavelength with at least eight nodes, it was chosen as q 1 32 (in terms of its wave number). This parameter was also kept constant when generating the sets of surfaces for training, testing, and validation. Thus, each surface is generated based on a pair of HPD and PS functions and five corresponding numerical parameters, except for the Gaussian which needs four. Figure 2 shows an example of a generated surface, using the bi-Gaussian HPD model and the self-affine PS model. The corresponding parameter settings are displayed to the right.
The parameter space for the surface dataset used for training was defined by four equidistantly spaced values for each of the seven parameters (q 1 was kept constant). In this way, a wide and dense dataset range was obtained. A surface dataset for testing is also needed, and it is important that it is different from, but still within, the same parameter space as the training set. Notice that the validation set is a subset of the training set. The training and test sets, for a pair of parameters (k and H), are schematically illustrated in Figure 3. As the figure shows, the parameters in the test set are shifted a half step to be placed in the void of the training set. This ensures that the test set is located at the maximum Euclidean distance to the training set. The parameter space for the training range is specified in the table shown to the right in the figure. The training set contains 1,536 unique surfaces and the test set contains 486 unique surfaces. From the training set, 20% of the surfaces are transferred to a validation set, which is used to detect overfitting during training.

Preprocessing and Contact Mechanics
The areal roughness parameters in Table 1 are calculated for all of the 2,022 surfaces, which are made dimensionless by scaling to exhibit unit rms roughness, that is, S q 1. This is done, in connection with the non-dimensionalization of the contact FIGURE 1 | Surface topography generation scheme.
FIGURE 2 | A randomized surface, generated by following the scheme in Figure 1, using the settings presented to the right. The dimensionless areal roughness parameters (for the scaled surfaces) in Table 1 are calculated according to ISO Central Secretary (2012). They are grouped as parameters of field type and of Bearing Area Curve (BAC) type. These parameters are the input for the ANN described in Section 2.3, with architecture illustrated in Figure 4. For more details on how to calculate the areal roughness parameters according to the standard, see Blateyron (2013). As a result of the contact mechanics simulations, the surfaces may be plastically deformed, and it is the hardness p p that limits the maximum pressure the surface can exhibit before it yields plastically. Therefore, four equidistantly spaced values for p p in the range [20, 100] were used for the training set and three different values were chosen for the test set, in the same way as described for H and k in Figure 3. Thereby, there are in total 7,602 contact mechanics calculations performed. Out of these, 6,144 were used to train the ANN and 1,458 were used for testing.
The output from the contact mechanics calculations are the in situ, areal roughness parameters, and the six contact mechanics parameters in Table 2. These are the real area of contact to nominal contact area ratio A r A e + A p , the elastic-A e and plastic contributions A p to it, the dimensionless maximum contact pressure p max p max /E, the dimensionless average interfacial separation u u/h r , and the dimensionless contact stiffness K Kh r /E.

The Artificial Neural Network
Here, the ANN architecture depicted in Figure 4, which is engineered to predict the contact mechanical response of surfaces represented by the areal roughness parameters (given in Table 1), will be described. The areal roughness parameters in Table 1 and the dimensionless hardness p p are used as input for the ANN and it outputs the corresponding, in situ, areal roughness parameters and the six contact mechanics parameters in Table 2. As emphasized with double borders in Figure 4, the ANN consists of five subnetworks. These all have four fully connected layers, but a different amount of neurons per layer. The arrows with continuous lines indicate connections that are fully connected with weights, whereas arrows with dashed lines indicate just passing the data from one part of the network to another. A regular MSE loss function was adopted for the training procedure.
The first subnetwork (1) , with 64 neurons, has rectifier (ReLU) activation functions (f (x) max(0, x)), and it takes the 14 areal roughness parameters in Table 1 and the value of the hardness as input. The purpose of this network is to process the areal roughness parameters, without the influence of the contact pressure, to extract suitable input for the second subnetwork (2) with 128 neurons. This subnetwork has sigmoid activation functions (f (x) 1/(1 + e − x )), and it is fully connected to the 14 areal roughness parameters, the hardness, the dimensionless contact pressure, and the output of the first subnetwork.
The input to the first and second subnetworks and their output are assembled into one vector with 64 values. This vector is the input to each of the three parallel subnetworks (3)-(5) , which have softplus (smooth rectifying)   activation functions (f (x) ln(1 + e x )) and 256 nodes each. As mentioned previously, the purpose of the ANN is to predict the 14 in situ areal roughness parameters as well as the six contact mechanics parameters in Table 2. To this end, each of the three parallel subnetworks output parameters grouped by its origin, that is, the five in situ areal roughness parameters of field type, the nine in situ areal roughness parameters related to the bearing area (or Abbott-Firestone) curve, and the six contact mechanics parameters.

RESULTS AND DISCUSSION
First, in Section 3.1, the performance of the ANN model will be evaluated with linear regression between predicted values and the output in the test dataset. Then, in Section 3.2, examples of how the predictions changes with the load will be presented and compared to the correct values.

Predicting Contact Mechanical Response
Herein, the test set, which contains 1,458 specimens that it has never seen before, is used to evaluate the ANN's predictive performance on surfaces for a whole range of loads. Depicted in Figures 5, 6 are linear regression of all the predicted parameter values and the R 2 -value, that is, where y is the target output, y is the predicted output, and y is the mean target output, is used as a measure of the accuracy. Overall, one can see that some parameters are predicted with extraordinary high accuracy, whereas a few are predicted with less precision. Figure 5 reveals that there is a systematic error for the predictions of S a and S q , which both have relatively low R 2 -values. The reason for the low R 2 -values is because the absolute majority of predictions (for both S a and S q ) are underestimated. Among the output shown in Figure 5, the one with highest accuracy is the mean quadratic slope parameter S dq . Visually, the bearing area curve parameter S mr1 shows a quite large spread, while the R 2 -value is rather high. This is caused by a relatively small percentage predictions with large errors.
There is much that can be said about the results presented in Figure 6. One thing, which is nearly impossible not to notice, is the regression for the dimensionless maximum pressure p max , with wide-spread and a low R 2 -value. The reason for this is that the p max is a local event, while areal roughness parameters are averaged in some sense. In other words, prediction of a local quantity based on average terms is a complicated task, and the low accuracy is, therefore, to be expected. The more important outputs A r , A re , A rp , u, and K are, fortunately, predicted with higher accuracy. For more details on the ANN's predictability, the reader is referred to next section.

Application
In this section, the accuracy of the predictions of the ANN for three different test specimens, taken from the test set that the network never has seen before, will be investigated over the whole range of loads considered when the test set was generated with the contact mechanics simulations. The test specimens are listed in Tables 3, 4, and they consist of the areal roughness parameters, corresponding to the surfaces topographies presented in Figure 7, combined with a value of the dimensionless hardness.
The predictions of the in situ dimensionless mean square height S q and skewness S sk are depicted in Figure 8, which also shows the correct values. This figure and Figures 9-12 share the same legend in which the lines (continuous blue, dashed red, and dotted turquoise) are for the predictions, and the correct values are represented with the markers (round blue, round red, and cross turquoise). The ANN predicts both the S q and S sk output parameters for Specimen 1 (bi-Gaussian and self-affine) with best accuracy. The lowest accuracy was observed when predicting S q for Specimen 2 (Gaussian and self-affine), and the lowest accuracy was observed when predicting S sk for Specimen 3 (Weibull and exponential). The predictions of the in situ dimensionless kurtosis S ku and mean quadratic slope S dq are depicted in Figure 9. For S ku , to the left in Figure 9, the accuracy is rather high and approximately the same for all three specimens. The variation of the in situ kurtosis is very complex, but still accurately captured by the ANN. To the right in Figure 9, it can be seen that the predictions of the in situ S dq are of very high accuracy, and this can be understood from the linear regression analysis presented in Figure 5.   Figure 7 and a value of the dimensionless hardness: part 1-field type parameters.  Figure 7 and a value of the dimensionless hardness: part 2-BAC type parameters and dimensionless hardness. The predictions of the dimensionless average interfacial separation u and real area of contact ratio A r are depicted in Figure 10. It is observed that the ANN very accurately predicts u for Specimen 2 over the whole range of loads tested. The accuracy for Specimen 1 is not so high at low loads but really good for moderate and high loads, and it is vice versa for Specimen 3. Overall, the ANN's accuracy in predicting u is good, which is required if the ANN would be employed in a mixed lubrication model like the one in Sahlin et al. (2010). As displayed in the right of Figure 10, the real area of contact ratio can be predicted with satisfactory accuracy for all but the lowest load, where it ideally should extrapolate A r → 0 as p c → 0. Better performance could (most likely) have been obtained, by extending the training set to include more results for lower loads. Note that this would also require a higher mesh density than the 256 × 256 used presently.
The ANN is trained such that the areal roughness parameters for p c 0 remains unchanged. The ANN is also trained such that the contact mechanics parameters are zero for p c 0, except for the average interfacial separation, which is specified as the surface's maximum peak height.
The predictions of the elastic part of the real area of contact ratio A e (left) and plastic part of the real area of contact ratio A p (right) are depicted in Figure 11. When looking at the predictions for A e , it is noticeable that there is a large error for Specimen 1; however, the other specimens are predicted with acceptable accuracy. From the results for A p presented to the right in Figure 11, it seems as if the ANN has qualitatively learned what the variation of A p would be. More precisely, that it is constant for low loads but that it starts to increase at some point. The reason for that it does not quantitatively capture the variation FIGURE 8 | Dimensionless root mean square height S q (left) and skewness S sk (right) for varying dimensionless nominal pressure p c , predicted (line) and real value (marker). correctly has probably to do with that a relative error for a large A p is much more significant than it is for a small A p , during the training procedure. Notice that the ANN predicts that Specimen 3 exhibits plastic deformation, but that the correct result is that the deformations are purely elastic for all loads considered (A p 0 is not displayed on the log-scaled axis). The predictions of the dimensionless maximum pressure p max (left) and contact stiffness K (right) are depicted in Figure 12. When looking at the predictions for p max , one can see that there is a quite large error. This was also brought up in connection to the presentation of Figure 6. Most surfaces will be plastically deformed, and it seems as it would be fairly easy for the ANN to learn that the p max will saturate at p p . Indeed, by looking at the predictions for Specimen 1 and 2, it is also clear that it has learned this. Specimen 1 that has the lowest p p is already plastically deformed at the smallest load in the range and Specimen 3 with the highest p p is not plastically deformed at all. Specimen 2 does, however, exhibit p max for an intermediate load in the range, and it can be observed that the ANN is able to predict that it will and that it is capable of capturing the position where it occurs. From Figure 12, it can also be observed that the contact stiffness, for all three specimens, can be predicted with quite high accuracy for moderate and high loads but that the accuracy decreases for lower loads.

CONCLUDING REMARKS
Two datasets containing a total of 2,022 different surface topographies were generated using the algorithm developed by Pérez-Ràfols and . Three different HPD functions and two different PS functions were obtained from FIGURE 10 | Dimensionless average interfacial separation u (left) and real area of contact ratio A r (right) for varying dimensionless nominal pressure p c , predicted (line) and real value (marker).
FIGURE 11 | Real area of elastic contact ratio A e (left) and real area of plastic contact ratio A p (right) for varying dimensionless nominal pressure p c , predicted (line) and real value (marker). Note that the correct values for A p for Specimen 3 is zero for all loads.
Frontiers in Mechanical Engineering | www.frontiersin.org May 2021 | Volume 6 | Article 579825 Gaussian, bi-Gaussian, and Weibull HPD functions and selfaffine and exponential PS functions, described with as few shapedefining parameters as possible. Fourteen areal roughness parameters were calculated for all surfaces in the dataset. Together with the surface indentation hardness and a given applied load, these 14 areal roughness parameters were used as input for the ANN. A numerical elastoplastic contact mechanics approach, in which the hardness limits the maximum pressure the surface can exhibit before it yields plastically, was then employed to perform simulations of pressing each of the generated surfaces against a flat rigid counter surface for a sequence of loads. Since four values for the hardness were considered for the training set with 1,536 different surface topographies and three were considered for the test set with 486 topographies, a grand total of 4 × 1536 + 3 × 486 7602 realizations were conducted. Out of these, 6,144 specimens were used for training and 1,458 were left for testing and validation. For each of the these specimens, 14 in situ areal roughness parameters and six contact mechanics parameters were calculated for the sequence of loads that was also used as input for the ANN.
An architecture for an artificial neural network (ANN), which consisted of five different subnetworks, was designed and trained on the dataset. Linear regression was applied, and the R 2 -value was used to appreciate the correlation between the network prediction and the correct data. A few parameters were almost perfectly predicted, whereas other were predicted with large errors. According to the R 2 -values, the most important parameters, that is, the real area of contact ratio A r , the dimensionless average interfacial separation u, and contact stiffness K were all predicted accurately by the ANN.
Summing up, the ANN can be used to roughly appreciate the in situ behavior of various kinds of surface topographies, if the areal roughness parameters, the indentation hardness, and the nominal contact pressure are known. Some parameters, that is, the real area of contact ratio, the dimensionless average interfacial separation, and contact stiffness can actually be predicted with high accuracy.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
KK has performed the main part of the development of the ANNbased model and has been coordinating the writing. RL contributed with expertise, to the discussions and was engaged in the writing. FR has been engaged in surface generation and contact mechanics and contributed to the writing. ML contributed with expertise, particularly regarding AI and machine learning, contributed to discussions and to the writing. AA has initiated the work and has been involved in all parts of it.