In this work we assess the applicability of the adjoint optimization technique for determining optimal surface topogra- phies of two surfaces in relative motion in presence of a thin lubricant ﬁlms that can cavitate. Among the existing numerical tools for topology optimization in engineering problems, the adjoint method represents a promising and versatile technique, which can also be applied to the ﬁeld of full ﬁlm tribology. In particular, the design of surfaces with complex textures can thoroughly beneﬁt from this method, as it allows dealing with a large number of degrees of freedom at low computational cost. We show that this optimization method can be successfully applied to cavitating lubricant ﬂows such as in pin-on-disc tribometers, giving the possibility to extend the results also to other typical applications such as journal and slider bearings. It is shown that the adjoint method can optimize the whole gap height distribution point by point in a more e cient way than traditional optimization approaches and parametric studies. In particular, thanks to the sensitivity analysis the adjoint method is able to ﬁnd the placement and depth proﬁle of each texture element. Keywords: surface textures, adjoint method, optimization, mass conserving cavitation 1. Introduction Improvements in the e ciency and reliability of tri- bological devices are often achieved through shape opti- mization of the sliding contacts either by modiﬁcations of the global contact geometry ( e.g.optimized Rayleigh step or slider bearing [1, 2]) or locally by the introduc- tion of surface textures [3]. The typical goal of the opti- mization process in hydro-dynamically lubricated con- tacts is to increase the separation between the sliding surfaces in order to reduce friction and wear. This goal is achieved by increasing the hydrodynamic pressure generated in the contact through the shape optimization of the surfaces. Shape optimization opens the way to potentially inﬁnite variations of surface modiﬁcations. For this reason, di erent optimization techniques have been proposed in literature, which can be categorized in two main groups: optimization techniques based on parametric studies and optimization techniques based on other more sophisticated optimization approaches. Email address: andrea.codrignani@imtek.uni-freiburg.de (A. Codrignani)The ﬁrst group comprises the vast majority of liter- ature, since the shape of the geometry modiﬁcations is often parameterized by using elementary shapes, both for textured or stepped surfaces [4]. This is because the complexity of a free-form optimization of the whole gap height distribution h(~x) is usually reduced through ele- mentary shapes or piecewise approximation of the in- vestigated geometry. For example Rohde [5] as well as Shen & Khonsari [6] maximized the load carrying ca- pacity of a stepped surface for a 2D slider bearing by parametrizing the width, the depth and the position of a Rayleigh step. Many other researchers focused on parametric stud- ies, in particular for the analysis of surface textures. Since the early works of Etsion and coworkers [3, 7] the parametric approach has been the widely used for the optimization of textures in several kind of applica- tion, such as 1D slider bearings [8, 9], 2D slider bear- ings [10], Piston-Liner Contacts [11, 12], thrust bear- ings [13] and pin-on-disc tribometers [14]. In some cases, the optimization of surface texture for 1D slider bearings has been also performed by means of analyt- Preprint submitted to Tribology International February 27, 2020 ical solutions [15, 16, 17]. It is important to notice that most of the studies on both textured and stepped surfaces always proved the well known theoretical re- sult that the stepped geometry ( i.e. a redesign of the global surface /contact geometry) performs better than any other textured case. In particular the work of Fillon and Dobrica [18] states the importance of texturing in the same region of the bearing where a Rayleigh step should be. In the second group more complex optimization strategies have been employed. For example Guzek et al. optimized a 2D parallel bearing with the aid of quadratic programming, hence exploiting the gradient of the cost function in order to reach the optimal design faster than with usual trial and error approaches [19]. Their method, however, focused on the optimization of predetermined texture topology such as rectangular or elliptical dimples. For this reason other methods based on the optimization of arbitrary shapes have been pro- posed. Among them Fesanghary and Khonsari [20] ap- plied a hybrid optimization algorithm based on sequen- tial quadratic programming and harmony search to a slider bearing with groove shapes approximated through splines. As a result, a higher load carrying capacity was obtained than through elementary groove shapes. This hints to the fact that an increase in the number of pa- rameters in the optimization algorithm, which in the ex- treme is equivalent to a point by point optimization of the whole geometry, will perform better than an opti- mization algorithm based on a reduced number of pa- rameters, such as an optimization of the surface topol- ogy described by a set of elementary shapes. Similar conclusions were found also in the work of Wang et al. [21] where elliptical air bearings are optimized by mod- ifying the whole surface in a continuous way thanks to an automated mesh generation. Other strategies for the free-form optimization of the gap height have been pro- posed by Maday using calculus of variations in the case of a 1D journal bearing [22, 23]. Beside the optimiza- tion based on the simple Reynolds equation, sequen- tial quadratic programming was also used by Alyaqout and Elsharkawy [24] for slider bearings in the frame of thermo-hydrodynamic lubrication by the coupling of the Reynolds and energy equations. More generally, optimization strategies share the common concept of sensitivity analysis, that is, look- ing for the direct relationship between variations in the design parameter space and the corresponding changes in the object function ( i.e.the function which describes the performance to be achieved). Sensitivity analysis is furthermore at the heart of optimization methods which fall under the category of Lagrangian based methodsand adjoint methods [25]. The main advantage of such methods is the possibility of optimizing the whole gap height distribution h(~x) without assigning any predeter- mined topology modiﬁcation. This is in contrast to what is done in a typical optimization study based on a para- metric geometry. The adjoint method takes advantage of the La- grangian formulation in order to deﬁne a second partial di erential equation which is referred to as the adjoint equation. The solution of the adjoint equation allows to compute the sensitivity in an e cient way, hence know- ing the impact of modiﬁcations of every single grid point on the objective performance. The adjoint method was ﬁrst introduced in control theory [26] and then suc- cessfully applied in the aerospace sector for the opti- mization of complex shapes under the fulﬁllment of the Navier-Stokes equations [27, 28, 29]. The ﬁrst usage of the adjoint method for tribological application was per- formed by Yoon and Choi [30] for the optimization of an air bearing for hard disks showing a huge reduction in computational cost compared to traditional paramet- ric approaches. Subsequently, the adjoint method was employed for the ﬁrst time by Van Ostayen [31, 32] for the optimization of the load carrying capacity of hydro- dynamic slider bearings with particular attention also to the shape optimization under unsteady conditions. In the current work we present an extension of the work by Van Ostayen and employ the adjoint opti- mization for the free-form optimization of both the macroscopic geometry and surface textures. To the au- thors’ knowledge, the optimization is carried out for the ﬁrst time under the constraint of the incompress- ible Reynolds equation with mass-conserving cavita- tion. Both the continuous and the discrete formulation of the adjoint method have been implemented under two di erent computational environments (ComsolR [33] and MatlabR ). The validation of the optimization algo- rithm is ﬁrst performed against analytical 1D solutions and subsequently between the two di erent implemen- tations because of the lack of analytical solutions in the presence of cavitation. In order to ease the occurrence of cavitation we choose to test the optimization algo- rithm on the geometry of a pin-on-disc tribometer. This choice has the additional advantage to provide useful insights for future experimental activities. The goal of the optimization is to achieve the max- imum ﬂy height for a given load since this results in an increase of the gap height and hence in a reduction of friction by bringing the transition from mixed lubri- cation to hydrodynamic lubrication to a lower sliding speed, as pointed out in several studies [6, 34]. More- over the adjoint method allows to freely ﬁnd the optimal 2 depth distribution inside the texture elements, giving the possibility to compute the optimal shape of texture ele- ments according to their position on the macro geometry in a much more e cient way than traditional methods [35]. 2. Numerical approach The present optimization problem deals with the thin lubricant ﬁlm between two surfaces in relative motion. The lubricant ﬂow is described by a set of Partial Dif- ferential Equation (PDE) which describes the mass and momentum conservation of a lubricant ﬂuid subjected to cavitation. The optimization problem treats the PDE set as a constraint. In this section we present the for- mulation of the adjoint method applied to the Reynolds equation with mass-conserving cavitation. 2.1. Governing equation Under the assumption of an incompressible Newto- nian ﬂuid and a su ciently regular gap height distribu- tionh(~x), the shear ﬂow between two sliding walls can be modeled through the incompressible Reynolds equa- tion in the following form [36]: g( ¯p;;h)r" h3 12r¯p+1 2~Vh(1)# =0 (1) Here, ¯ p(~x)=p(~x)pcavis the pressure relative to the cavitation pressure pcav,is the dynamic viscosity, ~V=fU;WgTis the upper wall velocity and (~x) is the cavity fraction. As proposed by Elrod and Adams [37] the cavity fraction can be expressed as a function of the local density (~x) of the ﬂow (~x)=1(~x) ref: (2) whererefdescribes the reference density of the lubri- cant. In this way the cavity fraction is deﬁned to be =0 where the ﬂow is completely in the liquid state and has a maximum value of =1 when the ﬂow is com- pletely in the vapor phase. On the other side, the relative pressure ¯ pis zero in the cavitation region and positive elsewhere. These two observations can be combined in a single condition which expresses the complementarity constraint of equation 1: (ppcav)=¯p=0: (3) Instead of solving equations 1 and 3 in an iterative way, Woloszynski et al. proposed an algorithm which di- rectly couples both equations in a single system [36].Besides ensuring mass conservation, the coupled solu- tion also has the advantage of being noticeably faster than previous iterative methods [38, 39, 40]. The ﬁrst step in the derivation of such an algorithm is to rewrite the constraint equation 3 using the Fischer- Burmeister equation, so that: f( ¯p;)¯p+q ¯p2+2=0: (4) After the discretization, the system of equations fandg can be combined in a single system, whose non-linearity can be addressed by applying the Newton-Raphson al- gorithm. The system of equations 1 and 4 is completed by assigning the pressure and cavity fraction values on the Dirichlet boundaries Dand by assigning the ﬂux ˆn[h3 12r¯p+1 2~Vh(1)]=qon the Neumann bound- aries N. The resulting normal force, which is also re- ferred to as the load carrying capacity, is given by: W=Z p(~x)pextd (5) where pextrepresents the external pressure acting on the other side of the wall. The tangential force at y=his computed as [36] ~FT=Z 266664h 2rp+(1)~V h377775d ; (6) where is the surface area in the xzplane (see Fig. 1). 2.2. Optimization through the adjoint method The optimization of hydrodynamic bearing has typi- cally two objectives: either the maximization of the load carrying capacity or the minimization of the friction co- ecient. The two objectives may lead to di erent op- timal shapes as highlighted among others by Rahmani et al. [41]. In this work we focus on the ﬁrst objective, since a higher load carrying capacity can also reduce friction indirectly due to the increase in the overall gap height [2]. For this reason we can formulate our opti- mization problem as follows: maximize h(~x)fobj=W=Z p(~x)pextd subject to g(p;;h)=0 and f(p;;h)=0 h(~x)>h0 for a given~V;;pcav:(7) where h(~x) is the generic gap height distribution and h0 is the minimal gap height as shown in ﬁgure 1. The 3 Figure 1: Sketch of the generic geometry and the coordinate system. optimization is constrained by the fulﬁllment of a par- tial di erential equation i.e.the Reynolds equation with cavitation, where pandare the unknowns and his the optimization parameter. The objective is to ﬁnd the val- ues of the variable ﬁeld h(~x) for all~xwhich optimize an objective function of the form fobj( ¯p;;h) under the constraint g( ¯p;;h)=0 and f( ¯p;;h)=0. This class of optimization problems can be more con- veniently rearranged as a minimization problem of a La- grangian functional in the following form: L( ¯p;; p;;h)=fobj( ¯p;;h) +Z  pg( ¯p;;h)+f( ¯p;;h) d (8) wherep(~x) is the adjoint pressure and (~x) the adjoint cavity fraction. Since the constraints g( ¯p;;h)=0 and f( ¯p;;h)=0 are always fulﬁlled, minimizing the ob- jective function fobjis equivalent to the minimization of the whole Lagrangian functional L( ¯p;; p;;h) in Eq. 8. Therefore, in the optimal point the gradient of the Lagrangian function is zero for all variables and we can derive the following three sets of equations: L pg( ¯p;;h)=0 L f( ¯p;;h)=09>>>>>>=>>>>>>;Reynolds equation with cavitation(9) L ¯plp( ¯p;; p;;h)=0 L l( ¯p;; p;;h)=09>>>>>=>>>>>;adjoint equation (10) L hd( ¯p;; p;;h)=0 decision equation :(11)The solution of the three sets of equations will give the optimal ﬁlm height h(~x) with the corresponding pres- sure ﬁeld ¯ p(~x) and mass fraction (~x), as well as the adjoint variable ﬁeld p(~x) and(~x). Equation 9 rep- resents the already shown Reynolds equation and the complementarity equation, while equation 10 is the set of adjoint equations which reads: 8>>>>>>><>>>>>>>:r h3 12rp! +0BBBBB@1pp p2+21CCCCCA=@fobj @¯p r 1 2h~Vp! +0BBBBB@1p p2+21CCCCCA=0: (12) Its boundary conditions are assigned in the same man- ner as for the Reynolds equation with cavitation, i.e.by assigning the value of the multipliers on Dand the ﬂux of the multipliers on N. The know term of the second line of equation 12 is zero since the objective function fobjdoes not depend on the cavity fraction . Further details of equation 12 are given in appendix. The third equation 11 is called decision equation which relates the direct variable and the adjoint variable to the variation of the gap height ¯hwhich minimizes the objective function. At this point one could follow two di erent strategies for the implementation of the opti- mization algorithm according to the derivation of the decision equation. The ﬁrst strategy consists of the derivation of the Lagrangian formulation with the objective function ex- pressed in Eq. 7. Since the objective function does not depend on the gap height, one can directly compute the decision equation 11 in the following way: d( ¯p;; p;;h)=3h2 12rprp+1 2rp~V(1)=0: (13) The inﬂuence of the Lagrangian multiplier of the cav- ity fraction is not considered since the complemen- tarity equation 4 does not depend on the gap height h. Equation 13 is an algebraic equation which introduces the concept of sensitivity, i.e.the direct inﬂuence of the design parameter hon the objective function fob j: S(~x)d fobj dh=3h2 12rprp+1 2rp~V(1):(14) whereS(~x) is the sensitivity [25, 42]. It is then clear that an iterative approach to the optimization problem can be pursued by updating the geometry by a quantity which is proportional to the sensitivity: hn+1(~x)=hn(~x)+"Sn(~x) (15) 4 where"is an under-relaxation coe cient determined empirically in order to ensure a fast convergence rate without instabilities. It should be noted that a smoothing process is re- quired after every update in order to increase the sta- bility and prevent sharp edges in the solution which could lead to unrealistic shapes. This has been real- ized through a Gaussian ﬁlter as proposed by Bletzinger [43] and St ¨uck and Rung [44] by applying a Gaussian smoothing in the following form: hs(~x)=Fsh(~x)=Z K(~r;)h(~x~r)d~r;8~x2 (16) where hsis the smoothed gap height, his the updated gap height from equation 15, Fsis the ﬁltering operator, ~rrepresents the local ﬁlter radius and is the standard deviation of the Gaussian ﬁlter kernel K. In this work is also referred to as the amplitude of the smoothing. The ﬁlter operator Kis normalized so that: Z K(~r;)d~r=18~x2 : (17) In this way the integral of the shape is not altered, while ﬁne scale oscillations due to the derivative computation are damped [44]. The ﬁlter is applied after updating the geometry in equation 15 and before applying the even- tual constraints. Such an iterative approach is also used, for example, in aeronautics for shape optimizations under the con- straint of the Navier-stokes equation [28]. This ﬁrst op- timization strategy has been implemented in MatlabR through a ﬁnite volume scheme and its detailed are shown in appendix A. The corresponding optimization process is shown in the right branch of the ﬂowchart in ﬁgure 2. The second method to implement the optimization strategy employs a smoothing term directly in the ob- jective function. This way the optimizer will ﬁnd the equilibrium between the maximization of the load car- rying capacity and a certain smoothing hence avoiding unfeasible shapes. This implementation was ﬁrst pre- sented by Van Ostayen for the optimization of dynami- cally loaded bearings in the hydrodynamic regime [32]. The smoothing term is introduced by modifying the ob- jective function in the following way: fobj=W+1 2xZ (r¯h)2d (18) whererepresents the smoothing amplitude, xis the local mesh size and ¯h=hh0is the geometry update as shown in ﬁgure 1. This way, the decision equation Figure 2: Flowchart of the optimization process. The computation of the sensitivity and of the update of the gap height are presented for both the continuous as well as the discrete implementation. 11 includes now the Laplacian of the update of the gap height as additional term: d( ¯p;; p;;h)=r xr¯h 3h2 12rprp +1 2rp~V(1):(19) This partial di erential equation can be solved for ¯h with natural boundary conditions. The implementation of this second strategy follows the path shown in the left branch of the ﬂowchart in ﬁgure 2. The second strategy presents some practical aspects from the computational point of view since the simul- taneous solution of the three partial di erential equa- tions (Eq. 9,10,11) can be e ciently implemented in commercial ﬁnite element software such as ComsolR [33]. This approach is also referred to as the one-shot approach since the solver minimizes the Lagrangian for- mulation by setting its total gradient to zero, as shown by Slawig in [45]. 5 0 0.25 0.5 0.75 100.511.52 x/lh/h 0 analytical Rayleigh step discrete adjoint continuous adjoint0.65 0.751.751.85Figure 3: Validation of the adjoint implementation with analytical and direct methods. Comparison between the analytical solution of the 1D Rayleigh step and the solution of the discrete and continuous adjoint method. The curves are non dimensionalized. The inlet pressure is set equal to the outlet one. Both the ComsolR and the MatlabR implementations also represent the two typical ways of addressing the adjoint method, namely the continuous and the discrete adjoint method respectively. In the continuous adjoint method the three equations 9,10 and 19 are ﬁrst analyt- ically derived and then discretized. In turn, the discrete method exploits the discretized form in order to build the solving system for the adjoint equation 10 and ﬁ- nally solves the decision equation in its algebraic form 14. A detailed derivation of the discrete adjoint method is given in appendix A. A quantitative comparison of the performance of the two approaches has been performed by Nadarajah and Jamson for the adjoint optimization constrained by the Navier-Stokes equation [46]. They show that the two approaches converge to the same solution for a su - ciently ﬁne mesh with minor di erences at the bound- aries where the discrete method shows a better conver- gence despite its slightly higher computational require- ments. The two di erent implementations of the opti- mization scheme for the present work are validated and compared quantitatively against each other in section 3. 3. Validation Due to the novelty of the proposed optimization method, we present a validation of both implementa- tions of the adjoint method in this section. A typical textbook example of an optimization problem for lubri- cant thin ﬁlms is the 1D Rayleigh step which maximisesthe load carry capacity of a 1D slider bearing. Figure 3 shows the comparison between the analytical curve and the two optimal solutions with the continuous and discrete adjoint method. As one can see the analytical shape is fully recovered by both methods, showing an inlet to outlet ratio of 1 :86 and a normalized step posi- tion at 72% of the bearing length [2]. In both optimiza- tion strategies the minimal ﬁlm height was constrained to be equal to 1 :0. If left unconstrained, the optimizer would reduce the ﬁlm height towards zero in order to obtain the maximum load. With respect to the 2D Rayleigh step without cavi- tation, Van Ostayen presented the optimized shape in [32] showing a very good agreement in terms of non dimensional load carrying capacity with the optimized solution of Rhode [32]. Figure 4: Schematic representation of the pin geometry in a pin-on- disc tribometer. The pin has a parabolic proﬁle whose height at the center is hs, while the minimum gap height at the center of the geom- etry is hc. The side of the pin hsideis 1 mm high. The domain is size isLxandLzin the x- and z-direction respectively. Cavitation does not occur in the optimal 1D Rayleigh step solution. In order to check the optimization meth- ods in a cavitated lubricating ﬁlm, we performed a fur- ther comparison with a more complicated geometry, namely a pin-on-disc tribometer. The pin used for this comparison has a circular base and a parabolic proﬁle as shown in ﬁgure 4. The chosen dimensions of the pin match those used in the experimental campaign by Braun et al. [34], where a radius Rpin=4 mm was cho- sen, while the parabolic proﬁle of the pin surface has a shoulder height of 1 m. The discrete method uses a uniform mesh consisting of 1025x1025 grid points cor- responding to a mesh size of 8 :78m. Further informa- tion about the mesh convergence are given in Ref. [14]. There is no analytical solution available for the compar- ison of such a geometry in presence of cavitation, but 6 −5 0 5−4−202·104 x[mm]sensitivityS[N/m]direct method discrete adjoint 0.40.81.2−2.8 −3.3 −3.8·104Figure 5: Validation of the adjoint implementation with analytical and direct methods. Sensitivity in the centerline at the ﬁrst iteration step for the a pin with parabolic shape. The sensitivity is shown with its dimensional value [N /m] since it relates the variation of the normal force (object function) to a variation in the geometry. it is possible to validate the ﬁrst step of the optimiza- tion loop through a comparison with a direct method to compute the sensitivity. In the direct method the total derivative of the Lagrangian formulation in Eq. 8 with respect to the design parameter can be obtained by intro- ducing a Dirac delta perturbation h(~x) in the reference geometry and measuring the e ect of it on the objec- tive function directly. This way, the sensitivity can be computed in the following way: S(~x)=d fobj dh=Wh(~x)Wref h(~x)(20) where Whis the load computed with the perturbed ge- ometry and Wrefis the reference load. Figure 5 shows the sensitivity evaluated in the cen- terline of the pin with both the direct method and the discrete adjoint method at the ﬁrst step of the optimiza- tion loop, when the geometry is still the reference one. As on can see, the two curves are in close agreement. The sensitivity is deﬁned such that a positive value of it means that the gap height should be increased while a negative value implies a reduction of the gap height. Moreover, it is interesting to notice that the second part of the pin surface, where the sensitivity is zero, corre- sponds to the cavitation region. This hints to the fact that there is no advantage to modify the surface in the cavitation region. The simple geometry of the parabolic pin also allows comparing the continuous and discrete adjoint results.hc 1m hs 1m hside 1 mm Lx=Lz 9 mm Rpin 4 mm  0:01 Pas U 0:1 m/s Pamb 101325 Pa Pcav 80000 Pa Table 1: Simulation parameters used in the results section. For this comparison the same uniform mesh is used for discrete method as for the sensitivity computation shown in Figure 5. On the other side, the COMSOL im- plementation of the continuous method employs a trian- gular mesh with a typical element size of 8 m. This is shown in ﬁgures 6a and 6b where the opti- mized geometry and its corresponding pressure proﬁle are shown for both the adjoint methods and the refer- ence case. The continuous and the discrete method are overall in good agreement. The small di erences in the corners of the two optimal solutions are due to the dif- ferent smoothing procedures as explained in section 2. 4. Results and discussion In this section we apply the optimization through the adjoint method to a pin-on-disc tribometer. Similar op- timization attempts documented in the literature treat cases without cavitation [32, 5, 6]. This speciﬁc experi- mental device is chosen since it is representative for typ- ical tribological applications in which cavitation occurs. Moreover, pin-on-disc tribometers are widely used for testing texture designs, hence easing the eventual exper- imental validation of the optimization results presented here. Figure 4 shows the reference geometry of the pin of a pin-on-disc tribometer. The surface has a parabolic pro- ﬁle in order to model the rounding e ects on the edge of the pellet which are typically created by the manufactur- ing processes such as polishing and grinding. Geomet- rical and operating parameters are presented in table 1. Although our COMSOL implementation of the adjoint method (continuos) is computationally faster than the MATLAB one (discrete), the results in this section are computed only by using the discrete adjoint method out of practical reasons in the data management. 7 (a) −4 −2 0 2 4012345 x(mm)h(µm)reference pin discrete adjoint continuous adjoint(b) −4 −2 0 2 412345·105 x(mm)p(Pa)reference pin discrete adjoint continuous adjoint Figure 6: Validation of the discrete and continuous adjoint implementations with the geometry of the parabolic pin in presence of cavitation. a, Comparison between the continuous and the discrete adjoint methods for pin with parabolic shape. The parabola is 1 m high as well as the gap height at the center. b, Comparison of the pressure proﬁle in the centerline of the pin. The comparison is carried out with the following operating parameters for both ﬁgures: U=0:1 m/s,=0:0028 Pas,Pcav=Pin=101325 Pa. 4.1. Sensitivity One of the ﬁrst relevant insights that we can deduce from the adjoint method is the sensitivity map of the surface, or, in other words, the direct inﬂuence of a gap height modiﬁcation on the load carrying capacity as pre- sented in equation 14. Figure 7a shows the sensitivity of the reference pin geometry and ﬁgure 7b shows the corresponding pressure distribution. The sensitivity is computed by evaluating the algebraic sensitivity equa- tion (14) for the initial geometry shown in Fig. 4. As one can see three regions can be clearly identiﬁed. The ﬁrst region covers the front part of the pin and is charac- terized by a positive value of the sensitivity. This means that in this region an increase of the gap height leads to an increase of the normal force. This is, therefore, the part of the pin surface which should be textured in order to achieve a higher load carrying capacity. It is also important to notice that any addition on top of the pin surface, i.e. a reduction of the ﬁlm height, would decrease the load carrying capacity, hence asperities el- ements [47] should not be added here. In the second region the sensitivity is negative, mean- ing that a reduction of the gap height would be bene- ﬁcial for the load carrying capacity. At this point it is up to the designer to choose if the surface can be ma- nipulated in both directions or eventually only in one (i.e.if the gap height hcan be increased or decreased). In the present work we focus only on surface modiﬁ-cations which are based on texturing techniques, hence allowing the surface to be modiﬁed only inwards e.g.by etching or by laser texturing. This is done by setting the additional constraint that the optimized gap height ¯hcan be only bigger than the reference one: ¯h(~x)>href(~x): (21) The last region of the sensitivity concerns the part of the domain where cavitation occurs. In this region the sensitivity is zero, meaning that any geometry modiﬁca- tion would have no consequences on the load carrying capacity (as long as it remains under a certain ampli- tude). This has been already shown in literature, for example by Gercha et al. who studied surface textures in a parallel bearing and showed how the texture had no impact if the whole ﬂuid ﬁlm is cavitated [48]. The same conclusion is reached in the work of Codrignani et al.[14] who investigated the e ects of a single dimple on a pin in a manner similar to the direct method shown in section 3. We note that during the optimization process the size and the position of the three regions may change. It could therefore become beneﬁcial to place texture el- ements in region with initial negative sensitivity after some iterations of the optimizer. 4.2. 2D Rayleigh step Figure 8a shows the pin geometry after convergence of the optimizer, hence until the relative residual of the 8 (a) (b) Figure 7: Analysis of the reference untextured pin. a, Sensitivity over the pin geometry at the ﬁrst optimization step. A positive sensitivity implies that the gap height should be increased in order to achieve a higher load carrying capacity. The vertical yaxis is 1000 magniﬁed in comparison to the other two ones. b, Pressure distribution of the parabolic pin. The cavitation region is delimited by a red line. The lubricant ﬂuid ﬂows from left to right in the x direction. gap height distribution rn h: rn h=Z jhn(~x)hn1(~x)j jhn1(~x)jd (22) drops below 105. The optimized surface consists of a single big pocket in the front part of the pin, the geom- etry of a 2D Rayleigh step. A similar geometry is typ- ically analysed in the literature for slider bearings, as described, among others, by Rohde [5], Shen & Khon- sari [6] and Fillon & Dobrica [18]. Thanks to the optimization the pin almost doubled its load carrying capacity, going from a normal force ofFN=8:64 N to 16:79 N. This is revealed by ﬁgure 9, that shows the convergence of both the normal force and the friction coe cient over the optimization steps. The ﬁnal pressure distribution is shown in ﬁgure 8b, with a maximum pressure of around 1.4 MPa while the un- optimized Rayleigh step yields 0.9 MPa. The adjoint formulation applied in this work assures that the normal force corresponds to the maximum, while, on the other side, it cannot be proven a priori that the friction coef- ﬁcient is globally minimized. Nonetheless, as already stated in section 2, an increase in the load carrying ca- pacity leads to an increase of the gap height which, in turn, corresponds to a friction reduction. Furthermore, optimizing the load carrying capacity and therefore in- creasing the ﬂy height for a given load, will result in an earlier transition from mixed lubrication to hydrody- namic lubrication. The discrete optimizer reached convergence after 1000 iterations which is as computationally expensiveas running 2000 simulations of the Reynolds equation, since each step consists of the solution of the Reynolds equation and the adjoint one. For the same level of ac- curacy, this is faster than describing the Rayleigh step geometry through a parametric study, since we are op- timizing the geometry point by point. In particular, the resolution of the system with a grid of 1025x1025 points took 10 hours on a single core of a laptop. The itera- tion number can be consistently reduced by solving the equations 9,10 and 11 in a single system, as presented in section 2 for the COMSOL implementation of the con- tinuous adjoint method. As a next step we would like to consider more com- plex constraints and their impact on the optimal solu- tion. In particular the geometry of the 2D Rayleigh step could introduce some experimental challenges if tested on a pin-on-disc tribometer because of its behaviour in the mixed lubrication region where contact could occur on the sharp edges of the Rayleigh step. In order to avoid this e ect we add a further constraint which will force the initial contact to occur in the central part of the pin. This can be done, for example, as shown in ﬁg- ure 10 where the a central stripe along the centerline of the pin is constrained to remain untextured. The intro- duction of such a constraint reduces the improvement in the load carrying capacity. In particular the resulting normal force is FN=14:42 N which is just two thirds of the performance obtained with the full 2D Rayleigh step. As one can see from ﬁgure 10 the sensitivity on the front part of the stripe is positive, meaning that the constraint prevents the increase of the gap height in a 9 (a) (b) Figure 8: Optimized pin geometry and pressure distribution with the 2D Rayleigh step. The optimization has the only constraint that the gap height cannot be reduced as stated in equation 21. a, Geometry of the optimized 2D Rayleigh step. The sensitivity is plotted over the surface of the pin. The vertical yaxis is 1000 magniﬁed in comparison to the other two ones. b, Pressure distribution over the optimitzed 2D Rayleigh step. The cavitation region is delimited by a red line. The lubricant ﬂuid ﬂows from left to right in the xdirection 1 250 500 750 1,000024·10−3 friction coeﬃcientplot 1 1 250 500 750 1,00005101520 iterationsnormal force FN(N) normal force friction coeﬃcient Figure 9: Convergence of the normal force and of the friction coe - cient during the optimization iterations. The simulations are stopped as soon as the residual on the gap height variation ¯hdrops below 105. Figure 10: Sensitivity over the optimized pin geometry with the ad- ditional constraint which eases the e ective texturation on a real pin- on-disc tribometer. The central stripe has a width of 30% of the pin radius.The vertical yaxis is 1000 magniﬁed in comparison to the other two ones. 10 region where this would be beneﬁcial for the load car- rying capacity. 4.3. Optimization of surface textures In order to apply the adjoint method to the optimiza- tion of surface textures we need to add two further con- straints: the texture pattern and the size of each element. Among the numerous shapes that texture elements can have we focus here only on elements which have a cir- cular base, where each grid point inside the texture ele- ment remains an optimization variable. The actual local depth of the optimal texture element will be automati- cally determined by the algorithm. The texture pattern is a pseudo-hexagonal one in which each texture ele- ment is equally spaced in both directions and each sec- ond row is shifted of half the distance in the z-direction as already used in [14]. Figure 11a and 11b show the optimized textured sur- face and the corresponding pressure distribution for a pseudo-hexagonal pattern made by circular dimples with a diameter D=200m. As we can see the tex- ture elements in the part of the domain with positive sensitivity were automatically selected to increase in depth. This is one of the biggest advantages of the ad- joint method, because usual parametric techniques place texture elements also in region where their presence is detrimental [4]. Moreover, the gap height distribution inside the texture elements is optimized point by point, meaning that also non circular shapes may be found. This can be seen in the elements on the lateral side of the pin which assume a rather elongated shape. This shape helps the lubricant ﬂow to reach the central part of the pin in order to enhance the load carrying capac- ity by bringing more ﬂuid where the gap height is lower. This could also explain why some elements between the fully circular ones and the elongated ones have a barrier inside them which splits the elements in two. All texture elements in this case have rather sharp edges because of the low amplitude of the smoothing =0:3 applied at each iteration step. By increasing the smoothing amplitude to =1:2 the shapes of the texture elements become more uniform as one can see in ﬁgure 12a. Nonetheless the elements on the side still assume an elongated form. The smoothing has hence a non negligible impact on the optimal solu- tion and thus on the value of the load carrying capacity. In fact the normal force of the slightly smoothed case isFN=9:1823 N while the heavily smoothed case de- livers a lower force FN=9:0623 N. Even though both cases yield an increase of the load carrying capacity, it is clear that an excessive smoothing could harm the op-timal performance, as reported also for other works on the adjoint optimization [44]. The di erence between the two values of the smooth- ing amplitude can be seen in ﬁgure 12b where the proﬁle at the centerline of the pin is shown for both cases. The less smoothed elements exhibits sharper edges which resemble those of a cylindrical dimple while the smoothed and less e ective dimples have a shape similar to a spherical one. This aspect should be considered in experimental activities because laser sur- face texturing is unable to create very sharp edges. This ﬁnding is in agreement with the results of other works which investigated the texture shape though parametric studies, such as Aggarwal & Pandey [49], Adjemout et al.[35] and Nanbu et al. [50]. The spacing lxandlzbetween texture elements is cho- sen so that the texture density is txt=15% if the whole pin surfaces is textured. In practice, since dimples ap- pear only whereS>0, the overall texture density is lower. This hints at the fact that the deﬁnition of texture densitytxtshould take into account the sensitivity anal- ysis. In fact the optimal value of texture density is sim- ply 100% of the region with positive sensitivity which corresponds to the 2D Rayleigh step shown in ﬁgure 8a. As ﬁnal comment to this subsection we note that spe- cial care has to be taken in the analysis of the validity of optimized shape, since texture elements with an exces- sive depth-to-diameter ratio may violate the applicabil- ity assumption of the Reynolds equations as discussed by Dobrica and Fillon [51]. The optimized texture in this section has a maximal depth-to-diameter ratio of 1=200 which is way below the limit of 1 =10 found in [51]. 4.4. Comparison between di erent optimization tech- niques for surface textures The great advantage of the adjoint method in com- parison to typical parametric studies is its ability to op- timize the whole surface point-by-point. In particular, in the case of surface textures, the shape of every single texture element can be optimized according to its posi- tion on the macro-geometry. In order to assess the e ec- tiveness of this optimization technique a comparison is carried out between the optimized surface presented in the previous section and the results of the same kind of texture pattern which is optimized through a parametric study with spherical dimples. Figure 13 shows a comparison of the load carrying capacity with the untextured reference case of the di er- ent constraints seen in the previous section and also with three parametric studies. At ﬁrst we focus on the impact 11 (a) (b) Figure 11: Optimization of texture on the pin surface. The texture has a pseudo-hexagonal pattern and the dimples are constrained to have a maximal diameter D=200m. a, Sensitivity over the textured pin surface with low smoothing amplitude =0:3. The vertical yaxis is 1000 magniﬁed in comparison to the other two ones. b, Pressure distribution. The cavitation region is delimited by a red line. (a) (b) −4 −3 −2 −1 0012 x(mm)h(µm) σ= 0.3 σ= 1.2 Figure 12: Analysis of the impact of the smoothing amplitude . a, Sensitivity over the pin geometry. With high smoothing amplitude =1:2. The vertical yaxis is 1000 magniﬁed in comparison to the other two ones. b, Comparison of the centerline gap height for di erent values of smoothing amplitude. The gap height is 1000 magniﬁed in comparison to length of the pin. 12 Figure 13: Performance of the adjoint optimization for di erent constraints. The ﬁgure shows also the comparison with parametric studies for textured surfaces done in this work. The normal force is displayed as percentage of the untextured case a. The operating parameters of this comparison are reported in Table 1. of the constraints on the load carrying capacity by com- paring the ﬁrst ﬁve columns. As one can see the best increase of the normal force is obtained by the Rayleigh step (column b), hence by optimizing each point of the surface under the only constraint that the surface can- not be extruded (equation 21). Column c concerns the case with an additional stripe in the middle of the pin (see ﬁgure 10) which acts as a further constraint which reduces the e ectiveness of the Rayleigh step. Columns d and e represent the adjoint optimization of the surface texture presented in section 4.3 which has texture ele- ments arranged in a pseudo hexagonal pattern. As we can see the introduction of such an invasive constraint has lowered the increase of the normal force, showing that the Rayleigh step is the best way to increase the nor- mal force in conformal contacts. Nonetheless columns d and e show a higher normal force than all the remain- ing columns which correspond to the parametric texture optimization. In particular we analysed in this work two recurrent cases in literature, namely the presence of tex- ture elements on the whole surface (last column) or only on the front part (second last column). Front textures areknown to perform better than total ones and the reason can be clearly seen from the sensitivity distribution in ﬁgure 7a where the sensitivity in the rear part of the pin is either negative or zero. Therefore, we also consider the case in which the parametric optimization is carried out by placing texture elements only in regions of pos- itive sensitivity, see ﬁgure 7a. The normal force of this case is shown in the third last column where we can ob- serve that thanks to the additional information provided by the computation of the sensitivity, this case is the best among the parametrically optimized ones. However, the texture optimized with the adjoint method performs bet- ter since it can ﬁnd the optimal shape and depth of ev- ery texture element according to their position on the surface. It is important to notice that, for the chosen operating parameters and geometry, the absolute improvement in load carrying capacity is small ( <7%) with every opti- mization techniques shown in ﬁgure 13. This could be due to the fact that the simulations are carried out in the hydrodynamic regime. An higher impact of the texture may be encountered in the mixed lubrication regime, as 13 foundboth experimentally and numerically in previous works [34, 52, 53]. As ﬁnal remark for this section we would like to un- derline that, although the practical realization of tex- tures with di erent element size may result challeng- ing [54, 55], a great improvement can be still achieved through the sensitivity analysis and assuming a constant shape for each texture element as shown in ﬁgure 13 column f. 4.5. Generalization of the optimized geometry for dif- ferent operating conditions As seen in the previous subsection, the adjoint method ﬁnds optimal texture elements which have dif- ferent shape and depth according to their position on the macro geometry. Nonetheless some characteristics of the optimized texture can be generalized as function of the operating conditions so that some general design rules can be deﬁned. This can be done for example by analysing the average depth of the texture elements as a function of the gap height for di erent combinations of viscosity and velocity as shown in ﬁgure 14. The dy- namic viscosity is varied in the range 0 :0028<< 0:28 Pas while the wall velocity ranges between 0 :01 and 1 m/s. The simulations are carried at imposed gap height from 0:5 to 10m. This results in a variation of the Sommerfeld number based on the gap height between 0:092106and 8:49106. As one can see the optimal depth varies linearly with the gap height. This linear behaviour has been found also by parametric studies in literature [14, 8, 53]. Nonetheless, the particularity of the adjoint method is that the texture elements appear only where the sensitivity is positive, hence avoiding to have texture elements which actually have a detrimental e ects on the load carrying capacity. 5. Conclusion We present for the ﬁrst time the implementation of the adjoint method for the optimization of surface tex- tures in hydrodynamically lubricated contacts with cav- itation and its validation. This method allows to ﬁne- tune each single texture element individually rather than through a trial and error approach with dimples of iden- tical geometry. The adjoint method is more e ective than traditional parametric approaches for the optimiza- tion of the texture shape. This is because it can de- termine where textures should be and how the texture shape should look like in a single simulation without the need to parametrize the texture with elementary shapes. This allows to carry out free-form optimizations in a computationally e cient manner. Figure 14: Average dimple depth as function of the gap height. Each point represent a di erent combinations of velocity and viscosity so that the following range of Sommerfeld number is swept: 0 :092106< S<8:49106. The simulations are carried out at constant gap height. We ﬁnd that the traditional approach to texture the full front of the pin does not lead to an ideal increase in load carrying capacity. The sensitivity analysis reveals that it is beneﬁcial to texture only part of the front por- tion of the pin. We carried out these optimizations as a function of operating conditions, conﬁrming a linear relationship between the gap height and optimal texture depth. Further development of this work should include also thermal e ects and the elasticity of the contact sur- faces [56, 57, 58, 59]. Acknowledgements The authors acknowledge support from the Euro- pean Research Council (ERC-StG-757343). The Au- thors wish to thank Dr. Johannes Schneider, Dr. Paul Schreiber, Philipp Mayer and Oscar Furst for their sug- gestions in the preparation of this work. 6. Appendix A, derivation of the discrete adjoint method In this appendix we show the detailed derivation of the discrete adjoint method for the Reynolds equation with mass-conserving cavitation. The Reynolds equa- tion 1 can be written in its linear formulation with the discretized unknown pand : g=Ap+Bc=0 (23) where gis a vector containing the residuum of the Reynolds equation. The discrtization is performed with 14 the ﬁnite volume method [60]. The complementarity constraint in equation 4 is discretized by deﬁning its residuum fin the following way: f=¯p+q ¯p2+2=0: (24) The nonlinear system can be solved with the Newton- Raphson method and its corresponding Jacobian matrix can be written in the following way [36]: J p ! ="Jf;pJf; Jg;pJg;# ¯p ! = f g! : (25) The computation of the Jacobian matrix is eased by the fact that its submatrices are already known from the dis- cretization of the Reynolds equation ( e.g.Jg;p=A). After having discretized the direct problem we can rewrite the Lagrangian formulation in its discrete for- mulation: L=fobj+T pg+T f (26) where pand are vectors of Lagrangian multipliers. Since gandfare zero everywhere by deﬁnition, one can choose pand freely so thatL=fobj. This allows us to rewrite the sensitivity as Sd fobj dh=dL dh= @fobj @¯p+T p@g @¯p+T @f @¯p!d¯p dh + @fobj @+T p@g @+T @f @!d dh +T p@g @h+T @f @h: (27) By setting the two terms in the brackets on the right hand side to zero we obtain the system of the adjoint equations:266664@g @¯p@f @¯p @g @@f @377775 p ! =0BBBB@@fobj @¯p 01CCCCA (28) where the term@fobj @is set to zero since the objective function does not depend on the gap height. The re- maining right hand side term@fobj @¯pcan be easily com- puted from the analytical deﬁnition of normal force in equation 5: @fobj @p=@W @p=@ @pZ (ppext)d =Z 1d :(29) Therefore@fobj @¯pis a 1D array which contains the area of the each single cell d . Equation 28 is a linear algebraic equation whose solution can be obtained per substitu- tion after having preconditioned the system in a similarmanner as presented by Wolozynski et al. [36] for the cavitation algorithm. The next step is to compute the sensitivity by solving the algebraic equation 27 which now reduces to: S=T p@g @h(30) since the term T @f @his zero because the complementar- ity constraint does not depend on the gap height. In or- der to compute the sensitivity we ﬁrst compute the ana- lytical derivative of the Reynolds equation with respect to the gap height, i.d.the term@g @h: @g @h=@A @h¯p+@B @h (31) where@A @hand@B @hcan be easily obtained from the imple- mentation of the Reynolds equation 23. The optimization can be now performed in an iter- ative way as explained in the ﬂowchart in Figure 2 in Section 2. 7.