Modeling a high-speed pin-on-disk experiment

Abstract

The purpose of this work is to analyze the heat transfer characteristics of Vascomax®C300 during high-speed sliding. This work extends previous research that is intended to help predict the wear-rate of connecting shoes for a hypersonic rail system at Holloman Air Force Base to prevent critical failure of the system. Solutions were generated using finite element analysis and spectral methods. The frictional heat generated by the pin-on-disk is assumed to flow uniformly and normal to the face of the pin and the pin is assumed to be a perfect cylinder resulting in two-dimensional heat flow. Displacement data obtained from the experiment is used to define the moving boundary. The distribution of temperature resulting from transient finite element analysis is used to justify a one-dimensional model. Spectral methods are then employed to calculate the spatial derivatives improving the approximation of the function which represents the data. It is concluded that a one-dimensional approach with constant heat transfer parameters sufficiently models the high-speed pin-on-disk experiment.

1. The pin-on-disk experiment

The pin-on-disk experimentation configuration, seen in Figure 1, was set up and conducted by the Air Force Research Laboratory at Wright-Patterson Air Force Base. The pin has a radius of 0.00635 meters and protrudes a length of twice the radius past the clamp as seen in Figures 1 and 2. Displacement was captured by a LDI-119 linear variable inductive transducer position sensor. Force was measured using a LC501 steel cantilever beam. Force and displacement were recorded every $5\times {10}^{-4}\mathrm{s}$ ; the raw data are shown in Figure 3. Temperature was recorded using type-J thermocouples spaced out in 1/8-inch intervals, offset by 30 degrees, and placed in holes drilled 1/5-inch deep. Type-J thermocouples are rated to 1030 K ± 0.75%. Temperature was recorded every $1\times {10}^{-2}\mathrm{s}$ ; the raw data are shown in Figure 4. Figure 4 shows that at 5.5 seconds thermocouple 1 exceeds 1100 K and there is a noticeable change in the behavior of the temperature profile. For this reason, quantitative comparisons between the models and data are only made from 0 to 5.5 seconds. The locations of thermocouple 1, 2, and 3 are shown as $T_1$ , $T_2$ , and $T_3$ in Figure 2. The pin is made from a discarded slipper shoe (Vascomax®C300) and the 12-inch diameter disk is made of the same material as the rail at the Holloman High Speed Test Track (HHSTT), AISI 4340 steel. The disk is cooled by a cold air gun during the test run to keep it from overheating, as shown in Figure 1.

Figure 1.

Pin-on-disk experiment configuration.

Figure 2.

Pin schematic.

Figure 3.

Raw force data.

Figure 4.

Raw thermocouple temperature data.

The raw data from the experiment exhibits noise and irregularities. The force was maintained manually pressing on the end of the cantilever beam. While the technician did their best to maintain a constant pressure on the pin, it is clear that vibration effects of the friction and the manual adjustments generated considerable noise to the forcing data as seen in Figure 3. The thermocouple closest to the disk (thermocouple 1) has an irregularity at 4 seconds that is unexplained. After 5.5 seconds, this thermocouple begins to read temperatures that exceed the rated threshold. Therefore, the data for thermocouple 1 cannot be trusted to be within the error tolerance after 5.5 seconds. At 7 seconds, there is an irregularity due to a drop in force and faulty readings by the thermocouple. The displacement, shown in Figure 5 is the smoothest, most regular, raw data. The smoothness of the data indicates that vibration effects did not play a large role in the collection of displacement data as it did in the forcing data. However, note in Figure 1 that the position sensor is not located in-line with the pin so that there may have been some loss of information through the connecting system.

Figure 5.

Raw displacement data.

2. Prologue

In 2008, a special run of the HHSTT was conducted. The sled system was highly instrumented to record measurements regarding the extreme aerodynamic forces the sled system undergoes. Previous work^1,2 used the 2008 test data as a baseline for the models in their research. Alban utilized a two-dimensional Finite Elements Analysis (FEA) model for the connecting shoe. Alban’s model predicted melt temperatures would be reached inside of the shoe (away from the leading edge). Alban et al. ¹ suggests that this would result in additional material wear that the model was not capable of removing. DeLeon looked to improve upon the wear removal modeling and abandoned the FEA approach, focusing instead on a one-dimensional Finite Difference Method (FDM) model. DeLeon et al. ² improved Alban’s model by considering material removal from the surface through the thickness of the shoe instead of only at the surface.

In 2019, a different approach was taken by conducting a pin-on-disk experiment to generate data for the research. ³ While the pin-on-disk experiment does not reach nearly the same velocity profiles as the sled system, data collection and modeling processes are greatly simplified. Traditional pin-on-disk experiments, such as performed by Ashby and Abulawi, ⁴ are conducted with velocities on the order of 1 to 10 meters per second. The pin-on-disk design used by Boardman et al. ³ has capability of reaching over 40 times that of the traditional test velocity, but experiments were conducted at a maximum of 240 meters per second.

This work continues the research using the data obtained by the pin-on-disk experiment. Two separate numerical solutions will be generated and discussed for the non-linear heat equation with a moving boundary. Data from the high-speed pin-on-disk experiment is used to compare each model to observation. Finite element analysis and spectral methods each provide unique approximations of, and insight to, the progression of, heat transfer through the pin. FEA will be used to show that a steady-state solution cannot effectively represent the data, leading to transient solution methods. Finite element analysis had not been previously applied to this particular pin-on-disk experimental setup, although finite elements had previously been applied to model slipper wear during a run of the hypersonic rail car. ¹ Transient FEA in two-dimensions shows convection negligibly contributes to the solution, effectively reducing the model to one-dimension. Using constant thermal properties, the transient finite element is shown to qualitatively represent the thermocouple data well. Investigation into whether variable parameters can provide a better model is performed with spectral methods. Using variable definitions for the parameters in the thermal diffusivity relation ( $\kappa =k/\rho c$ ) to represent the heat transfer is shown to not improve the model. Spectral methods were implemented to achieve high accuracy on a coarse mesh by virtue of spectral accuracy. Spectral methods have not been previously applied to any aspect of the hypersonic rail system or this particular pin-on-disk experiment. The spectral method produced the best model of the experimental data.

3. Assumptions concerning the tribology of the system

Tribology is the study of the interaction between surfaces in contact. Friction refers to the force that resists sliding. Despite having been studied since prehistoric times, friction has not received a full physical explanation, but there are two basic issues at play: the nature of the friction force and the energy dissipation mechanism. ⁵ In other words, we look to understand how much heat will enter the system due to friction, and how that heat will conduct through the system. Tables 2 and 3 provide a summary of thermal properties of the pin (Vascomax® C300) and the disk (AISI 4340) that are used in this study.

Table 2.

AISI 4340 steel material properties.

Property	Value
Density	7750 kg/m3
Heat transfer coefficient	72.14 $\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$
Thermal conductivity	41.8 $\mathrm{W}/(\mathrm{m}·\mathrm{K})$
Specific heat	505 $\mathrm{J}/(\mathrm{kg}·\mathrm{K})$

Table 3.

Finite elements analysis material removal.

Step	Displacement atstep $i,d_i$ $\left(\mathrm{meters}\right)$	Total Time(seconds)
Initial	.000000	0
1	.000000	0–1
2	.00005	1–2
3	.000125	2–3
4	.00035	3–4
5	.0007	4–5
6	.0011	5–6
7	.0015	6–7
8	.0020	7–8
9	.0022	8–9

Flux density, $q$ , is the amount of heat that enters the system due to friction per unit area. Friction is treated as a heat source at the point of contact. Flux density, or heat flux, is defined as the product of the partition function, ⁶ the coefficient of friction, pressure, and velocity divided by the contact area, ⁴ as shown in Equation (1):

$$q=\beta \mu V/A$$

(1)

The Partition Function, $\beta$ , defines the proportion of the heat generated that will enter the pin. The partition function becomes almost constant at high speeds where the Peclet number is greater than 10. ⁶ The Peclet Number is defined by $\mathrm{Pe}=\mathrm{VR}/2\kappa$ . ⁵ The Peclet number is a dimensionless quantity that is used to characterize temperature evolution. When the Peclet number is low, the temperature distribution can be represented by a stationary source; when it is high, the heat across the plane of the source can be assumed to generate uniformly. ⁶ Using the material properties given in Tables 1 the Peclet is of the order of ${10}^4$ which is large and leads to the assumption of uniform heat generation across the face of the pin. The value of the partition coefficient can be shown, as is done in Figure 6, to depend only on the Biot number of the disk. ⁶ Using Table 2, the Biot number of the disk is 0.53 since the Biot number is defined by $\mathrm{Bi}=\mathrm{hR}/k$ . ⁶

Table 1.

Vascomax® C300 material properties at room temperature.^7,8

Property	Value
Density	8000 kg/m3
Melt temperature	1685 K
Thermal conductivity	30.807 $\mathrm{W}/(\mathrm{m}·\mathrm{K})$
Specific heat	858 $\mathrm{J}/(\mathrm{kg}·\mathrm{K})$

Figure 6.

Laraqi partition function.

The coefficient of friction, $\mu$ , at high speeds is inversely proportional to velocity. ⁴ The coefficient of friction in steel sliding conditions is approximated by Equation (2) where $\tilde{V}$ is the normalized velocity given by Equation (3). ⁴ The normalized velocity is dimensionless and can be calculated using the data from Tables 1 and 2 using the diffusivity relationship.

$$\mu =0.78-0.13{\log }_{10}\tilde{V}$$

(2)

$$\tilde{V}=\frac{\mathrm{VR}}{.5({\kappa }_{\mathrm{disc}}+{\kappa }_{\mathrm{pin}})}$$

(3)

Thermal diffusivity, $\kappa$ , thermal conductivity, $k$ , density, $\rho$ , and specific heat, $c$ , are related by the diffusivity relationship. Experimental data was not found for these material properties for Vascomax®C300 in most of the temperature ranges experienced during the pin-on-disk experiment. Previous work ³ used data from Vascomax®250 to approximate thermal conductivity using the log-fit extrapolation seen in Figure 7. For specific heat function, the function seen in Figure 8 has been proposed. ² Using the thermal diffusivity relation, the function for thermal diffusivity is seen in Figure 9. ² Prior to developing this thermal diffusivity function, the constant values in Table 1 were used ¹ resulting in $\kappa =4.4883\times {10}^{-6}$ which is plotted along with the function in Figure 9. Density is assumed to be constant in all cases.

Figure 7.

Thermal conductivity of Vascomax® 250.

Figure 8.

Specific heat approximations of Vascomax® 300.

Figure 9.

Thermal diffusivity approximations of Vascomax® 300.

The clamp is assumed to be a heat sink that maintains ambient temperature. We use the same convection coefficient for the surrounding air as was previously assumed to be 25 $\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$ . ³ The pin is assumed to be a perfect cylinder which wears away in uniform sheets as to always maintain a perfect cylinder during wear. Surface zones are assumed to be removed from the disk and pin nearly instantaneously, and the materials to be Gaussian isotropic; these homogeneity assumptions result in thermal properties that are independent of position. ⁵

4. Finite element analysis

While developing the FEA model the flowchart in Figure 10 was followed. ⁹ First, a mathematical model was developed and the empirical data were obtained for subsequent comparison to the FEA results. The model was then built using FEA software ¹⁰ with initial discretization chosen based upon the thermocouple and displacement locations. We use the displacement data at 1 second intervals to define a non-linear problem with a moving boundary for the FEA by moving the flux boundary inward and removing the elements to model the material being worn away. Several outputs are displayed for examination after solutions were found. The mesh was refined until the refinement no longer significantly altered the results.

Figure 10.

FEA flowchart.

The assumption of uniform heat flux across the face of the pin is sufficient to reduce the pin model to two dimensions. As a result, every plane which is normal to the pin face and cuts through the center of the pin face, as shown in Figure 11, will have the same temperature profile. As a result, the pin-on-disk can be viewed as a two-dimensional plane heat transfer problem. Redrawing the pin from Figure 2 under all the assumptions produces the schematic shown as the two-dimensional plane object in Figure 12 where the three thermocouples have been projected onto the same plane.

Figure 11.

Planes through the center, and normal to, the face of a cylinder.

Figure 12.

Plane model of the pin.

In order to double check that the model has been correctly set up in the software, and in order to properly serve the first three boxes of our flowchart, a one-dimensional steady-state model was considered and compared against the computer model. Starting from conservation of energy the basic finite elements matrix equation is given in Equation (4). ¹¹ The heat rate (heat flow per unit time) is given in Equation (5). ¹¹

$$\left[\mathit{K}\right]\left\{\mathit{T}\right\}=\left\{\mathit{Q}\right\}$$

(4)

$$Q=-\mathrm{kA}\frac{\mathrm{dT}}{\mathrm{dx}}$$

(5)

Reducing Figure 12 to 1d and solved the steady-state problem using four elements the problem would reduce to Equation (6) containing five nodes. Assuming that the lengths of the elements are equal, so that the interior nodes represent the three thermocouples; constant thermal conductivity; $T_4$ is held to ambient temperature $T_{\infty }$ ; $Q_0$ is the frictional heat flux source q (from Equation (1) times $A$ , and $Q_1=Q_2=Q_3=0$ (no generation or loss); and, let the unknowns be $T_0$ , $T_1$ , $T_2$ , $T_3$ and $Q_4$ . Then Equation (6) is completely defined, but the matrix is singular. Since $T_4$ is known, one method to remedy this is to impose the boundary condition and move the known quantities to the right-hand side as in Equation (7). ⁹ Now the system can be solved by performing Gaussian elimination ¹² on Equation (7) to find the temperatures given in Equation (8). With all the temperatures known, $Q_4$ can be readily solved by going back to Equation (6), as shown in Equation (9).

$$\left[\begin{array}{ccccc}1\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& 2\hfill & -1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -1\hfill & 2\hfill & -1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -1\hfill & 2\hfill & -1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -1\hfill & 1\hfill \end{array}\right]\left\{\begin{array}{c}{T}_0\hfill \\ \multicolumn{1}{c}{T_1}\\ \multicolumn{1}{c}{T_2}\\ \multicolumn{1}{c}{T_3}\\ \multicolumn{1}{c}{T_4}\end{array}\right\}=-\left\{\begin{array}{c}\mathrm{qA}\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{Q_4}\end{array}\right\}\frac{L/4}{\mathrm{Ak}}$$

(6)

$$\left[\begin{array}{ccccc}1\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& 2\hfill & -1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -1\hfill & 2\hfill & -1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -1\hfill & 2\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\left\{\begin{array}{c}{T}_0\hfill \\ \multicolumn{1}{c}{T_1}\\ \multicolumn{1}{c}{T_2}\\ \multicolumn{1}{c}{T_3}\\ \multicolumn{1}{c}{T_4}\end{array}\right\}=-\left\{\begin{array}{c}\mathrm{qL}/4k\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{T_{\infty }}\\ \multicolumn{1}{c}{T_{\infty }}\end{array}\right\}$$

(7)

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\left\{\begin{array}{c}{T}_0\hfill \\ \multicolumn{1}{c}{T_1}\\ \multicolumn{1}{c}{T_2}\\ \multicolumn{1}{c}{T_3}\\ \multicolumn{1}{c}{T_4}\end{array}\right\}=\left\{\begin{array}{c}{T}_{\infty }-4\mathrm{qL}/4k\hfill \\ \multicolumn{1}{c}{T_{\infty }-3\mathrm{qL}/4k}\\ \multicolumn{1}{c}{T_{\infty }-2\mathrm{qL}/4k}\\ \multicolumn{1}{c}{T_{\infty }-\mathrm{qL}/4k}\\ \multicolumn{1}{c}{Q_4}\end{array}\right\}}\\ \multicolumn{1}{c}{\approx \left\{\begin{array}{c}2750\hfill \\ \multicolumn{1}{c}{2135}\\ \multicolumn{1}{c}{1520}\\ \multicolumn{1}{c}{910}\\ \multicolumn{1}{c}{296}\end{array}\right\}\mathrm{Kelvin}}\end{array}$$

(8)

$$Q_4=-\frac{\mathrm{Ak}}{L/4}(T_4-T_3)\approx 755\mathrm{Watts}$$

(9)

The part in Figure 12 was first partitioned by displacement data at 1 second intervals and thermocouple locations producing the partitions in Figure 13. Partitioning at thermocouple locations ensures that when the part is meshed there will be a node located at the same geometric location as the thermocouples. Therefore, the values at these nodes represent the interpolated approximation of data from the experiment. The values from the data for the displacement partitions in 1 second intervals are given in Table 3. Partitioning at displacement locations accounts for the wear: eight steps were created and after each step the model is changed so that that the elements in the lowest partition are removed from the model to represent the wear. The shape functions were generated using four-node linear heat transfer bilinear quadrilaterals when defining the mesh for the part.

Figure 13.

Partitions of the plane model.

The FEA model was generated using the constant material properties that were defined using Table 1. The heat flux properties are summarized in Table 4 using values from averaging data, or form the assumptions section, as applicable. The plate is given a thickness of $0.00997455\mathrm{m}$ so that the surface area acted upon by $q$ is the same as the pin face. This creates an issue for the convection boundary with the air though. The surface area of the sides of the cylinder is $2R\left(2\pi R\right)=5.067\times {10}^{-4}$ , but the area on the sides of the plate is $2\left(2R(.00997455)\right)=2.534\times {10}^{-4}$ which is smaller by a factor of 2. Therefore, a surface film condition of $50\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$ , instead of 25 as shown in Figure 12, is applied to the sides of the plane to correct the heat rate.

Table 4.

Constant flux parameters.

$\beta$	0.12
$\mu$	0.20
$F$	656.76 N
$V$	47.83 m/s
$A$	0.0001267 m2

Steady-state solutions were obtained to validate the computer model set up versus the results of Equation (9) and to test for convergence. A coarse mesh containing 96 elements and 117 nodes, a fine mesh containing 728 elements with 783 nodes, and a control mesh of 899 elements with 960 nodes were used to check convergence of the solution. Convergence of the mesh is achieved for the fine mesh of 728 elements. The maximum and minimum values along the bottom boundary were used to determined convergence. The difference in maximum and minimum temperatures between the control and fine mesh are within .0001% of each other, as shown in Figures 14 and 15. The fine mesh is used for all subsequent analysis.

Figure 14.

End state of the steady-state FEA model with 793 nodes.

Figure 15.

End state of the steady-state FEA model with 960 nodes.

For transient solutions the model seen in Figure 12 is used. Transient steps were implemented in Abaqus ¹⁰ using 1 second increments per step, with a moving heat flux boundary corresponding to the displacement data in Table 3. At the end of each 1 second interval the bottom-most remaining partition is removed and the flux is moved to the new bottom-most surface to model the wear. This creates a two-dimensional non-linear, moving boundary, transient FEA model of the pin.

The initialized two-dimensional FEA model is shown in Figure 16. The elements surrounding the three thermocouple locations have been highlighted for reference. Figure 17 shows the removal of the bottom seven partitions and the final state of the model at 8 seconds. The nearly horizontal isotherms were noted, and investigation into the contributions of the convection with the surrounding air was conducted by removing the surface films from the model, effectively creating a one-dimensional representation. The temperature profiles of the two-dimensional and one-dimensional model are shown in Figure 18. The average relative difference between models was found to be less than 1% and so the convection was determined to be negligible in this study.

Figure 16.

Initial state of the two-dimensional FEA model (t = 0 s).

Figure 17.

End state of the two-dimensional FEA model (t = 8 s).

Figure 18.

Temperature profiles of thermocouples for 2-D and 1-D FEA models.

The temperature profile of the one-dimensional model is plotted against the data in Figure 19. The model generally trends well with the data, nearly tracing the data for thermocouple 2, but that the model exceeds $T_{\mathrm{melt}}$ at the surface. Exceeding the melt temperature at the surface does not make physical sense, so to correct the partition function was chosen to be $\beta =0.11$ , instead of $0.12$ as was determined in the assumptions section.

Figure 19.

Temperature profiles of thermocouples for 1-D FEA model and the data.

There are many factors involved that could have been changed to lower the temperature at the surface for the model, but this change was taken because it is supported by the way the partition function is trending in Figure 6 and the large Peclet number for this experiment. The temperature profile for this alternate one-dimensional model is given in Figure 20. While thermocouple 2 is no longer traced by the data, the average relative difference between the model and the data is smaller and the temperature at the surface no longer exceeds melt. Table 5 provides a summary of FEA computed average relative difference at thermocouples 1, 2, and 3 labeled here as T₁, T₂, and T₃, respectively.

Table 5.

Summary of finite elements analysis relative differences.

Set 1	Set 2	$T_1$ relativedifferences (%)	$T_2$ relativedifferences (%)	$T_3$ relativedifferences (%)	Avg. relativedifferences (%)
Steady state	Data	200.2	202.3	155.0	185.8
Transient 2-D	Data	11.7	1.0	8.3	7.0
Transient 1-D	Data	13.0	1.8	9.0	7.9
Transient 1-D Alt	Data	6.3	2.7	7.0	5.3
Transient 2-D	Transient 1-d	1.2	1.1	0.6	0.9

Figure 20.

Temperature profiles of thermocouples for the alternate 1-D FEA model and the data.

5. Spectral methods

After determining it appropriate to reduce the model to one-dimension, we decided to implement spectral methods primarily for efficiency. Diffusion of heat through an initially constant temperature cylinder is a problem that results in smooth data on a simple domain. For problems on simple domains, with smooth data, spectral methods provide good approximations. ¹³

The general form of the heat diffusion equation is given by Equation (10). ¹⁴ Transforming to cylindrical coordinates, ³ and reducing to one-dimension, results in Equation (11). Redefining the length by Equation (12), and the temperature by Equation (13), results in the scaled heat equation given by Equation (14). The MATLAB cheb command generates a derivative matrix, $D_N$ , by interpolating at $N+1$ Chebyshev points. ¹³ By rescaling the heat flux by Equation (15), ³ and using the spectral differentiation matrix, $D_N$ , Equation (14) can be written in the form of Equation (16) with the boundary conditions given by Equations (17) and (18).

$$\begin{array}{c}\multicolumn{1}{c}{\rho c\frac{\partial T}{\partial t}=\nabla ·(k\nabla T)=k(\nabla ·\nabla T)+\nabla k\nabla T}\\ \multicolumn{1}{c}{=k{\nabla }^{2}T+\left(T\right)\frac{\mathrm{dk}}{\mathrm{dT}}{|\nabla T|}^2}\end{array}$$

(10)

$$\frac{1}{\kappa }\frac{\partial T}{\partial t}=\left(\frac{{\partial }^{2}T}{\partial z^2}\right)+\frac{T}{k}\frac{\mathrm{dk}}{\mathrm{dT}}\left[{\left(\frac{\partial T}{\partial z}\right)}^2\right]$$

(11)

$$x=\frac{z}{2R};0\le x\le 1$$

(12)

$$T_s=\frac{T-T_0}{T_{\mathrm{melt}}-T_0}$$

(13)

$$\frac{R^2}{\kappa }\frac{\partial T_s}{\partial t}=\frac{{\partial }^2{T}_s}{\partial x^2}+\frac{T_s}{k}\frac{\mathrm{dk}}{d{T}_s}\left[{\left(\frac{\partial T_s}{\partial x}\right)}^2\right]$$

(14)

$$q_s=\frac{\mathrm{Rq}}{k(T_{\mathrm{melt}}-T_0)};\mathrm{at}x=0$$

(15)

$$\frac{\partial T_s}{\partial t}=\frac{\kappa }{R^2}[{\left(D_N\right)}^2{T}_s+\frac{T_s}{\kappa }\frac{\mathrm{dk}}{d{T}_s}\left[{\left(D_N{T}_s\right)}^2\right]]$$

(16)

$$T_s=0;\mathrm{at}x=1$$

(17)

$$T_s=\frac{1}{{(D_N)}_{N+1,N+1}}\left[-q_s-{(D_N)}_{N+1,1:N}{\left(T_s\right)}_{1:N}\right];atx=0$$

(18)

The schematic for the scaled one-dimensional spectral methods model of the pin is shown in Figure 21. The system given in Equations (16)–(18) is completely defined so that the derivative in time can be solved using finite-difference methods. For initial solutions, the MATLAB built-in stiff differential equation solver ode15s ¹⁵ was called since the heat equation is known to exhibit stiffness. ¹⁶ The moving boundary was defined in the same manner as in the FEA section using the displacement function given by the data in 1 second increments, provided in Table 3. However, the derivative matrix, $D_N$ , always assumes the data are in the domain [−1,1].

Figure 21.

Scaled one-dimensional model of the pin.

In order to satisfy the domain requirements and implement the solution, we rescale each $d_i$ from Table 3 into the rescaled relative displacement points, $r_i$ in Table 6. The rescaled relative displacement points were found by first finding the relative displacement, $d_i^{*}=d_i-d_{i-1}$ , then finding the rescaled relative displacement by, $r_i=2d_i^{*}/(1-d_i)$ . After each 1 second run the temperature is interpolated using the chebfun routine ¹⁷ at the points $\stackrel{\rightharpoonup }{C}\left(2-r_i\right)/2-r_i/2$ , where $\stackrel{\rightharpoonup }{C}$ is the vector of Chebyshev points, $C_j=\cos (j\pi /N)$ . ¹³ So that the solutions are computed on the [−1,1] Chebyshev grid and then re-scaled back onto the $x\in [0,1]$ grid.

Table 6.

Spectral method material removal.

Step	Relative rescaleddisplacement at step $i,r_i$	Total time(seconds)
Initial	.000000	0
1	.000000	0–1
2	.0078	1–2
3	.012	2–3
4	.0366	3-4
5	.0582	4–5
6	.069	5–6
7	.0714	6–7
8	.0936	7–8
9	.038	8–9

The same constant properties were used as in the one-dimensional FEA case provided in Tables 1 and 3. The rescaled displacements, which model the wear of the pin, and the modeled temperature through the pin are visualized by Figure 22. The resulting temperature profiles of the thermocouples and surface temperature are shown in Figure 23. The average relative difference between the spectral methods model with constant parameters and the data are 6.6%, which is an improvement from the FEA models that also used properties from Tables 1 and 3.

Figure 22.

Wear and temperature profiles through the scaled 1-D pin using spectral methods.

Figure 23.

Temperature profiles of thermocouples, surface, and the data using spectral methods.

Variable parameters were also considered by using the fifth order averaged interpolant for the forcing data (shown in Figure 24) and the functional estimate of the thermal conductivity seen as the curve fit of Figure 7. ² Thermal diffusivity was held constant at $4.4883\times {10}^{-6}$ , which is the result obtained using the data from Table 1. The average relative difference between the variable parameter model and the data at thermocouples increased to 6.9%. Since the variable parameter model fits the data slightly worse compared the constant parameter model, it is determined that considering the variations of parameters is not advantageous to this spectral methods model.

Figure 24.

Polynomial interpolant of the raw force data.

Recall that the FEA model was improved by decreasing the partition function, but note that the spectral method does not benefit in the same way. While the FEA model was seen in Figure 19 to be overshooting the data, Figure 23 shows the spectral method model to be generally undershooting the data. To find the model parameters that would make a better fit the data two main factors were considered: how much heat enters the system through friction which was controlled by the value of the product of $\mu \&\beta$ which directly effects the heat flux in Equation (1); and how that heat diffuses through the system by controlling the value of $\kappa .$ Since it was determined that the variations of parameters are not advantageous to the model, and to greatly simplify the discovery process, $\mu \beta$ and $\kappa$ were considered constants.

A two-dimensional optimization routine was created using the MATLAB function fmincon ¹⁸ to determine the optimal values of $\mu \beta$ and $\kappa$ . The optimal values found were $\mu \beta =0.0541$ and $\kappa =2.5877\times {10}^{-6}$ . The resulting temperature profiles at the thermocouples, plotted in Figure 25, have an average relative difference to the data of 5.2%, the best of any method here-in. However, the optimized spectral method has a surface profile, which can also be seen in Figure 25, that exceeds the melting temperature and therefore does not make physical sense.

Figure 25.

Temperature profiles of thermocouples, surface, and the data using optimized spectral methods.

Finally, a convergence study was conducted to investigate the stiffness of Equation (16). Table 7 shows the exponential increase of time-steps and overall time required to solve for the first 1 second of the solution using the MATLAB solver ode45. ¹⁹ Equation (16) is determined to be stiff which results in an out of memory error by $N=160$ . Solution convergence using the solution at $N=150$ as the reference is plotted in Figure 26. From Figure 26, we see that near full convergence can be obtained by using $N=36$ .

Table 7.

Time-steps and time required to solve equation (16) using ode45.

N	Time-steps	Time (s)
5	53	0.01
10	93	0.01
20	1041	0.04
40	16317	0.75
80	260437	9.38
160	N/A a	N/A a

a

MATLAB Error: Out of memory.

Figure 26.

Convergence of ode45.

All previous solutions presented used ode15s with $N=150$ to avoid the stiffness and ensure convergence. However, this approach was determined to not be required since the average relative difference between the ode15s solution using $N=150$ and the ode45 solution using $N=36$ is less that 1% and is plotted in Figure 27. Table 8 provides a summary of spectral method computed average relative difference at thermocouples 1, 2, and 3 labeled here as T₁, T₂, and T₃, respectively. Where ode15s with $N=150$ took about 1.7 s to compute the full solution, ode45 with $N=36$ took 0.5 s to compute the full solution.

Table 8.

Summary of spectral method relative differences.

Set 1	Set 2	$T_1$ relativedifferences (%)	$T_2$ relativedifferences (%)	$T_3$ relativedifferences (%)	Avg. relativedifferences (%)
Constant parameter	Data	8.4	5.8	5.7	6.6
Variable parameter	Data	7.8	8.6	4.3	6.9
Optimized parameter	Data	2.7	12.7	0.3	5.2
Constant parameter	Variableparameter	5.5	3.2	0.3	3.4
ode15s (N = 150)	ode45 (N = 36)	1.4	0.7	0.5	0.9

Figure 27.

Temperature profiles of thermocouples for the ode15s with N = 150 and ode45 with N = 36 Models.

6. Conclusions

Two separate numerical solutions were generated to solve the non-linear heat equation with a moving boundary, and empirical evidence from a high-speed pin-on-disk experiment was used to provide validation. Finite element analysis and spectral methods each provided unique opportunities to inspect the progression of heat transfer through the pin. The frictional heat source was assumed to be generated uniformly across the pin face and the opposite end of the pin was assumed to maintain the temperature of the ambient air. Finally, assuming uniform convective heat loss off the sides of the pin, the model reduced to two dimensions. Finite element analysis was used to first show that the steady-state solution could not effectively represent the empirical evidence (Equation (8) and Figure 14) so that a transient solution would be required. The transient finite element analysis model then provided insight that the pin experiences negligible heat loss with the surrounding air (Figure 18). Negligible convection effectively reduces the model to one-dimension. The one-dimensional finite element model qualitatively represents the thermocouple data well, as exhibited by the tracing of thermocouple 2 in Figure 19. However, the model over-represents the temperature profile along thermocouples 1 and 3. The finite element analysis was completed using only constant parameters, and analysis into the effect of variable parameters was desired. With the reduction to one-dimension, it was believed that switching to a different solution method would provide a more efficient form to interrogate variable parameters.

Spectral methods, formed via the differentiation matrix, were implemented. As an added benefit, solution times using spectral methods were reduced to about 1 second. The solution profiles of the temperature at the thermocouples using spectral methods qualitatively represent the empirical data well, as exhibited by Figure 23. While thermocouple 2 was virtually traced by the FEA method, the overall representation of the three thermocouples is modeled best by the spectral methods. The introduction of variable parameters does not improve the model qualitatively or quantitatively when using spectral methods, so it is concluded that including variable parameters is not advantageous for this pin-on-disk experiment. Optimizing the parameters provided a $\kappa$ that is about half of the assumed constant and a $\mu \beta$ that is about twice that of the assumed values.

Finite elements analysis and spectral methods both provided good qualitative temperature profiles for the thermocouple temperatures using one-dimension and constant parameters. Spectral methods provided the best quantitative measure of fit by producing an average relative error of 6.6% versus 7.9% using finite elements. Implementing spectral methods for the space derivative also provided the ability to avoid the stiffness presented by the mathematical model; achieving spectral accuracy with a course mesh does not present an issue for the computational efficiency or accuracy of explicit methods for the time derivative.

For future work, extended FEA models and radial basis functions are suggested. The assumptions that were made to reduce the cylindrical pin to a two-dimensional plane model were uniform heat generation across the face of the pin and uniform radial convection with the surrounding air. By extending the finite element model to three dimensions it may be possible to determine if any aerodynamic effects of the spinning wheel significantly affect the model; and, if non-uniform frictional heat generation ⁵ plays a significant role. The finding that reduced the two-dimensional plane model to one-dimension was that radial convection was negligible in the finite element analysis. Radial basis functions maybe used to generate a derivative matrix with high accuracy and can be easily extended to higher dimensions ²⁰ so that Equation (16) can be solved in two or three dimensions using derivative matrices similarly. Radial basis function-finite-difference (RBF-FD) methods can reduce the derivative matrix to a sparse matrix in the interest of efficiency. ²⁰ Since the size of the matrix will grow exponentially in higher dimensions, it may be of benefit to explore if the computational cost can be reduced by finding an appropriate number of nearest neighbors while adding negligible errors to the solution.