Heat Transfer Analysis for Non-Contacting Mechanical Face Seals Using the Variable-Order Derivative Approach

: This article presents a variable-order derivative (VOD) time fractional model for describing heat transfer in the rotor or stator in non-contacting mechanical face seals. Most theoretical studies so far have been based on the classical equation of heat transfer. Recently, constant-order derivative (COD) time fractional models have also been used. The VOD time fractional model considered here is able to provide adequate information on the heat transfer phenomena occurring in non-contacting face seals, especially during the startup. The model was solved analytically, but the characteristic features of the model were determined through numerical simulations. The equation of heat transfer in this model was analyzed as a function of time. The phenomena observed in the seal include the conduction of heat from the ﬂuid ﬁlm in the gap to the rotor and the stator, followed by convection to the ﬂuid surrounding them. In the calculations, it is assumed that the working medium is water. The major objective of the study was to compare the results of the classical equation of heat transfer with the results of the equations involving the use of the fractional-order derivative. The order of the derivative was assumed to be a function of time. The mathematical analysis based on the fractional differential equation is suitable to develop more detailed mathematical models describing physical phenomena.


Introduction
Over the last few years, there have been many studies using fractional-order differential and integral operators to generalize classical differential and integral calculus with the aim of further understanding the nature of complex systems. Currently, attempts are being made to apply fractional calculus to solve various physical, mechanical, biological, or chemical problems. While classical integer-order operators are dependent only on the local behavior of the function, fractional-order operators accumulate all the information about the function. Another fundamental feature of fractional derivatives is that they are defined along a segment, not at a point, as is the case with classical derivatives. This feature ensures a more accurate and effective analysis of different phenomena. One of the shortcomings that most models have to overcome is that, mathematically, velocity is an instantaneous velocity defined at a point. Thus, it can be seen from the literature that differential calculus is used in the heat transfer theory.
It has been found, for example, that variable thermal conductivity is a key physical property of materials, especially when it is dependent on temperature. Variable thermal conductivity is of significance in a wide range of applications, including modern physics and mechanical engineering. It is taken into consideration primarily to determine the effect of temperature on the performance of machine elements. Special attention is paid to this property when rapid changes in temperature occur, as they may affect the operation of mechanical assemblies and subassemblies, and consequently the whole machine. The relationship between variable thermal conductivity and fractional differential calculus has been investigated extensively for problems related to fluid mechanics [1,2], viscous elasticity [3], relaxation [4], thermal elasticity and conductivity [5,6], control theory [7,8], and many others [9][10][11]. Over the years, researchers have extended the classical theory of elasticity and introduced the theory of thermoelasticity.
The mechanical properties, especially the relaxation, of additively manufactured materials were analyzed, for example, in [4]. A mathematical model was developed using fractional differential calculus to further fit the relaxation curves with the experimental data.
Warbhe et al. [5] used the quasi-static approach to the fractional-order theory of thermoelasticity to solve a two-dimensional problem for a thin circular plate with zero temperature across the lower surface and constant and evenly distributed temperature across the insulated upper surface. The integral transform technique was employed to solve the physical problem, while the displacement potential function was applied to calculate thermal stresses. Povstenko and Kyrylych [6] studied the interface between two solids. They considered generalized boundary conditions of a non-ideal thermal contact to solve the heat conduction equation of the fractional order as a function of time using the Caputo derivative.
The extensive monograph by Kaczorek [7] presented the theoretical and practical aspects of the application of fractional calculus to fractional discrete-time linear systems.
The research by Nowacki, described for instance in [12][13][14], was a breakthrough in this field. His findings are still valid and thus widely cited, as they can be applied to solve a large number of theoretical and practical engineering problems.
Knowledge of all theories available in this area is crucial to understand the deformations of elastic materials. Their thermal and mechanical behavior needs to be taken into account whenever thermoelastic deformations result in leaks, as is the case with noncontacting face seals, where the gap height may range from one to several micrometers.
Fractional differential calculus has proved well-suited to predict many physical phenomena associated with elastic media, e.g., thermal conductivity, heat transfer, and viscoelasticity. The classical differential equations describing physical phenomena are modified using the variable-order derivative (VOD) time fractional approach. Recently, there has been much research focusing on extending the applications of fractional differential calculus, e.g., [15][16][17][18]. These studies provide an insight into many important aspects associated with heat conduction. In [15], for example, Povstenko presented fundamental solutions to the time-fractional advection diffusion equation for two cases: the plane and the half-plane. He also considered Cauchy problems, inverse source problems, and Dirichlet problems. The solutions were expressed in terms of Bessel integral functions combined with Mittag-Leffler functions. In his other articles [16,17], Povstenko discussed solutions for heat transfer in composite media, in which he employed the time-fractional heat conduction equation with the Caputo derivative of fractional order to describe heat transfer in both constituent materials (0 < α ≤ 2 and 0 < β ≤ 2, respectively). The problem was solved under ideal contact conditions. This suggests that the temperature of the two materials and the heat fluxes were the same.
Liu et al. [19] proposed an approach to solving the time fractional nonlinear heat conduction equation. Similar solutions were proposed by Koca and Lotfy [20,21].
Another issue is an inverse problem in heat transfer. Liu and Feng [22], for instance, solved it using the fractional differential equation for a time-dependent derivative. The solution to the 2D time-fractional inverse problem of diffusion was based on an improved kernel technique. Maximum a posteriori estimation was required. Convergence was achieved using the regularization term and deriving a priori probability.
The primary aim of this article was to show how the results obtained by solving the classical equation of heat transfer differ from those obtained for the same physical phenomenon by employing the fractional differential equation. The study involved using the derivative order as a time-dependent function. The results from the two approaches were compared graphically. The analytical solutions of the main physical quantities were developed using the Marchi-Zgrablich transform, the finite Fourier cosine transform, and, finally, the Laplace transform.

Mathematical Model
The classical Fourier heat equation can be written as: where q-heat flux vector, ϑ = T − T o and λ-thermal conductivity. When combined with the energy conservation principle, the heat equation can be used to calculate the local deformation: with ρ being the density and C p the specific heat capacity. Many researchers have examined the application of differential calculus or integral calculus involving the use of derivatives of fractional orders. Fractional calculus is a natural extension of the notions of differentials and integrals used in classical differential calculus and integral calculus, respectively. Recently, there has been much interest in the use of fractional differential calculus to explain and model many physical phenomena. For example, fractional differential calculus has been used to design and/or model PI λ D µ controllers, heat conduction [23][24][25], thermoelasticity [26], complex nonlinear systems [7], supercapacitors, electrical and mechanical systems [27][28][29], electrical filters, dielectric relaxation, diffusion, and viscoelasticity [30].
The time fractional heat conduction equation for the rotor can be written as [31]: For the case considered here, it was assumed that α = α(t). The function α(t) was defined as a function of time (see Figure 1). The function describing the change in the coefficient α(t) can be written using the following formula: The transition to the time fractional heat conduction equation with the boundary conditions, was described, for example, in [32] as: where L is the thickness of the sealing rings. The Caputo derivative of the fractional order can be defined as described in [15]: where Γ(α) is the gamma function, with the initial conditions being: The heat conduction equation for the stator is thus rewritten as: The heat transfer model analyzed here also needs to take into consideration the heat flux generated in the gap, which can be described with a relationship based on the simplified energy equation. A similar method was used in [32]: The velocity of the fluid particles ν ϕ is linearly variable. Ranging from zero on the stator surface to the value of (ω r) on the rotor surface, it can be described as: The energy Equation (8) was solved taking into account relationship Equation (9) to calculate the distribution of temperature in the gap.
Considering sealing rings with unmodified face surfaces and neglecting other factors affecting their geometry, e.g., mechanical deformations, we can assume that the height of the gap is constant: h = const. The gap height considered by other researchers is that discussed, for example, in reference [33].
As dynamic viscosity is a temperature-dependent parameter, here it can be defined using the formula proposed by Li et al. [1]: The average temperature of the medium in the gap can be determined from: Equation (10) is a function describing the distribution of dynamic viscosity µ(r) in the radial direction.
The seal performance is largely affected by the forces acting on the stator and the rotor. The closing force produced by the spring ensures the leak tightness of the device in the OFF mode and the contact of the stator and the rotor. However, when the device operates (is ON), the opening force is generated by the pressure of the medium in the gap. The distribution of pressure can be described using the one-dimensional Reynolds equation: where the boundary conditions are: Equation (12) is the simplified Reynolds equation describing the changes in pressure of the incompressible medium in the gap only in the radial direction.

Boundary Conditions
The schematic diagram in Figure 2 illustrates the non-contacting face seal with the flexibly mounted stator (1) and the rotor (2) connected to the shaft (6) of the turbo machine. The gap between the stator and the rotor, with a height ranging from one micrometer to several or even more than ten micrometers, is filled with a sealing medium, e.g., water. The distribution of temperature on the stator and rotor surfaces is dependent mainly on the heat flux in the gap.  In this analysis, it is assumed that the inner and bottom surfaces of the rotor (r i ) and the inner and top surfaces of the stator are not in contact with the surroundings; thus, heat transfer for these surfaces takes the general form: ∂ϑ ∂n = ∂θ ∂n = 0. In a real system, these surfaces are in contact with other elements of the seal characterized by different physical properties. Under certain conditions, heat transfer taking place between these elements can be assumed to be negligible. Another reason to introduce the boundary conditions ( ∂ϑ ∂n = 0) for the surfaces is the considerable simplification of the calculations for the analytically solved model.
As the gap is limited by the ring faces, the boundary conditions for the rotor and the stator are as follows: For the outer surface of the rotor (r o ), heat transfer is assumed to be by convection, and it can be expressed by: where α f is the convection coefficient. All the above boundary conditions are necessary to solve the system of three differential equations.

Problem Solution
For the case considered here, a cylindrical coordinate system corresponding to the geometry of the physical model was used, and thus Equation (3) can be written as: Equation (18) was solved by adopting the Marchi-Zgrablich transform, the finite Fourier cosine transform, and the Laplace transform. The general form of the finite integral transform by Marchi-Zgrablich [34] is: The inverse integral transform proposed by Marchi-Zgrablich can be expressed as: where a n = f p (n) Thus, Equation (18) is written as: The finite Fourier cosine transform is described in the general form as: where m = 1, 2, 3, . . .
The inverse transform is given in the general form as: Applying the finite Fourier cosine transform to write Equation (23) and assuming the boundary conditions given in Equation (14), we obtain: With the Laplace transform [25] being: we can rewrite Equation (27) as: Once the transformations are completed, Equation (18) has the following form: where m 2 π 2 + L 2 k 2 n L 2 = ω 2 and ω 0 After the three inverse integral transforms are employed, the equation describing the temperature distribution for the rotor takes the final form: where E α,β (z) is the function of the Mittag-Leffler type [35,36]: For α, β = 1, E α,β (z) = e z , and thus Equation (31) is solved using the classical heat equation with predetermined initial and boundary conditions. S p (λ, α, k n · r) = (−k n Y 1 (k n r i ) + α Y 0 (k n r o ) − λ k n Y 1 (k n r o ))J 0 (k n r) +(k n J 1 (k n r i ) − α J 0 (k n r o ) + λ k n J 1 (k n r o ))Y 0 (k n r) Equation (33) is dependent on the boundary conditions assumed for the rotor or stator model. The distribution of temperature in the stator is defined as: where θ f is the excess temperature of the medium in the gap, calculated as θ f = T f − T o .

Operating Conditions
The parameters characterizing the material used for the rotor and the stator are provided in Table 1. The parameters concerning the seal geometry and performance are shown in Table 2. It is assumed that the surrounding fluid is water with a temperature of 20 • C. The analytical solution to the mathematical model considered here shows the distribution of temperature in the cross-section of the seal when α in Equation (18) has a fractional order value of 1 ≤ α ≤ 2. Moreover, the parameter α is time-dependent, which is illustrated graphically in Figure 1.

Results and Discussion
The simulations were performed to determine the influence of the factor α on the distribution of temperature in the seal (Figure 4).   (Figure 4a), the temperature of the rotor increases gradually because of the heat flux generated in the gap. The greatest difference between the diagrams is observed for the time range 0.1-1.25 s. With the assumption that α(t), Figure 1 suggests that if values higher than t 0 = 1 s are used, the order of the derivative tends toward unity.
The heat flux provided by the sealing ring face at a given moment of time for t = 0.05 (s) causes a local increase in temperature of 2.1 • C, which may result in thermal deformations of the sealing ring in contact with the gap. In the case of the fractional differential equation (for t = 0.05 (s)), the equation is hyperbolic in nature, as illustrated in Figure 5b. The solution based on the variable-order derivative (VOD) time fractional model provides a better correlation between the predicted data and the observed results than the classical approach.

Validation of Results
The results were validated using Ansys Workbench software ( Figure 6).   From Figure 7, it is apparent that the results obtained using both calculation methods are in good agreement.

Conclusions
This article has presented a mathematical model to determine the heat transfer phenomena occurring in the non-contacting face seal used in a turbo machine. The problem was solved analytically by applying three integral transforms. The resulting relationships described the changes in temperature in the face seal cross-section. The boundary conditions were taken into consideration. The classical Fourier law was extended by applying the time fractional derivative where α = α(t). This model can be used to describe unstable heat transfer conditions or thermal shock.
Even small local changes in temperature of an order of 2-10 • C can cause thermoelastic deformations of the stator and the rotor. As a result, there is a change in the gap geometry. During the first moments of the startup, unstable conditions can be observed because the seal operates under dry friction conditions; thus, the heat flux generated between the rotor and the stator is much greater than during the non-contact operation of the seal assumed in this model. The deformations may contribute to a greater leak.
It should be mentioned that the mathematical analysis based on the fractional differential equation is suitable to develop more accurate models to analyze similar physical phenomena.