Quantifying the effect of asperity size and shape on thermal contact conductance of metal-metal contacts through direct numerical simulations

International Journal of Heat and Mass Transfer 115 (2017) 336–346

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Quantifying the effect of asperity size and shape on thermal contact conductance of metal-metal contacts through direct numerical simulations Navni N. Verma, Sandip Mazumder ⇑ Department of Mechanical and Aerospace Engineering, The Ohio State University, Columbus, OH 43210, USA

a r t i c l e

i n f o

Article history: Received 21 April 2017 Received in revised form 24 July 2017 Accepted 11 August 2017

Keywords: Thermal contact conductance Interface TCC Scale-resolved Numerical Asperity shape Asperity size

a b s t r a c t The thermal contact conductance (TCC) between two metallic rough surfaces is dictated both by the size and shape of the asperities on the surface. In this study, scale-resolved direct numerical simulation of thermal transport across the interface was conducted to isolate and quantify their effects. The mean roughness height (indicator of asperity size) and the mean absolute slope (indicator of asperity shape) of the asperities—two common surface topography descriptors—were treated as parameters, in addition to the contact pressure. Microscale models of the interface geometry were generated by statistically reconstructing the surface pairs from the aforementioned topography inputs, followed by generation of unstructured meshes that resolved all fine-scale features of the interface, including the air pockets. Steady state conjugate heat conduction computations were conducted, and the TCC values were extracted. The extracted TCC values were found to be in good agreement with experimental measurements. Examination of trends revealed the TCC to be a strong function of asperity size and a weak function of asperity shape. The extracted TCC data was finally reduced and curve-fitted by correlating it to two non-dimensional parameters: the number of contacts and the mean absolute slope. The accuracy of this correlation (curve-fit) was tested, and the TCC values predicted by it were found to be well within the statistical uncertainty of the raw DNS results. The curve-fit represents a simple, yet valuable tool for future researchers to estimate the TCC of metal-metal interfaces of any topography characterized by common topography descriptors. Ó 2017 Published by Elsevier Ltd.

1. Introduction Engineered surfaces, when viewed under the microscope, exhibit inherent surface roughness. When two surfaces are placed in contact, tiny air pockets are trapped within the interface and the actual contact area is a small fraction of the nominal (or apparent) contact area, even at high contact pressure [1]. The air pockets trapped at the interface hinder heat transfer by conduction between the surfaces since the thermal conductivity of air is three orders of magnitude smaller than that of most metals. On account of this thermal conductivity mismatch, heat is constricted to flow through regions where the two surfaces actually come in contact. This resistance to heat flow gives rise to a sharp temperature drop across the interface. The effective thermal conductance of an inter⇑ Corresponding author at: Department of Mechanical and Aerospace Engineering, The Ohio State University, Suite E410, Scott Laboratory, 201 West 19th Avenue, Columbus, OH 43210, USA. E-mail address: [email protected] (S. Mazumder). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.08.031 0017-9310/Ó 2017 Published by Elsevier Ltd.

face, which is the inverse of this resistance, is often referred to as the thermal contact conductance (TCC). Metal-metal contacts, which are the central focus of this study, are omnipresent in engineering applications. Hence, accurate quantification of the TCC of metallic interfaces is critical to the effective analysis and modeling of heat transfer in machine components and their subsequent design. Over the past several decades, a number of simplified models for predicting the TCC of metallic interfaces have been developed. In the vast majority of these models, the TCC has been computed as the net sum of the conductance of the discrete contact spots. Greenwood [2] proposed a method to compute the TCC as the sum of the parallel conductance of all the contact spots with an assumed surface height distribution at the interface. Cooper et al. [3] proposed a similar methodology, in which the single circular contact hypothesis was extended to multiple contacts with a distribution of asperity heights derived from surface topography measurements. Mikic [4] proposed a simple correlation of the TCC of conforming rough contacts to the root-mean-square roughness

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

337

Nomenclature ai A d f h k Nc P Q Ra Ry y T DT

constants in curve-fit [Eq. (5)] apparent area of contact (m2) separation distance (m) cumulative distribution function thermal contact conductance or TCC (W/m2/K) thermal conductivity (W/m/K) number of contacts per one thousand asperity pairs probability heat transfer rate (W) arithmetic mean roughness (m) random number surface profile height (m) temperature (K) temperature difference (K)

and average of the absolute value of the asperity slope. The aforementioned works assume that all the contact spots at the intersection of perfectly aligned asperities are circular with a constant radius. In reality, contact spots of varying shapes and sizes reside at the interface. Furthermore, detailed distribution of contacts in real interfaces is not commonly measured. Only surface topography descriptors, such as average surface roughness heights and their root-mean-square variance (or standard deviation), are known. Black et al. [5] extended the aforementioned models to include the effect of contacting spheres of varying radii. Hong et al. [6] proposed an integrated thermo-mechanical model, which accounted for partial contact between asperities. Singhal et al. [7] proposed a coupled thermo-mechanical model, in which actual surface profile data, in conjunction with deformation analysis, was used to predict the actual contact area responsible for heat transfer. They approximated the contact between two rough surfaces with the contact between a single rough surface with equivalent characteristics and a perfectly smooth surface. Fundamentally, the TCC depends on the actual area of contact when two surfaces are placed against each other. The contact spots are usually discrete spots that result when the tallest asperities of the two surfaces touch and deform under the influence of an external load (or pressure). The actual contact area depends on the load applied across the interface and also the structure of the asperities on the two surface pairs. As the load is increased, new asperities may come in contact, which, in turn, may increase the contact area and affect the heat flow across the interface. Thus, it is fair to hypothesize that the contact area primarily depends on the number of contacts, which in turn, depends on the applied load. This hypothesis has been shown to hold water in a previous study [8], which shows that the number of contacts is the dominant factor determining heat transfer across the interface. If the number of contacts was the only factor determining the TCC, then, it would imply that when surface pairs with the same arithmetic mean roughness heights (Ra ) and standard deviation (ra ) are brought in contact, the TCC would remain the same. However, such an implication does not agree with reality. Experimental measurements show that the TCC depends not only on the height distributions of the asperities, but also on the shape of the asperities [9–12]. This is to be expected since the shape of the asperities that come in contact would definitely have a bearing on the contact area through which heat transfer occurs. In order to prove that the height distribution of the asperities (which translates to number of contacts) is not the only factor determining the TCC, following a previous study [8], surface pairs with the same average and rootmean-square roughness heights were statistically reconstructed, and the resulting TCC computed. For this particular study, a mini-

Greek symbols standard deviation of surface heights (m) profile heights measured from a reference plane (m) d thickness (m) Da mean absolute slope

ra l

Abbreviations 2D two-dimensional 3D three-dimensional DNS direct numerical simulation TCC thermal contact conductance

mum of 10 ensembles were used. In other words, 10 different surface pairs were reconstructed statistically from the same average roughness height and root-mean-square values. Fig. 1(a) shows two such ensembles. In both ensembles, the number of contacts is the same per unit area (3 in the section shown in the figure). For each ensemble, scale-resolved direct numerical simulation of heat conduction across the interface was conducted, and the TCC was computed following the procedure outlined in a previous study [8]. From the 10 ensembles, an average TCC and its standard deviation were also computed. The number of ensembles was then progressively increased from 10 to 35, and the average and standard deviations were recomputed. Fig. 1(b) shows the standard deviation of the TCC—normalized by the average and expressed as a percentage—as a function of the number of ensembles. If the height distribution of the asperities (average and root-mensquare heights) was the only parameter governing the TCC, then, the standard deviation of the TCC is expected to asymptote to a value of zero as the number of ensembles is increased since each ensemble is generated using the same two inputs. The fact that Fig. 1(b) shows that the standard deviation of the computed TCC does not decrease beyond a certain value as the number of ensembles is increased, indicates that even though the asperities are generated from the same roughness height distribution, since they have fundamentally different shapes in each case, the resulting TCC values are not perfectly correlated to the roughness height distribution alone. In order to account for the effect of asperity shape on the TCC, a parameter to characterize asperity shape is needed. A surface texture parameter commonly employed to measure the horizontal variation of a profile is the mean absolute slope, Da , of the surface profile. While there are several parameters that can be used to describe the shape of asperities [13], there are some compelling reasons that make the mean absolute slope the most favorable choice for the purpose of this study. First, the mean absolute slope is one of the few topography parameters that can be directly processed from the measured profile data and is available as output from most surface profile measuring equipment. Second, the vast majority of experimental investigations conducted to measure the TCC and analyze factors that influence the TCC, typically characterize surface texture using the mean absolute slope [14–21]. Fig. 2 illustrates two surfaces statistically reconstructed from the same height distribution, but with different mean absolute slopes. It is clear that the surface with the higher mean absolute slope value has more slender asperities. Since the height and width (or base area) of an asperity together define its slope, there are two ways of varying the mean absolute slope of surfaces. The first is to vary the width of the asperities,

338

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

Fig. 1. (a) Two ensembles statistically reconstructed from the same roughness height distribution and (b) variation of computed TCC as a function of the number of statistical ensembles.

Fig. 2. Illustration of two surfaces statistically reconstructed from the same height distribution, but with different mean absolute slopes. Both surfaces are sampled from a Gaussian height distribution and have an arithmetic mean roughness height, Ra , of 5.22 lm. The surface on the left has a mean absolute slope, Da , of 4.03, while the one on the right has a mean absolute slope of 1.64.

while keeping the heights statistically the same. Zhang et al. [9] measured TCC values for two surfaces with the same arithmetic mean roughness height (or Ra ) values but different mean absolute slope (or Da ) values, which resulted from the different shape and material of the grains of the two surfaces. They found that the TCC increases linearly with the mean absolute slope of surfaces subjected to the same load and arithmetic mean surface roughness height. TCC predictions of their numerical model also showed good agreement with experimental results. The second way to change the mean absolute slope is to change the asperity heights by changing the arithmetic mean roughness height, Ra . The effect of dissimilar Ra values on the TCC has been investigated by several groups [10–12,14]. All studies reported that at a given applied pressure, the TCC is higher at lower values of Ra . Based on these studies, it appears that the method used to vary the mean absolute slope (or asperity shape) has contrasting effects on contact conductance. This is probably because the slope is inherently connected to

the height of the asperity, and its effect cannot be isolated in experimental measurements. The objective of the present study is to isolate and quantify the effects of asperity height and asperity shape through spatial scaleresolved DNS. DNS involves several key steps. First, rather than make assumptions with regard to the shape and height of the asperities, the topography of the interface is first statistically reconstructed from commonly measured surface roughness descriptors, namely the mean roughness height (Ra ), the rootmean-square height (ra ), and the mean absolute slope (Da ). Next, the interface geometry and all relevant length scales—both in the metal and air—are spatially resolved using an unstructured mesh. Finally, the steady state heat conduction equation is solved numerically on this unstructured mesh, and the temperature and heat fluxes obtained from such calculations are ultimately postprocessed to extract the TCC. In order to make the extracted TCC values accessible and worthwhile to future researchers, a reduced

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

model is also pursued, wherein a simple correlation between the TCC and the aforementioned two parameters is proposed. 2. Theory and research method 2.1. Statistical reconstruction of surfaces Over the past several decades, many different models have been proposed for describing the surface roughness. Broadly, these models may be classified as: (1) deterministic models, in which the surface profile data is extracted directly from stylus-based roughness measuring instruments [22,23] and (2) statistical models, which use appropriate asperity height distribution functions to describe the surface topography [24,25]. Fractal models, which account for the multiscale nature of surface roughness, have also found limited use [26–28]. Although deterministic models are preferred since they represent the topography exactly, these models require finescale topography measurements in order to generate surface images. Such data can only be collected using high-resolution profilometry, or digital microscopy, followed by processing of the measured data to fit the deterministic model. Unfortunately, such high-resolution data is rarely collected in the industry, and almost never documented. Manufacturing methods of most engineered surfaces are standardized across the industry. Generally, the only information available for describing the surface topography of machine components are a couple of numerical values, such as the arithmetic mean roughness, Ra , and/or the root-mean-square roughness height, ra . If the asperity heights follow a Gaussian distribution, as is often the case [29,30], the arithmetic mean roughness is related to the root-mean-square roughness height through the simple relation: ra ¼ 1:25  Ra . In the vast majority of cases, Ra is the only parameter available to describe a surface, and consequently, a statistical model is necessary for surface reconstruction. In this study, to preserve compliance with the data available from the industry, a statistical approach that geometrically reconstructs surface topographies using Ra values and appropriate height distribution functions was developed. A schematic of an interface between two surfaces and the associated surface topography parameters used in the geometric reconstruction is shown in Fig. 3. The arithmetic mean roughness, Ra , characterizes the vertical deviation of the surface heights relative to a mean height, and is mathematically written by [31]

Ra ¼

1 A

Z 0

A

jy  ljdA

ð1Þ

where A is the sampling area (length in a two-dimensional geometry) of the profile, and y is the profile height sampled from a reference plane, as illustrated in Fig. 3. The mean, l, is computed from

339

the average of the profile heights measured from the reference plane. Other than the arithmetic mean roughness and the root-meansquare height (which is related to the arithmetic mean roughness for Gaussian surfaces), another property that may vary from surface to surface is the mean absolute slope. As discussed earlier, the mean absolute slope may be altered either by changing the arithmetic mean roughness itself or the width of the asperity bases. The latter is preferable for modeling purposes since it can alter the mean absolute slope independently. In this study, the mean absolute slope was computed by fitting a spline to the sampled roughness heights, and then computing the derivatives by differentiating the spline equations directly. The surface pairs reconstructed in this study are the same ones considered by Litke [32] in their experimental study. In their study, only the Ra values of the surface specimens used in the experiments and the stylus diameter of the profilometer were specified. The specimens were prepared using bead-blasting. The profilometer stylus diameter provides an estimate of how often the heights are sampled. For example, no peak or valley that is within a distance of 4 lm can be resolved by a profilometer stylus that has a diameter of 4 lm. This is because, for contact-type roughness measurement instruments, the surface profile is generated from the locus of points detected by a stylus. As a result, the spatial resolution of the roughness features is limited by the diameter of the stylus. Here, a nominal sampling distance of 4 lm was first used. Subsequently, smaller sampling distances were also considered. For the same Ra value, the smaller the sampling distance, the more slender the asperities, and the larger the mean absolute slope Da , as shown in Fig. 2. In other words, the sampling distance was used as an indirect means of independently varying the mean absolute slope, which, as discussed earlier, is the parameter of choice for characterizing the shape of the asperity. The number of asperities that actually come in contact when two surfaces are brought together is strongly dependent on the Ra values of the two surfaces and the separation distance under consideration. The separation distance, d, is the distance between the mean planes of the two surfaces in contact, as illustrated in Fig. 3. It was also treated as a parameter in this study. The separation distance is a direct manifestation of the pressure (load) applied across a contacting interface. The prediction of the separation distance as a function of the load is a complex contact mechanics problem that requires high-fidelity nonlinear finiteelement analysis including stiction, material flow, and surface deformation, and is well beyond the scope of this particular study. For the purpose of this study, the separation distance values for the reconstruction of the interface were extracted from the results of a simplistic surface deformation model developed by Singhal [33]. Table 1 shows the separation distance values extracted from Singhal [33] corresponding to the interface pressure range considered

Fig. 3. Illustration of the interface between two surfaces showing relevant surface topography descriptors.

340

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

Table 1 Theoretically estimated number of contacts for two surfaces with Ra values of 5.02 and 5.22 mm as a function of the pressure applied across the interface. Pressure (MPa)

Separation distance d (mm)

Probability of contact (%)

Number of contacts with 104 samples

5.282 4.843 4.401 3.968 3.498 3.080 2.649 2.199 1.763 1.330 0.892 0.461

25.45 25.78 26.15 26.55 26.92 27.51 28.16 28.95 29.43 30.94 32.41 34.70

0.247 0.220 0.193 0.168 0.147 0.119 0.093 0.069 0.058 0.032 0.017 0.006

25 22 19 17 15 12 9 7 6 3 2 1

in the experimental analysis [32]. For the validation study considered here (Section 3.1), the data from Table 1 are used directly with the understanding that this may be a source of error. For the data shown in Table 1, the surface specimens used have Ra values of 5.02 (top) and 5.22 mm (bottom). When two rough surfaces are brought within a prescribed separation distance, d, there may or may not be any contact. Whether they come in contact or not will depend on the value of d and the Ra values of the two surfaces. If d is significantly smaller than the sum of the Ra values of the two surfaces, some contact is expected. As to how many contacts actually occur can only be determined by statistically reconstructing the two surfaces and searching for overlaps. Clearly, this requires well-conceived computational algorithms that will be elaborated upon shortly. For a back-of-theenvelope estimate, however, a simple theoretical approach may be adopted. For a given d, contact between two asperities will occur when the sum of the heights of a pair of asperities is equal to or greater than d. This means that the probability of contact between two surfaces may be expressed as

Pðy1 þ y2 P dÞ ¼ 1  Pðy1 þ y2 < dÞ ¼ 1  f ðdÞ

ð2Þ

where y1 and y2 are asperity heights following a Gaussian distribution function, with means and variances ðl1 ; r2a;1 Þ and ðl2 ; r2a;2 Þ respectively. By definition, the distribution of the sum, ðy1 þ y2 Þ, will also be Gaussian, with mean ðl1 þ l2 Þ and variance ðr2a;1 þ r2a;2 Þ. As a result, Pðy1 þ y2 < dÞ, which is the probability that the sum of two surface asperities is less than the separation distance, can be computed from the cumulative distribution function, f , of the sum. Such back-of-the envelop calculations are critical since they provide a rough estimate of the number of asperity pairs that must be reconstructed in an actual numerical calculation in order to get statistically meaningful results. For example, Eq. (2) was used to compute the probability of contact for various prescribed separation distances. The probability was then multiplied by the number of samples (= 104) and rounded to integer values to estimate the number of contacts, which is shown in Table 1. It is clear that for large d, even 104 samples (or asperity pairs) are not sufficient to produce a statistically meaningful number of contacts. For example, at a separation distance of 30.94 mm, only 3 asperity pairs come in contact out of the 104 pairs sampled with a standard deviation of 1 (not shown in the table). If 104 asperity pairs were to be reconstructed geometrically, the resulting mesh would be too large for scale-resolved direct numerical simulations to be tractable. To avoid this problem, in this study, computations were restricted to 103 asperity pairs. However, as just discussed, since 103 asperity pairs are not sufficient for statistically meaningful results, such calculations were only conducted for smaller separa-

tion distances, and an indirect approach, as described shortly, was adopted for larger separation distances. The topography of real surfaces is complex. One of the tasks of this study was to develop an efficient algorithm to reconstruct surface topography resembling real surfaces using the aforementioned profile descriptors. The commercial mesh generation software CFD-GEOMTM was used in the present study for reconstruction of the surfaces. This particular software has geometry creation tools that allow interfacing with Python scripts. This feature enables parameter-based geometry creation and subsequent parametric variation of the geometry. A Python script was written to enable automatic reconstruction of surfaces from profile descriptors. The asperity heights were sampled from a Gaussian deviate. This was done by discretizing the height distribution ranging from 6ra to þ6ra into 10,000 intervals, and then computing the cumulative distribution function, f , for each interval, i, using

f i ¼ 1 þ erf

  yi  l pffiffiffi ra 2

ð3Þ

where yi is the height of an asperity in the i-th interval, and l and ra are the mean and standard deviation, respectively, of the surface height distribution. Since the cumulative distribution function, f represents the fraction of surface heights lying at or below a specific height as a function of that height, for each asperity, a random number Ry was drawn and the criterion f i < Ry < f iþ1 was checked over the discretized height distribution. If the above criterion was satisfied, a point at a height corresponding to the height of the i-th interval, yi , was created with reference to a mean line. Following the generation of the heights, the surface was constructed by fitting a spline through all the points. As discussed earlier, the mean absolute slope was directly computed from the spline. Upon completion of the reconstruction of the two surfaces separated by a distance, d, the Python script detects the overlapping asperity pairs and deletes the overlapping regions, thereby creating the contact plane (line in 2 D). Both 2 D and 3 D interface models were created. Fig. 4(a) illustrates the reconstructed 2 D interface (only a small portion is shown) resulting from two surfaces with Ra values of 10 mm and zero separation distance. The Python script for the 2 D model not only reconstructs the surface topographies but also detects the overlapping asperities and creates the contact spots, thereby automating the complete process of geometry creation. A 3 D interface model with the same parameters is illustrated in Fig. 4(b). Unlike for 2 D interface models, in the case of 3 D models, the contact of the intersecting asperities could not be automated completely with the Python script, and some steps of the reconstruction process were executed manually, making 3 D surface reconstruction extremely tedious and time-consuming. For both 2 D and 3 D, the interface geometry created using the Python-script based CFDGEOMTM software were imported into the mesh generation software ICEM-CFDTM, which was then used to generate unstructured meshes. This choice was prompted by the fact that the mesh quality of the meshes generated by ICEM-CFDTM was found to be superior to those generated by CFD-GEOMTM. Fig. 5(a) shows the unstructured mesh generated for the 2 D interface model shown in Fig. 4(a). It is comprised of 400,288 triangular elements. The tetrahedral mesh generated for the 3 D interface model shown in Fig. 4(b) is illustrated in Fig. 5(b). The 3 D volumetric mesh is comprised of a total of 7,491,069 tetrahedral elements. It is clear from Fig. 5 that the mesh is fine enough to resolve all the relevant length scales. Hence, the method is referred to as scale-resolved DNS. 2.2. Calculation of thermal contact conductance (TCC) Following mesh generation, the steady state heat conduction equation was solved. This was performed using the finite-volume

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

341

Fig. 4. Reconstructed interface with Ra ¼ 10 mm and d ¼ 0: (a) 2D and (b) 3D.

based commercial code ANSYS-FluentTM. A finite-volume based code was chosen over a finite-element based code because the finite-volume method guarantees energy conservation irrespective of mesh size [34]. A small temperature difference of 10 K was applied between the top and bottom leads. Small temperature difference avoids potential nonlinearities associated with temperature dependent thermal conductivity. Since the length scales associated with the air pockets are extremely small, the resulting Grashof numbers are miniscule, and natural convection within the air pockets may be safely neglected. Radiation fluxes were also estimated, and found to be less than 1% of the total heat flux in the worst-case scenario. Hence, only conduction through the air pockets was ultimately considered. Fig. 6 shows the temperature distribution for a 2D model of contact between two aluminum surfaces with a separation distance, d ¼ 15 lm between their mean planes. The top surface has an Ra value of 5.02 mm and is held at 308 K, while the bottom surface has an Ra value of 5.22 mm and is held constant at 298 K. In the top and bottom homogenous metal blocks, the heat flows uniformly. Across the interface, however, the heat flow is impeded by the air gaps and the heat gets redirected through the contact spots. Thermal conductivity values of 190 W/m/K and 0.0263 W/m/K were used for aluminum and air, respectively. Once the solution to the steady state heat conduction equation was obtained, the total heat transfer rate, Q , from the top to the bottom surface was recorded. Under the assumption of onedimensional transport of heat across the interface, the TCC of the interface may be derived from an equivalent thermal resistance formulation to yield

h¼

DT Q =A

1  kdtt  kdb



ð4Þ

b

where A is the apparent area (length in 2 D) of contact. The top and bottom blocks have thicknesses dt and db (Fig. 6), and thermal con-

Fig. 5. Unstructured mesh for (a) 2D and (b) 3D (only one cutting plane shown) interface models.

ductivity values of kt and kb , respectively. The thicknesses are measured from the mean plane of the two interfaces. Since the surfaces are statistically reconstructed, in order to compute the mean TCC and its standard deviation, five ensembles were computed for each case considered in this study. 3. Results and discussion One of the difficulties associated with performing scaleresolved DNS of heat conduction is the large number of cells (mesh count) necessary to resolve all length scales. This difficulty is particularly apparent for 3D structures where, typically, 7–8 million cells are necessary to sufficiently resolve just 100 asperity pairs. Furthermore, as discussed earlier, since several ensembles must be computed for a single set of conditions, 3D calculations quickly become prohibitively time-consuming. Thus, prior to full-fledged validation and parametric studies, a study that compares 2D and 3D results for a small subset of cases was first undertaken to quantify differences in the predicted TCC values, if any. Fig. 7 shows the mean values of TCC (= h), and the statistical errors computed from the DNS of both 2D and 3D interface models, for a mean interface temperature range of 300–1000 K. For this particular study, the simulations were conducted for a steel-steel pair, with Ra values

342

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

Fig. 6. Typical computed steady state temperature distribution for a 2D interface between two rough aluminum surfaces.

3.1. Model validation

Fig. 7. Comparison of the TCC values extracted from DNS of 2D and 3D interface models with Stainless Steel (AISI-304) surfaces in contact, both with Ra ¼ 10 mm, and air as the interstitial fluid, for a mean interface temperature range of 300– 1000 K.

of 10 mm for both surfaces and zero separation distance between them. The thermal conductivity of materials was considered to temperature dependent in these simulations. For both steel (AISI304) and air, thermal conductivity increases with increase in temperature. The thermal conductivity of steel increases from 14.9 W/ m/K at 300 K to 25.4 W/m/K at 1000 K and the thermal conductivity of air increases from 0.0263 W/m/K at 300 K to 0.0667 W/m/K at 1000 K [35]. This increase in thermal conductivity is responsible for increase in the TCC at elevated temperatures. For each 2D case, the standard deviations were found to be within 11% of the mean values. Due to the complexity involved in the mesh generation process for 3D interface models, only 5 ensembles were computed for each case and the standard deviations were found to be within 17% of the mean values. Furthermore, to keep the mesh count tractable in 3D, only 100 asperities in each surface were considered. While both 2D and 3D simulations predict the same trend in the variation of TCC with interface temperature, the values computed from 3D models are consistently slightly higher—about 21% —than those computed from 2D models. These results imply that 2D models may be used for further studies with the understanding that they underpredict results of 3D models by about 20%.

For validation of the model, experimental measurements conducted by Litke [32] were considered. Litke [32] conducted experiments for various types of interfaces. The parameters varied in their experiments included contact pressure, arithmetic mean roughness and the material comprising the surface pair. For the purpose of this validation study, the contact between two aluminum surfaces with Ra values of 5.22 mm (bottom) and 5.02 mm (top) was considered. Although the contact pressure was varied between 0.5 MPa and 5.5 MPa in the actual experiment, this parametric variation cannot be captured directly in the simulations. As discussed in Section 2.1, this would require a nonlinear contact mechanics model, which is beyond the scope of the present study. Instead, to reconstruct the interface shape, values of separation distance corresponding to the experimental pressure range were obtained directly from the results of the deformation model reported by Singhal [33] (Table 1). Thus, in the present study, the separation distance is a parameter rather than the contact pressure. Furthermore, as discussed in Section 2.1, since the experimental pressure range produced very few contact spots, an indirect approach was used to validate the results. First, instead of simulating the experimental range (25.45–35.7 mm for separation distance, d), a slightly different separation distance range (14.55–27 mm) was simulated. Considering a lower value of d allowed reconstruction of interface models with statistically significant number of contacts using only 1000 asperity pairs, while still keeping the computational model tractable. Following solution of the heat conduction equation, the TCC values were extracted following the procedure discussed in Section 2.2. One of the critical pieces of information missing in the data reported by the experimental work of Litke [32] is the shape of the asperities on the surfaces. Neither the mean absolute slope, Da , nor any other shape descriptor was reported. Therefore, for the present study, we were forced to treat the mean absolute slope as a parameter. Four different values of Da were considered. Fig. 8(a) shows curve fits to the extracted TCC values for all values of Da and d considered in this validation study, along with the experimentally measured TCC values. Several important observations may be made from the results shown in Fig. 8(a). First, the agreement between the experimental TCC values and numerically predicted TCC values are quite good. This is particularly true if one considers the fact that 2D results underpredict 3D results by about 20%, as discussed earlier. Second, the effect of asperity shape is clearly evident. The TCC values are higher when the mean absolute slope is higher. This implies that more slender asperities will result in better heat conduction if the number of contacts is the same. This is because slender asper-

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

343

Fig. 9. Comparison of the numerically computed TCC with experimentally measured TCC [32] for intermediate and large separation distances.

Fig. 8. Comparison of the numerical and experimental results [32] for four different values of the mean absolute slope for the full separation distance range: (a) TCC values and (b) contact area.

ities result in larger contact area. To prove this point, the actual surface area of contact was computed for the four mean absolute slope values depicted in Fig. 8(a). Fig. 8(b) shows the contact area (normalized by the planar area of the interface) in each case. Indeed, as the asperities become more slender, the contact area increases. Third, the variation of the TCC with the separation distance is far more pronounced than its variation with the mean absolute slope. As discussed earlier and in a previous study [8], for a given Ra , the separation distance directly correlates with the number of contacts. Conversely, for a given separation distance, the number of contacts directly correlate with the Ra values of the two surfaces. To summarize, these results prove the hypothesis that the TCC is a function of both Ra and Da , i.e., both the arithmetic mean roughness height and the mean absolute slope, and that the dependence on mean absolute slope (an indicator of asperity

shape) is weaker. This issue will be further elaborated upon in Section 3.2. Although Fig. 8(a) appears to indicate almost perfect agreement between the computed values of the TCC with experimentally measured values, a closer examination indicates otherwise. Fig. 9 depicts the TCC values for a smaller range of separation distance, i.e., it depicts only a small portion of the data shown in Fig. 8(a). At intermediate and small separation distances (left end), depending on the shape of the asperities (i.e., value of Da ), the numerically obtained TCC values either overpredict or underpredict the experimentally measured values. At large separation distances (right end), no matter what the Da value, the numerical results underpredict the measured results. This may be partly due to the fact that the 2D results consistently underpredict the 3D results, as discussed earlier. In other words, in the case of 3D models, which would be a more accurate representation of the asperities, better agreement is likely. One of the notable differences between the numerical and experimental TCC values is how the TCC changes with separation distance. While both results show monotonic decrease in the TCC with d, the experimental results show a distinct reduction in the rate of change of TCC at the left end. This discrepancy may be explained as follows. In reality, when two surfaces are compressed, the asperities deform, and the material flows laterally to fill the air pockets. Initially, when the separation distance is reduced, the porosity of the interface decreases rapidly. Eventually, as the porosity approaches zero, the effect of air pockets become less and less important, and the rate of change of TCC becomes small. In the numerical model, however, no lateral flow of material is accounted for. As d is decreased, and more overlap occurs between asperities, the extra material is simply removed. Thus, the nonlinear change in porosity is not accounted for. The replication of the actual physics requires a coupled thermomechanics model, which, as discussed earlier, is beyond the scope of the present study. Nonetheless, this discrepancy warrants further investigation and may be cited as a candidate for future work. 3.2. Model reduction In Figs. 8 and 9, the TCC was expressed as a function of the separation distance, d. This is only because these TCC values were computed for a prescribed fixed value of the arithmetic mean

344

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

roughness height, Ra , while the separation distance was varied. Conversely, similar data could have been generated by fixing d, and varying Ra if one disregarded the fact that changed Ra may also have changed Da . Ultimately, the effect of d and Ra can be collapsed into a single parameter, namely the number of contacts, which, as shown in a previous study [8], is the most influential parameter for the TCC. Based on this premise, for a given material combination, the TCC may be expressed as a function of 2 independent variables as follows: h ¼ hðN c ; Da Þ, where N c is the number of contacts per 1000 asperity pairs (or non-dimensional number of contacts). Fig. 10(a) shows the mean TCC values obtained from all cases for which DNS calculations were conducted in this study. In order to generate the data for this plot, 5 ensembles were used for each case. Fig. 10(b) shows the corresponding standard deviation in the data. Clearly, when the number of contacts is small, the standard deviation in the data is large—sometimes as large as 40%. This

is because when the separation distance is large, the number of contacts is only a handful. In such cases, two different ensembles may produce—as an example—either 2 contacts or 3 contacts. However, this one extra contact can cause a large change in the computed TCC value. When the number of contacts is large, such variations have little impact on the computed TCC. Hence, for large number of contacts, the standard deviation is relatively small. The data presented in Fig. 10(a) cannot be used directly by future researchers. To make this data useful to the community at large, it was fitted using least square analysis, and a twodimensional fit was obtained. The TCC value using this fit is written as:

h ¼ ða1 þ a2 Da ÞN3c þ a3 D2a N2c þ a4 ðD2a þ N2c Þ þ a5 Da N2c þ a6 Nc Daa7 þ a8 Da þ a9 Nc

ð5Þ

where the constants a1 through a9 are tabulated in Table 2. The curve fitting toolbox in MatlabTM was used to fit the data. The inputs provided to the toolbox are the model equation and the initial values for the coefficients. From these inputs, the estimates for the coefficients are iteratively computed using the nonlinear least squares method which minimizes the summed square of residuals. The least squares formulation uses the preconditioned conjugate gradient method to solve the system of equations. Polynomial models (quadratic and cubic) were used as the starting point to fit the data which were then customized to better capture the functional dependence of TCC on Nc and Da . Several custom nonlinear models were tried and compared based on the local errors between the raw data points and predicted values, the monotonicity of the fit, the L2 norm of the fit, R-squared value of the fit, and the width of interval for the fitted coefficients. Fig. 11(a) shows a contour plot of the fitted data, while Fig. 11(b) shows the error between the fitted function [Eq. (5) or Fig. 11(a)] and the raw data [Fig. 10(a)]. It is seen that the maximum error in the fit is about 20%, which is considerably smaller than the standard deviation of the raw data itself. Therefore, the fit was deemed acceptable. As mentioned earlier, the raw data used to generate the curvefit given by Eq. (5) was obtained using DNS results from two contacting aluminum surfaces with Ra values of 5.22 mm (bottom) and 5.02 mm (top). It was also hypothesized that the actual Ra and d values are not what matters, and that the combined effect of these two may be collapsed into a single parameter, namely the nondimensional number of contacts, N c . To test this hypothesis, a sample calculation was performed, in which a completely different (from the ones used to generate the curve-fit) set of Ra , d, and Da values were chosen as inputs: Ra ¼ 7 lm, d ¼ 23 lm, and Da ¼ 3:3. The contacting surfaces were assumed to be statistically identical. First, DNS calculations were performed for these set of inputs using 5 ensembles. The DNS results were generated in order to serve as the benchmark result for testing of the accuracy of the curve-fit. From the DNS calculations, the mean TCC value was found to be 207.6 kW/m2/K, with a standard deviation of 27.7 kW/m2/K. Next, using the same inputs, the number of contacts

Table 2 Constants in the curve-fit [Eq. (5)] for TCC shown in Fig. 11(a).

Fig. 10. TCC extracted from DNS calculations for various values of the number of contacts and the mean absolute slope: (a) mean values and (b) standard deviation.

Constants in Eq. (5)

Value

Confidence bounds

a1 a2 a3 a4 a5 a6 a7 a8 a9

0.000581465710256 0.000024720202887 0.004797036796298 0.003580025784628 0.045203135751905 0.000328196832913 4.516130493687395 0.153724251910059 5.484847201963203

(0.001554, 0.0003909) (0.0002634, 0.0003128) (0.006843, 0.002752) (0.06975, 0.07691) (0.02474, 0.06566) (0.005517, 0.006174) (4.207, 13.24) (2.066, 1.759) (4.402, 6.568)

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

345

4. Summary and conclusions The TCC of an interface depends both on the size and the shape of the asperity pairs that come in contact, as well as the applied pressure. The effect of these parameters are difficult to isolate and quantify except using a high-fidelity model—one in which all length scales are resolved including the asperities that come in contact and the trapped air pockets. In this study, such a highfidelity model—referred to as scale-resolved DNS—is used to delineate the effect of size and shape of the asperities on the TCC. The model used here has the limitation that it only addresses the thermal aspects of the problem at hand, and assumes that the contact pressure (load) versus displacement behavior of the interface is already known. The microscale interface geometry was statistically reconstructed from two commonly used surface topography descriptors: (1) the arithmetic mean roughness height, Ra , which is an indicator of the asperity size and (2) the mean absolute slope, Da , which is an indicator of asperity shape. The reconstructed geometries were discretized using unstructured meshes, and the steady-state heat conduction equation was solved using the unstructured finite volume method with an applied temperature difference. The resulting heat fluxes were used to compute the TCC. First, TCC values were computed as a function of the mean separation distance (applied pressure) for a fixed arithmetic mean roughness height. The height distribution was assumed to be Gaussian, as is the case for most surfaces. Next, the same computations were repeated for different mean absolute slope values. The predicted TCC values were found to be in good agreement with experimentally measured values. Closer examination of the TCC values extracted for the entire parameter space considered in this study revealed that the effect of the arithmetic mean roughness height and the applied pressure (separation distance) could be characterized by a single parameter, namely the number of contacts. Essentially, the combination of the asperity heights and the separation distance dictates how many asperities may finally come in contact. Based on this finding, the extracted TCC data was further correlated with two nondimensional parameters: (1) the number of contacts per one thousand asperity pairs and (2) the mean absolute slope. It was found that the TCC is a strong function of the number of contacts (asperity size) and a weak function of the mean absolute slope (asperity shape), which corroborates previous findings [8]. Subsequently, the raw data was fitted to a nonlinear function using leastsquare analysis. Errors in the fit were quantified, as well. This curve-fit provides a straightforward pathway for estimating the TCC of any metal-metal interface whose surfaces are described using commonly used surface topography descriptors.

Fig. 11. Extracted TCC values after curve-fit [Eq. (5)]: (a) fitted TCC values, and (b) error in fit (from raw data shown in Fig. 10).

per 1000 asperities (N c ) was calculated using the joint probability density functions of the two surfaces, as described in Section 2.1. This yielded N c ¼ 32. The curve-fit provided by Eq. (5) was next used to compute the TCC, and yielded a value of 260 kW/m2/K. Thus, the error (between DNS results and the curve-fit) in the predicted TCC is 52.4 kW/m2/K (or 25%), which appears to be quite large. However, this error is well within the 3r range (124.5– 290.7 kW/m2/K). Furthermore, the original curve fit also had a maximum error of 20%, as shown in Fig. 11(b). Therefore, it is fair to conclude that the curve-fit provided by Eq. (5) is universally applicable—within the bounds of statistical uncertainty and fitting errors—to any surface topology input. In summary, Eq. (5) represents a valuable tool for estimation of the TCC of metal-metal contacts of arbitrary surface topology.

Acknowledgments ESI Group is acknowledged for providing licenses to their commercial software CFD-ACE+TM. The Simulation Innovation and Modeling Center (SIMCenter) at the Ohio State University is acknowledged for providing financial support. Conflict of interest statement We have no conflict of interest to disclose. References [1] C.V. Madhusudana, Thermal Contact Conductance, Springer-Verlag, New York, 1996. [2] J.A. Greenwood, Constriction resistance and the real area of contact, Brit. J. Appl. Phys. 17 (12) (1966) 1621. [3] M. Cooper, B. Mikic, M. Yovanovich, Thermal contact conductance, Int. J. Heat Mass Transfer 12 (1969) 279–300.

346

N.N. Verma, S. Mazumder / International Journal of Heat and Mass Transfer 115 (2017) 336–346

[4] B. Mikic´, Thermal contact conductance; theoretical considerations, Int. J. Heat Mass Transfer 17 (1974) 205–214. [5] A.F. Black, S.V. Garimella, Characterization of rough engineering surfaces for use in thermal contact conductance modeling, J. Thermophys. Heat Transfer 20 (4) (2006) 817–824. [6] J. Hong, J. Peng, B. Li, An integrated mechanical-thermal predictive model of thermal contact conductance, J. Heat Transfer 135 (4) (2013), 041301–041301. [7] V. Singhal, P.J. Litke, A.F. Black, S.V. Garimella, An experimentally validated thermo-mechanical model for the prediction of thermal contact conductance, Int. J. Heat Mass Transfer 48 (2005) 25–26. [8] N.N. Verma, S. Mazumder, Extraction of thermal contact conductance of metalmetal contacts from scale-resolved direct numerical simulation, Int. J. Heat Mass Transfer 94 (2016) 164–173. [9] X. Zhang, P.Z. Cong, M. Fujii, A study on thermal contact resistance at the interface of two solids, Int. J. Thermophys.: J. Thermophys. Prop. Thermophys. Appl. 27 (3) (2006) 880–895. [10] R. Xu, L. Xu, An experimental investigation of thermal contact conductance of stainless steel at low temperatures, Cryogenics 45 (2005) 694–704. [11] Q. Tang, J. He, W. Zhang, Influencing factors of thermal contact conductance between TC4/30CrMnSi interfaces, Int. J. Heat Mass Transfer 86 (2015) 694– 698. [12] R. Dou, T. Ge, X. Liu, Z. Wen, Effects of contact pressure, interface temperature, and surface roughness on thermal contact conductance between stainless steel surfaces under atmosphere condition, Int. J. Heat Mass Transfer 94 (2016) 156–163. [13] E.S. Gadelmawla, M.M. Koura, T.M.A. Maksoud, I.M. Elewa, H.H. Soliman, Roughness parameters, J. Mater. Process. Technol. 123 (2002) 133–145. [14] M. Rosochowska, R. Balendra, K. Chodnikiewicz, Measurements of thermal contact conductance, J. Mater. Process. Technol. 135 (2) (2003) 204–210. [15] S.M.S. Wahid, C.V. Madhusudana, Thermal contact conductance: effect of overloading and load cycling, Int. J. Heat Mass Transfer 46 (2003) 4139–4143. [16] C.L. Yeh, C.Y. Lin, Thermal contact resistance correlation for metals across bolted joints, Int. Commun. Heat Mass Transfer 30 (7) (2003) 987–996. [17] P. Misra, J. Nagaraju, Test facility for simultaneous measurement of electrical and thermal contact resistance, Rev. Scient. Instrum. 75 (8) (2004) 2625–2630. [18] M. Bahrami, J.R. Culham, M.M. Yananovich, G.E. Schneider, Review of thermal joint resistance models for nonconforming rough surfaces, Appl. Mech. Rev. 59 (2006) 1–12.

[19] Y.-R. Jeng, J.-T. Chen, C.-Y. Cheng, Thermal contact conductance of coated surfaces, Wear 260 (1) (2006) 159–167. [20] S. Woodland, A.D. Crocombe, J.W. Chew, S.J. Mills, A new method for measuring thermal contact conductance—experimental technique and results, J. Eng. Gas Turb. Power 133 (7) (2011) 71601. [21] V. Gopal, M.J. Whiting, J.W. Chew, S.J. Mills, Thermal contact conductance and its dependence on load cycling, Int. J. Heat Mass Transfer 66 (2013) 444–450. [22] H.A. Francis, A finite surface element model for plane-strain elastic contact, Wear 76 (1982) 221–245. [23] M.N. Webster, R.S. Sayles, A numerical model for the elastic frictionless contact of real rough surfaces, J. Tribol. 108 (3) (1986) 314–320. [24] Y.Z. Hu, K. Tonder, Simulation of 3-D random rough surface by 2-D digital filter and Fourier analysis, Int. J. Mach. Tools Manuf. 32 (1) (1992) 83–90. [25] J.-J. Wu, Simulation of rough surfaces with FFT, Tribol. Int. 33 (1) (2000) 47–58. [26] A. Majumdar, B. Bhushan, Fractal model of elastic-plastic contact between rough surfaces, J. Tribol. 113 (1) (1991) 1–11. [27] C.Q. Yuan, J. Li, X.P. Yan, Z. Peng, The use of the fractal description to characterize engineering surfaces and wear particles, Wear 255 (1) (2003) 315–326. [28] J. Lopez, G. Hansali, H. Zahouani, J.C. Le Bosse, T. Mathia, 3D fractal-based characterization for engineered surface topography, Int. J. Mach. Tools Manuf. 35 (2) (1995) 211–217. [29] J.A. Greenwood, J.B.P. Williamson, Contact of nominally flat surfaces, Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 295 (1442) (1966) 300–319. [30] T. Cui, Q. Li, Y. Xuan, Characterization and application of engineered regular rough surfaces in thermal contact resistance, Appl. Therm. Eng. 71 (1) (2014) 400–409. [31] B. Bhushan, Introduction to Tribology, Wiley, New York, 2013. [32] P.J. Litke, Experimental Determination of Thermal Contact Conductance M.S. Thesis, Purdue University, West Lafayette, Indiana, 2002. [33] V. Singhal, Prediction of Thermal Contact Conductance by Integrated Thermal and Surface Deformation Analysis M.S. Thesis, Purdue University, West Lafayette, Indiana, 2001. [34] S. Mazumder, Numerical Methods for Partial Differential Equations: FiniteDifference and Finite-Volume Methods, first ed., Academic Press, New York, 2016. [35] F.P. Incropera, D.P. Dewitt, T.L. Bergman, A.S. Lavine, Introduction to Heat Transfer, sixth ed., Wiley, New York, 2011.