Identification of Roughness with Optimal Contact Response with respect to Real Contact Area and Normal Stiffness
Additive manufacturing technologies are a key point of the current era of Industry 4.0, promoting the production of mechanical components via the addition of subsequent layers of material. Then, they may be also used to produce surfaces tailored to achieve a desired mechanical contact response. In this work, we develop a method to prototype profiles optimizing a suitable trade-off between two different target mechanical responses. The mechanical design problem is solved relying on both physical assumptions and optimization methods. An algorithm is proposed, exploiting an analogy between genetics and the multiscale characterization of roughness, where various length-scales are described in terms of rough profiles, named chromosomes. Finally, the proposed algorithm is tested on a representative example, and the topological and spectral features of roughness of the optimized profiles are discussed.
The role of the interface between different material constituents/phases is a main research topic in this current era of Industry 4.0, which is radically changing perspectives of industry about the technological manufacturing of components and/or materials [1, 2]. In engineering applications, one often requires the design of surface textures able to satisfy desired target responses and/or to allow rapid manufacturing and morphological changes. This can be useful, e.g., for the in-line control of mechanical components [3, 4]. Moreover, the use of additive manufacturing technologies is currently promoting a detailed and unified specific description of roughness, able to take into account several length-scales in its representation [5, 6].
Nowadays, most approaches to design roughness are based either on mimicking natural surfaces [7, 8], or on performing parametric studies via artificially generated surfaces [9, 10]. In this context, the present work proposes a method to design roughness with the aim of matching two given mechanical contact target responses as close as possible. These two targets are associated with a frictionless elastic normal contact problem. In the article, multiscale roughness is modeled by a finite series of cosinusoids which, using a biological analogy, is called the surface roughness genome. Then, the constituent waves identifying various length-scales of roughness are properly selected and mixed based on mechanical considerations and suitable optimization tools. The resulting profile achieves a mechanical response that optimizes a suitable multiplicative trade-off between the two different mechanical targets.
The article is structured as follows. In Section 2, we summarize the contact problem and we introduce the surface roughness genome. In Section 3, we propose an algorithm to generate prototype profiles with the aim of reproducing two given target mechanical contact responses as close as possible. In Section 4, the effectiveness of the proposed algorithm is tested on a representative numerical example. Finally, Section 5 concludes the paper and discusses possible extensions.
2. Contact Problem and Roughness Model
This study refers to the modeling and optimization of the mechanical behavior of two contacting rough surfaces. First, the contact problem and the roughness model are presented. Then, the mechanical interaction of different roughness length-scales is used to identify two categories of roughness. Finally, an algorithm is proposed to optimize roughness, expressing the objective function of the optimization problem as a suitable trade-off between two different target responses. To model contact, the frictionless elastic normal case is considered, taking as targets a desired stiffness-load curve, , and a desired contact area-load curve, .
2.1. Contact between Two Rough Surfaces
Following the seminal work of Johnson , under the assumption of linear elasticity, the nonconforming contact problem between two rough surfaces is equivalent to a contact problem between a rigid composite rough surface and a half-plane depending on effective elastic parameters; see also Barber  for details and proofs. Normal contact is herein controlled by imposing an approaching far-field closing displacement , whose value is computed from the tallest summit of the rough surface, whose height is denoted as . The deformation of asperities produces a normal contact traction distribution in the contact area . The average total contact pressure is computed as the sum of all the contact forces acting in this area divided by the nominal contact area. The derivative of the contact pressure with respect to the imposed displacement defines the contact stiffness .
As an example, in Figure 1, which refers to a two-dimensional case, the topography of a numerically generated rough surface is presented (a), together with its deformed configuration during contact with a half-plane (b). In the figure, denotes the value of the dimensionless displacement.
In the work, roughness is described in terms of various parameters (see Section 2.2). For each choice of the set of parameters, the contact problem is solved using the Boundary Element Method (BEM); see Bemporad and Paggi  for a comparison of various techniques which can be used to impose the satisfaction of unilateral contact constraints.
2.2. Roughness Description over Multiple Length-Scales
Referring now for simplicity to one-dimensional rough profiles, their topography is herein described in terms of the real part of the Multivariate Weierstrass-Mandelbrot function (MWM). Ausloos and Berman  proposed the MWM function to describe stochastic processes with a larger number of features than in the original 2D formulation early proposed by Mandelbrot .
A rough profile modeled via the MWM function has then a topography depending on several parameters:
In the absence of any subscript, the base 10 is used to compute the logarithm. In (1), the parameter is a positive integer, the parameters , , and the parameter . These parameters are called genes, exploiting an analogy with biology. Similarly, the expression surface roughness genome refers to the overall ensemble of genes providing the realization of a profile over several different length-scales.
The indices and are integers, which determine, respectively, the longest and the shortest wavelengths contributing to the definition of the rough profile . Such a rough profile is realized over an observation length , and sampled in nodes according to a resolution . Therefore, in this model, roughness is realized by varying the index between and , together with the other parameters. Taking into account the reference wavelength , the indices and are expressed as follows:In (2), the notation denotes the largest integer smaller than or equal to its argument. The number of terms in the inner summation in (1) is . The reader is referred to Majumbdar and Bhushan ; Wang and Komvopoulos ; Wu [18, 19] for further details on the spectral properties relating (2) to (1).
The particular combination of genes associated with each fixed value of the index identifies a rough profile, whose expression readsand which is named chromosome. This profile is associated with the features of roughness defined at the fixed reference length-scale . In this way, it is possible to describe the features of a single length-scale of roughness in terms of a chromosome and obtain the complete profile summing up chromosomes:Without loss of generality, the value is assigned to the chromosome with the longest wavelength, which refers to the coarsest realization of the profile.
2.3. Macroscale and Microscale Roughness
The nonlinear interaction of multiple length-scales of roughness provides the collective mechanical response of a rough profile. For example, the thermal/electric contact conductance is intimately related to the stiffness-load curve , as addressed by Barber . Moreover, as investigated by Paggi and Barber , the stiffness-load curve is mainly affected by the long wavelengths contributions. On the other hand, the contact area morphology is determined by the short wavelengths composing the power spectral density (PSD), as confirmed by many authors; see, e.g., Barber ; Paggi and Barber ; Yastrebov et al. .
In the following, the parameters that identify the stiffness-load curve are isolated and emphasized. Here, roughness is characterized as a result of two contributions, according to the role of each chromosome in the definition of the evolution. Considering a given observation length-scale, the associated rough profile is split into two rough contributions and :
The profile , which represents the macroscale roughness, is obtained by summing up the chromosomes whose evolution correlates with the one of the complete profile with a correlation coefficient .
The profile , which represents the microscale roughness, is associated with the remaining set of chromosomes, having .
This distinction is shown in Figure 2. The macroscale contribution is depicted with a red line and the microscale one with a blue line. Their sum gives the complete rough profile (black line). The stiffness-load curve and the area-load curve of those profiles are shown in Figure 3.
Regarding the curves, the ones provided by macroscale roughness and by the complete profile almost overlap. The same observation can be done for the evolution for small pressures. However, the two curves diverge for N/m.
From the operative stand point, chromosomes associated with macro-roughness are selected as follows. First, the evolution of a rough profile is computed by solving the frictionless normal elastic contact problem with imposed displacements in its peak-valley amplitude. To solve such a problem, the Nonnegative Least Squares method proposed in Bemporad and Paggi  is used. Then, the evolution is computed for each of the chromosomes in the power density spectrum of the profile, and the correlation coefficient between and the curve associated with the complete profile is computed. Chromosomes with , i.e., the ones leading to macroscale roughness, are included in a set of chromosomes . The remaining chromosomes are associated with microscale roughness.
3. Optimal Design of a Rough Profile
In the following, we consider the problem of designing a rough profile in such a way that its stiffness-load curve and contact area-load curve are close, respectively, to a given target stiffness-load curve , and a given target contact area-load curve . To do that, we propose an algorithm, called Mixed Chromosomes Cross-Over (M-CCO). Loosely speaking, starting from a given genome database, the algorithm mixes chromosomes belonging to genomes in the database leading to macro-roughness (which make it possible to achieve the target) with chromosomes belonging to genomes in the database leading to micro-roughness (which make it possible to achieve also the target). Then, their parameters (genes) are optimized to refine and further improve matching with the target contact responses. In Section 3.1, we describe the M-CCO algorithm. Then, in Section 3.2, we report some details on another algorithm, in this case taken from the literature, which is exploited in an inner loop of the proposed M-CCO algorithm.
3.1. Proposed Algorithm: Mixed Chromosomes Cross-Over (M-CCO)
Before going into the details of the proposed M-CCO algorithm, we introduce some notations and concepts. First, we define the similarity score which quantifies, for each index , how much the target mechanical response, , and the one associated with the -th rough profile (genome) in the database are similar. For , the target and the -th response are, respectively, and (where the index refers here to the -th genome in the database). For , they are, respectively, and . In (6), denotes the -norm, computed on the points used to discretize each mechanical response. The -th mechanical response coincides with the target one when holds.
To compute consistently the similarity scores in (6), we have to assume that the two mechanical responses are evaluated on the same range of pressures (from 0 to the same maximum pressure), i.e., that . To do that, for each genome in the database, its amplitude gene is rescaled aswhere the step in (7) is possible because the gap function in the BEM formulation scales linearly with the height field, and so does the pressure; see Johnson . The maximum pressure level is achieved by solving the contact problem with a far-field displacement equal to the target profile amplitude. The same remark holds for the maximum pressure level of the -th genome. Moreover, all the genomes in the database refer to the same reference length .
The steps of the M-CCO algorithm are reported in Algorithm 1 and are described as follows. Firstly, the similarity score in (6) is computed separately between each target ( or ) and the corresponding mechanical response ( or ) associated with each genome in the database, obtaining the vectors and with components and (Steps )-().
Then, the macroscale roughness is firstly identified, finding all the genomes with a similarity score larger than . At this point, these genomes are reduced, limiting to their chromosomes leading to macro-roughness. In other words, for each -th selected genome, only the chromosomes that affect the target significantly (i.e., the ones for which their stiffness-load response and have correlation coefficients larger than 0.95) are kept. The set of reduced genomes obtained in this way is named (Step ).
Similarly, the genomes with a similarity score larger than are determined. Then, they are reduced, keeping only their chromosomes leading to micro-roughness. In other words, for each -th selected genome, only the chromosomes that affect the target significantly (i.e., the ones for which their contact area-load response and have correlation coefficients larger than 0.95) are kept. The set of reduced genomes so obtained is named (Step )).
The resulting sets and contain information about different sets of chromosomes. All possible new genomes obtained by summing a reduced genome from the set and one from are stored in the set of genomes (Step ). For each such genome , the contact mechanics problem is solved via BEM, determining its and curves (Step (). Moreover, its amplitude gene is rescaled in order to satisfy the maximum pressure requirement, according to (7).
For each genome , two similarity scores are now computed, and . They refer, respectively, to the similarity of each target response ( or ) and the related mechanical response curve ( or ) associated with that genome. A mixed similarity score is also computed, by multiplying these two similarity scores, i.e., (Steps )-().
Now, the three genomes in with the largest mixed similarity score are stored in the set (Step (14)). Then, for each genome in , the similarity with the target function (the mixed similarity score) is increased by optimizing the values of the genes belonging to the chromosomes related to macro-roughness. As an alternative, the weighted sum of two objectives (i.e., of the squares of the two similarity scores) could be optimized, as done in [22, 23] for different problems. This approach resembles some regularization techniques applied in machine learning problems; see, e.g., Gnecco and Sanguineti [24, 25]; Gnecco et al. . However, for such techniques, the optimal selection of the regularization parameter may not be straightforward. Nevertheless, the weighted product of the two objectives can be a valid alternative. See Gnecco et al. ; Gnecco  for some numerical examples.
The square of the mixed similarly score is maximized by applying the Globally Convergent Method of Moving Asymptotes (GCMMA) algorithm [28, 29] (Step (16)); see Section 3.2 for further details. Then, starting from the optimized genome, the GCMMA is applied to its remaining chromosomes, which are associated with microroughness (Step (17)), considering the same objective function of Step (11).
Finally, the new genome with the maximum obtained value of the similarity score is identified (Step (19)).
3.2. GCMMA and Its Use inside the M-CCO Algorithm
The Globally Convergent Method of Moving Asymptotes (GCMMA) [28, 29] is an algorithm to solve nonlinearly constrained optimization problems. It replaces each such problem with a sequence of approximating optimization subproblems. These are easier to solve, and are still nonlinearly constrained.
The term moving asymptotes refers to the functions which are used in the formulations of the approximating subproblems. Such asymptotes usually change from one iteration of the GCMMA to the successive one. The GCMMA is globally convergent in the sense that, for each possible initialization of the vector of parameters to optimize, the algorithm is proved to be convergent to a stationary point of the initial problem.
In Steps (16) and (17) of Algorithm 1, the GCMMA is applied by taking as objective function not the mixed similarity score , but its square. This is done in order to get a smoother objective function. The GCMMA iterative solution is obtained after a number of steps, starting from an initial choice for the vector of parameters to optimize. For the simulations, has been chosen, in order to keep the computational time per iteration small.
To apply the GCMMA algorithm inside Algorithm 1, a MATLAB implementation has been made, based on the MATLAB functions mmasub.m and subsolv.m, developed by Svanberg himself. Technical details on the implementation of these two MATLAB functions are provided in Svanberg . At each step of the GCMMA, the partial derivatives of the objective function with respect to all the optimization variables are computed numerically, via a small perturbation for each such variable. Moreover, the amplitude parameter is rescaled according to Eq. (7).
The GCMMA algorithm is applied with several different initializations, in order to increase the probability of obtaining a good constrained local maximizer of the original objective function. Finally, the best solution, i.e., the one associated with the greatest (square of the) mixed similarity score obtained during the various optimizations, is produced as output, for both Steps (16) and (17) of Algorithm 1.
To save computational time, the number of variables in each optimization problem is reduced as described in the following. The parameters (genes) , and determine the frequency spectrum, i.e., the interaction among several length-scales. Hence, they are not modified during the optimization, but they are fixed to their original values. The parameter is rescaled at each optimization step in order to satisfy the pressure requirement, i.e., applying Eq. (7). So, only the phases are considered as parameters to optimize. In the optimization, their values are constrained in the range between of their initial values, in order to preserve the main features of the original chromosomes. Moreover, for both Steps (16) and (17) of Algorithm 1, only the phase genes associated with the subset of chromosomes related, respectively, with macro-roughness and micro-roughness, are optimized.
4. Numerical Experiments
The M-CCO algorithm is now tested, designing roughness in a realization length µm. First, an artificial genome database is generated to emulate a priori known database of real rough profiles (see Section 4.1). The two targets and are therefore selected not to match any of the mechanical responses associated with one of the rough profiles in the database. The two selected targets are shown as red lines in Figures 4(a) and 4(b), respectively. Numerical results about the application of the M-CCO algorithm to this specific case are reported in Section 4.2.
4.1. Genome Database
For all the profiles in the database, the amplitude gene is fixed to . Moreover, the main wavelength is set equal to µm. Twenty different pairs of values for and have been generated according to a Sobol’ sequence . To generate the phase matrix , has been imposed for all the genomes. Three phase vectors of thirty elements each have been generated, with values from 0 to for their components. Such values have been also determined using a Sobol’ sequence. Combining the pairs with the phase vectors, a database made of genomes has been obtained.
Each profile in the database has been discretized using nodes. To formulate and solve each frictionless normal contact problem, the elastic modulus has been set equal to MPa. Then, the mechanical response has been computed in equally displaced far-field displacements starting from the tallest summit of the rough profile, proceeding likewise in Johnson ; Bemporad and Paggi .
4.2. Numerical Results
The three best genomes obtained from the M-CCO algorithm are summarized in Table 1. The table reports, for each such -th genome, its mixed similarity score , and the similarity scores and of its mechanical responses with respect to the two targets. Also, it reports the genomes used in the initialization of the GCMMA algorithm in Steps (16) and (17) of Algorithm 1. Such genomes are used to represent, respectively, macro- and micro-roughness. The first index refers to the pair (see Figure 5(a)), whereas the second one to the vector of phases (see Figure 5(b)).
The best approximation of the target contact responses is provided by the genome , whose mechanical responses and are shown by blue curves with round markers in Figures 6(a) and 6(b), respectively. Looking at the contact area-load curves reported in Figure 6(b), one can notice that the mechanical responses of all these three solutions diverge consistently from the target curve in the range of contact pressures between N/m and N/m.
As far as the stiffness-load curve is concerned (see Figure 6(a)), the and curves show a very good overlap with the target . The second solution is the only one associated with a different genome representing macro-roughness, and its trend is not close to in the range of contact pressures between N/m and N/m.
The curve almost overlaps for N/m. Indeed, the genomes and are composed by the same starting sets of chromosomes representing macro-roughness. Also the curves and are very similar (see Figure 6(b)), and quite close to the target , apart from an intermediate range of values for the contact pressure.
For both genomes and , the same pair is used to represent micro-roughness, but a different vector of phases is used in the GCMMA initialization in Step (17) of Algorithm 1. This highlights the dominant role of macro-roughness in the frictionless elastic normal contact response.
The three rough profiles obtained from the M-CCO algorithm are visualized in Figure 7, which emphasizes the roles of macro-roughness (see Figure 7(a)), and micro-roughness (see Figure 7(b)). Each rough profile in Figure 7(c) is obtained by summing these component profiles (represented by the same colors and line types).
In Figure 7(a), the macro-roughness of each solution is depicted. The dashed blue and the dash-dotted red profiles, which are associated, respectively, with the genomes and , have a similar macro-roughness topography, with small differences, e.g., in the range between N/m and N/m.
As observed before, these solutions are originated from the same set of chromosomes as regards the macro-roughness level. Moreover, they have also quite close values of similarity scores (see Table 1).
The micro-roughness of genomes and is shown in Figure 7(b). The figure shows that they have similar micro-roughness, as they are associated with the same pair (). The resulting topographies of the genomes and are shown in Figure 7(c) and are also quite similar. Moreover, these two genomes have quite close mechanical responses (see Figure 6).
The topography of genome is quite different from the ones of and , since its macro-roughness component differs significantly.
Finally, the spectral decomposition of the three genomes is shown in Figure 8 (in the figure, has the same dimension as ). The difference in amplitude is due to different values for the amplitude gene , which are obtained by applying Eq. (7) with respect to the threshold pressure.
All the results of Fast Fourier Transform (FFT) filtering show a power spectral density of the complete profile containing several peaks. For all the three genomes, the first part of the PSD does not present any evident peak. For the genome , the FFT filtering peaks are close to the discrete spectrum of the genome (see the squares in the central part of the spectrum in Figure 8(b)). Moreover, the squares in Figure 8(b), which represent the micro-roughness, are six, whereas the peaks observed from the FFT filtering are five.
Looking at the solutions and , which show similar mechanical responses and topographies, one can notice that they present also similar PSDs in the low-frequency part of the spectrum. However, their high-frequency behaviors are different, because they derive from different genomes.
In this paper, the Mixed Chromosomes Cross-Over (M-CCO) algorithm has been proposed to optimize a roughness topography to match two different target mechanical contact responses, namely the normal contact stiffness and the real contact area ; both are expressed as functions of the average contact pressure . The algorithm has been inspired by an analogy with genetics, according to which the M-CCO mixes the macro-roughness and the micro-roughness of two different “genomes”. Then, an application of the algorithm to a representative numerical example has been described, showing very promising results. One can conclude that two particular classes of genomes can be identified and combined to match both the and targets. Particularly, chromosomes associated with macro-roughness are fundamental for an optimal design.
Still, the results show that there is room for possible improvement, which could be obtained, e.g., by future variations of the M-CCO algorithm. Particularly, the algorithm could be improved in the selection of the chromosomes to combine. Future work will also concern the introduction of additional interactions between the rough surfaces in contact, such as friction, adhesion, wear, and so on. Such mechanical interactions are difficult to predict using BEM, and it is probably necessary to move to the Finite Element Method (FEM) to tackle nonlinear multi-field coupled contact and fracture problems with complex and realistic geometries. For example, the case of elasto-plastic contact is certainly a fundamental aspect to be investigated, and plastic deformation can be taken into account by using ad hoc computational models.
This work is a first step in the generation of surface topographies achieving desired target mechanical responses. The problem investigated in the article and the M-CCO algorithm are potentially useful for the production of surface pressure sensors in sealing applications. In this case, the contact has to be assured for a contact pressure larger than some threshold. At the same time, the electrical conductivity is associated with the curve. Then, the curve might change its trend suddenly above the contact pressure threshold. Finally, the curve may be chosen in such a way to improve sealing.
Another important potential application refers to the problem of surface morphing. In this case, roughness may be modified in time based on external stimuli. Morphing is a quite novel concept, which has been recently applied to smart structures and devices. However, its application to surfaces is a very new research area.
The problem considered in the work is important also in view of the fact that, nowadays, additive manufacturing technologies simplify the realization of topographies, enabling the production of surface prototypes. 3D printing is one of the most fascinating and robust technology, due to its very fast progress and impact on several industrial sectors, which promotes technology transfer and patent applications. Applications of the proposed methodology to surface printing could open new research areas.
The MATLAB data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
M. Paggi acknowledges the support of MIUR-DAAD to the Joint Mobility Program 2017 “Multi-Scale Modeling of Friction for Large Scale Engineering Problems”. G. Gnecco acknowledges the support of FFABR (Fondo per il Finanziamento delle Attività Base di Ricerca) from the Italian Ministry of Education, University and Research (MIUR).
M. Brettel, F. Friederichsen, M. Keller, and M. Rosenberg, “How virtualization , decentralization and network building change the manufacturing landscape: an Industry 4.0 perspective,” International Journal of Mechanical, Aerospace, Industrial and Mechatronics Engineering, vol. 8, pp. 37–44, 2014.View at: Google Scholar
F. Borodich and O. Savencu, “Hierarchical models of engineering rough surfaces and bio-inspired adhesives,” in Bio-inspired Structured Adhesives: Biological Prototypes, Fabrication, Tribological Properties, Contact Mechanics, and Novel Concepts, L. Heepe, L. Xue, and S. Gorb, Eds., vol. 9 of Biologically-Inspired Systems, pp. 179–219, Springer, Cham, Switzerland, 2016.View at: Google Scholar
K. Cummins, “The rise of additive manufacturing. The Engineer,” 2010, https://www.theengineer.co.uk/issues/24-may-2010/the-rise-of-additive-manufacturing/.View at: Google Scholar
M. Sherge and S. Gorb, Biological Micro- and Nano- Tribology & Nature’s Solutions, Springer, Berlin, Germany, 2001.
M. Nosonovsky and B. Bushan, Multiscale Dissipative Mechanisms and Hierarchical Surfaces: Friction, Superhydrophobicity, and Biomimetics, Springer, Hoboken, NJ, USA, 2008.
K. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, UK, 2003.
B. B. Mandelbrot, Fractals: Form, Change and Dimension, W. H. Freeman, San Francisco, Calif, USA, 1977.View at: MathSciNet
K. Svanberg, “MMA and GCMMA two methods for nonlinear optimization,” Technical Report, KTH, Stockholm, Sweden, 2007.View at: Google Scholar
H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, SIAM, Philadelphia, PA, USA, 1992.