Abstract
In current theory models for rotating heat pipes, the temperature field of the vapor phase is often supposed to be homogenous, and as a result of such simplification, the experiment result of the heat transfer performance for high rotating speed has some discrepancy with that predicted by theory models. In this paper, the analytical solution of the vapor flow in rotating heat pipes was obtained on the hypothesis of potential flow and with the method of variable separation. A specific rotating heat pipe was examined under three kinds of boundary conditions: linear distribution, uniform but asymmetric distribution, and uniform as well as symmetric distribution of heat load. The flow field was calculated, the Coriolis force is estimated, and it is found that (1) for a rotating heat pipe with high speed or large heat load, it is necessary to consider the Coriolis force, as its magnitude can be that of the gravitational acceleration; (2) the maximum Coriolis force is located at the vapor-liquid interface at the evaporator and condenser sections, and the directions in these two sections are opposite; (3) the Coriolis force is closely related with working conditions and working fluids, and it decreased with working temperature and increased with the heat load; and (4) the maximum viscous shear stress is located at the adiabatic section.
1. Introduction
A rotating heat pipe is a special kind of heat pipe which rotates with the object it cools, and the liquid fluid in it flows back under the centrifugal force. The structure of a rotating heat pipe includes a cylindrical rotating heat pipe [1], a parallel axis rotating heat pipe [2, 3], a radial rotating heat pipe [4], a truncated cone rotating heat pipe [1], and an even fan-shaped rotating heat pipe [5]. On account of the stable and strong driving ability caused by the centrifugal force, the advantage of heat pipes to get more even temperature field can be brought into full play by rotating heat pipes, and an important application of it is the thermal and cold protection for rotary machines, such as generators [6], rotor of gas turbines [4, 5, 7], aeroengine nose cones [8–10], milling machines [11, 12], and drills [13]. For such equipment, due to friction, to electromagnetic induction, and to being heated or even cooled [8–10], the temperature of the rotating shaft, blade, or nose cone is often uneven, causing thermal stress [4], deformation, and even failure of the machine. Thus, thermal and cold protection is required; however, there are some problems with the currently used thermal or cold protection method more or less, like economic or security ones, for instance, hydrogen cooling in generators [14] and film cooling in gas turbines. Therefore, the passive dredging thermal protection scheme of rotating heat pipes is now gradually revealing its advantage.
Current researches on rotating heat pipes contain two directions: experiment and numerical simulations, the former focuses on the relationship between thermal resistance and parameters like rotation speed and liquid charging rate, such as the research of Xie et al. [15], Ponnappan et al. [16], Song et al. [17], and so on. Ponnappan et al. [16] researched on the case of high rotation speed and found that the variation of thermal resistance was opposite to the result got by the Nusselt condensation model. While numerical simulation mainly studies on the four processes of evaporation, condensation, vapor flow, and liquid flow happened in the rotating heat pipe. According to recent researches, the thermal resistances in the four processes are of different magnitudes. The resistance in the evaporation process is usually small and can be ignored, while for the condensation process, the modified Nusselt model is widely adopted which is now deemed as sophisticated, and when it comes to vapor and liquid flow, the vapor phase is usually simplified as a system with uniform temperature and pressure, and the interior boundary condition of the liquid flow can be set based on such hypothesis. Thus, liquid flow can be solved under different types of models, according to their dimensions. For example, Daniels and Al-Jumaily [18] and Uddina et al. [19] ignored the inertia term of the liquid phase and obtained its temperature distribution by the zero dimension model. Song et al. [17], Bertossi et al. [1], and Hassan and Harmand [2, 20] considered the liquid phase as a one-dimensional flow; Li et al. [21] treated it as a two-dimensional flow. However, to discuss the performance of the heat pipe clearly under different working conditions, the effect of the vapor phase should be taken into full consideration. Faghri et al. took full researches on who resolved originally the velocity field of the vapor phase in rotating heat pipes in [22] and who got the complex distribution of vapor phase parameters in high rotating speed. In their article, the effect of the radial Reynolds number Rer on the vapor flow was intensively discussed, and it is found that when Rer is small, the circumferential velocity is proportional to the radius and the field takes on the feature of plane flow; when the Rer is large enough, the linear relationship between the circumferential velocity and the radius was broken, with the velocity becoming larger at the evaporation section and smaller at the condensation section, which indicated some circumferential forces existing and becoming obvious at high rotation speed, and it is just the Coriolis force. An Indian scholar, Solomon et al. [23], pointed out that vapor flow in the rotating heat pipe is affected by the Coriolis force to some extent, and this is the first literature to date pointing out the existence of the Coriolis force in rotating heat pipes.
In summary, it can be seen that vapor flow in a rotating heat pipe is very complicated, and its flow pattern has a significant impact on the whole performance of a rotating heat pipe; however, theoretical models aiming for rotating heat pipe simulating and mechanism revealing are insufficient. It is worth studying whether the heat transfer characteristics of a rotating heat pipe found by Ponnappan et al. at high rotation speed are affected by the vapor flow pattern, and the dynamics and thermodynamics law of vapor flow of a rotating heat pipe at high rotation speed or at high heat load is also needed to be further clarified. Thus, in this article, it is aimed to obtain the flow of vapor in rotating heat pipes by an analytical method based on the hypothesis of potential flow and especially explore the Coriolis force inflicted to the vapor phase.
2. Establishment of Physical Model for Calculation
For the sake of convenience of theoretical research and viability of the analytical method to get the solution of vapor flow, a two-dimensional rotating heat pipe model was established, which is also the mainstream model for rotating heat pipe research, and the structure of which is shown as the following figure: the rotating heat pipe has a rotating shell which was installed on the object it cools, heat was sucked on the left side—the evaporator, and released on the right side—the condenser, and the middle part of the pipe corresponds to the adiabatic section. The liquid film is attached to the inner wall of the shell, under the function of the centrifugal force. It is advisable and reasonable to suppose that the centrifugal force is large enough, so the variation of the thickness of the liquid film along the direction of the axis can be neglected. The vapor phase flow is driven by the phase change at the liquid-vapor interface, which constitutes the boundary condition of the vapor phase flow. The length of the evaporator is denoted as l1, and the length of the condenser is denoted as l2. The heat flux at the evaporator is denoted as q1, and the heat flux at the condenser is denoted as q2. Accordingly, the velocity boundary condition of the vapor flow at the evaporator is denoted as u1, and that at the condenser is denoted as u2. For the sake of convenience, the inner diameter of the heat pipe is denoted as 2b, and the radius of the liquid-vapor interface is denoted as a, so the thickness of the liquid film is b-a.
3. Mathematical Descriptions of the Model and Simplification of the Governing Equations
To solve the vapor flow in the physical model shown in Figure 1, it should be based on the NS functions coupled with the energy function and state function under rotating coordinates. In order to obtain the analytical solution, the flow of the vapor is supposed to be the potential flow; under which hypothesis, the above functions can be simplified to a great extent. The universal form of NS functions is the vector form, as shown below:
This form is regardless with the coordinates, and for the case of rotating coordinates, which is a noninertial system, the function that the relative velocity satisfies has the same form as (1), with the difference of the mass force item, which can be expressed as in which the mass force is an inertial force, the first item is the Coriolis force, and the second item is the centrifugal force. Between the two forces, the former is a nonpotential force and the direction of it is circumferential, while the latter is a potential force and the direction of it is radial. In most literatures, the Coriolis force is often neglected, and only the centrifugal force is taken into account, such treatment can make the model simple and, to some extent, can result in with acceptable accuracy. So firstly, the Coriolis force is also neglected by former researchers, and the mathematical model can be simplified accordingly, and then the estimation of the Coriolis force is carried out based on the obtained solution. For the potential hypothesis, the viscous item in (2) can be neglected, and therefore, the governing equations can be reduced into
In summary, the hypothesis applied for this problem is listed below explicitly: (1)The flow is a two-dimensional plane flow, and the circumferential velocity is neglected.(2)The viscous items concerning the parameter viscosity are all neglected, so the flow is the potential flow.(3)The vapor phase is incompressible, and its state parameters are all saturated ones.(4)Estimation of the Coriolis force is based on the potential flow, but in the solving process, the Coriolis force is neglected, so that their coupling relationship is simplified.
So supposing the potential function of the velocity is denoted as ψ, then it should satisfy the Laplace function, as shown below:
4. Solving the Velocity Field of the Vapor Phase
The calculating region of the mathematical problem above is written as
The left boundary condition of the vapor flow is
And the right boundary condition is
The boundary condition at the symmetry axis is
The boundary condition at the vapor-liquid interface is
The notations z1 and z2 are defined as
The governing equation (4) can be denoted below in rotating coordinates:
For the two-dimensional plane flow problem discussed in this paper, the second item in (11) can be neglected, and the governing equation can be turned as
This equation can be solved by the variable separation method, and its analytical solution can be expressed in (13); the deduction process of which is illustrated in the appendix.
It is taken as the form of a series, in which βm is the eigenvalues and can be calculated as below:
Am is the coefficient of each item in the series. I0(βmr) is the zero-order modified Bessel function. The partial differentials of the function ψ to r and z are the two velocity solutions, ur and uz, which are shown below:
The coefficient Am is determined by the boundary condition (9), and for different boundary conditions, they have different values. The derivation process is illustrated in the appendix, and the results are shown below.
4.1. Linear Distribution of the Boundary Velocity
When the distribution of heat flux is linear, the boundary condition of the radial velocity of the vapor flow is linear, and it can be described as
And the solution of the coefficients Am is shown below:
4.2. Uniform Distribution of the Boundary Radial Velocity
When the distribution of heat flux is uniform at the evaporation and condensation sections, the boundary condition of the radial velocity of the vapor flow is also uniform, and it can be described below:
And the solution of the coefficients Am is shown below:
4.3. Uniform and Symmetric Distribution of the Boundary Radial Velocity
One special case of the uniform distribution of boundary condition is the symmetric case, in which the length of the evaporator and that of the condenser are equal and the magnitude of radial velocity U1 and that of radial velocity U2 are equal too, as shown below:
For this case, the solution of the coefficients is shown below, which means that the even terms are all zero.
So the final solution of the velocity potential can be expressed as
And the velocities ur and uz can be expressed as
5. Results and Discussion
5.1. The Potential Solution of the Flow Field for Three Cases of Boundary Condition
Directed by the analytical solution obtained above, a specific heat pipe was calculated; the structure parameters of which were set as follows: the length L is 200 mm, the inner diameter of the heap pipe is 40 mm, and the thickness of the liquid film is 4 mm. In order to compare the influence of boundary condition on the distribution of vapor flow, especially the Coriolis force, three cases for different heat loads are calculated which are (1) uniform as well as symmetric distribution, (2) linear distribution, and (3) uniform but asymmetric distribution of heat load.
The calculated flow field for case (1) is drawn as Figure 2, in which the velocity field of the liquid film is also exhibited, which is obtained with the same method as the vapor flow. In this case, the length of the evaporator and that of the condenser are both 40 mm. In the figure, the black lines denote the streamline; the vectors denote the direction and the magnitude of the velocity at each point, while the colored lines stand for the contour line of the velocity potential. In the figure, subfigure (a) corresponds to the velocity field of the liquid film at the evaporator; subfigure (b) corresponds to the velocity field of the liquid film at the condenser section; subfigure (c) corresponds to the velocity field of the vapor phase at the evaporator; and subfigure (d) corresponds to the velocity fields of the vapor phase at the condenser section.
(a)
(b)
(c)
(d)
Figure 3 exhibits the wholesome vapor flow field of the rotating heat pipe for the boundary condition as case 2, in which the heat load at the evaporator and the condenser is linearly distributed. The length of the evaporator is 40 mm, and along the axial direction of z, the heat source at the evaporator increases linearly; the length of the condenser is 120 mm, and along the axial direction of z, the heat sink at the condenser decreases linearly.
Figure 4 exhibits the wholesome vapor flow field of the rotating heat pipe for the boundary condition as case 3, in which the heat load at the evaporator and the condenser is uniformly distributed. The length of the evaporator is 40 mm, and the length of the condenser is 120 mm. As the length of the evaporator and that of the condenser are not identical, therefore, the magnitudes of the heat load as well as the boundary velocity at both sides are not equal.
5.2. Discussion of the Coriolis Force
Once the solution of the velocity field is obtained, the Coriolis force can be discussed, which is closely related with the radial velocity, and is calculated by the following:
In the equation, fc denotes the Coriolis force per unit mass, so its dimension is that of acceleration, m/s2, ur denotes the radial velocity of the working fluid, and ω denotes the absolute angular velocity relative to the ground; as the flow is supposed to be a plane flow, it is equal to the angular velocity of the shell. It can be deducted that the Coriolis force is proportional to the rotating speed of the rotating heat pipe, but it is also noticeable that when the circumferential velocity is taken into account, the angular velocity of the working fluid and that of the shell are not equal, which is just the result of the Coriolis force.
In this article, a specific working condition is calculated, in which the rotating speed n is 3000 rpm, the structure parameters are the same as that discussed above, the working fluid is R134a, which is a common refrigerant, and the working temperature of the vapor phase T is 50°C, so its calculating parameters are presented in Table 1.
The heat load conditions are also set for three cases as before; for the uniform distribution cases, the heat source is set to be 5 W/cm2, and for the linear distribution case, its total heat load is set to be equal to the other two cases. The distribution of the Coriolis force along the axial direction is exhibited in Figure 5. In this figure, subfigure (a) exhibits the distribution of the Coriolis force in the vapor phase in the case with linear distribution of heat load; subfigure (b) exhibits the distribution of the Coriolis force in the vapor phase in the case with uniform but asymmetric distribution of heat load; and subfigures (c) and (d) exhibit the distribution of the Coriolis force in the vapor phase and the liquid film in the case with uniform and symmetric distribution of heat load.
(a)
(b)
(c)
(d)
It can be seen that the direction of the Coriolis force at the evaporator and that at the condenser are different, the Coriolis force at the adiabatic section is negligible, and as the radial coordinates increase, the magnitude of the Coriolis force increases, so the maximum value of the Coriolis force is located at the interface of the liquid-vapor phase. And it can also be observed that the magnitude of the maximum Coriolis is that of the gravitational acceleration g. As the rotating velocity 3000 rpm is a typical value, such as the dynamos and steam turbines in power plants and the heat load is not an extreme value, it is reasonable to consider the Coriolis force in rotating heat pipes. From subfigures (c) and (d), it can also be concluded that the Coriolis force in the vapor phase is much larger than that in the liquid phase and, in the liquid phase, the Coriolis force can be negligible.
Figure 6 exhibits the distribution of the Coriolis force along the radial direction. In this figure, subfigure (a) exhibits the distribution of the Coriolis force in the vapor phase in the case with a linear distribution of heat load; subfigure (b) exhibits the distribution of the Coriolis force in the vapor phase in the case with uniform but asymmetric distribution of heat load.
(a)
(b)
In this figure, z = 30 mm is located at the evaporator section, z = 50 mm is located at the adiabatic section, z = 160 mm is located at the condenser section, and z = 45 mm is located at the boundary of the evaporator section and the adiabatic section. It can be seen that, at the most region of the evaporator and the condenser sections, the relationship between the Coriolis force and the radial coordinate is nearly linear, and the more uniform the boundary condition is, the stronger the linearity is. But it is also noticeable that this solution is based on a two-dimensional plane flow; if the circumferential velocity is taken into consideration, further investigation into the degree of this linearity is needed.
The magnitude of the Coriolis force is influenced by working conditions such as the heat load q and the working temperature of the vapor phase T. What is more, it has a great relationship with the working fluid itself. Figure 7 exhibits the maximum Coriolis forces at different working conditions for two specific working fluids: R134a (a) and water (b), and boundary condition distribution for both working fluids is as case (1): uniform and symmetric.
(a)
(b)
From Figure 7, it can be seen that the magnitude of the Coriolis force decreases with working temperature and increases with the heat load. It can also be found that the Coriolis force is larger when the working fluid is water, compared with R134a. As the heat loads calculated are possible ones, it can be concluded that it is possible for the Coriolis force to be 1~10 times the gravitational acceleration.
5.3. Viscous Shear Stress at the Interface
From the obtained axial velocity, the shear stress can be discussed. Because of the potential flow, hypothesis neglects the viscous items in the governing equation, so at the interface of the liquid and vapor phases, the calculated velocity cannot be continuous, as shown in Figure 8.
(a)
(b)
From this figure, the viscous shear stress can be evaluated by the axial velocity difference of the vapor and liquid phases across the interface. Subfigure (a) corresponds to the case of linear distribution boundary condition, and subfigure (b) corresponds to the case of the uniform but asymmetric distribution boundary condition. The total heat loads for the two cases are equal, and the heat flux of the evaporator for uniform distribution case is 5 W/cm2. For the two cases, the length of the evaporator and that of the condenser are identical, which means that z = 30 mm is located at the evaporator section, z = 45 mm is located at the adiabatic section, and z = 140 mm and z = 160 mm are located at the condenser section. It can be seen that the axial velocity of the liquid phase is approximately zero and the maximum axial velocity difference is located at the adiabatic section. So the viscous shear stress at this region needs to be further investigated based on a more accurate model, in which the viscous items should be taken into account.
It can be deducted that if the nonslip boundary condition is considered, the gradient of the axial velocity along the axial direction will shift a lot, so according to the continuity equation, which is written below, the gradient of the radial velocity along the radial direction will shift a lot, and the linearity relationship between the Coriolis force and the radial direction indicated by the potential flow solution will change, and a more exact relationship needs to be explored.
6. Conclusions and Outlooks
The analytical solution of the vapor flow in rotating heat pipes was obtained in the hypothesis of potential flow and with the method of variable separation, which is a solution of a series form. A specific rotating heat pipe was examined under the three kinds of boundary conditions: linear distribution, uniform but asymmetric distribution, and uniform as well as symmetric distribution of heat load. The flow field was calculated, and the Coriolis force is estimated for different boundary conditions, different working fluids, and different working conditions, and the following conclusions were arrived: (1)When neglecting the Coriolis force and the viscous item in the governing equations, it is viable to obtain the analytical solution of the working fluid flow, especially for the vapor phase, with the method of variable separation, and such method is applicable for any cases of heat load distribution.(2)The magnitude of the Coriolis force can be that of the gravitational acceleration; for rotating heat pipes with high speed and large heat load, the Coriolis force needs to be considered.(3)Along the axial direction, the maximum Coriolis force is located at the vapor-liquid interface at the evaporator and condenser sections, and the directions of the Coriolis force at these two sections are opposite, while at the adiabatic section, the Coriolis force is negligible.(4)The Coriolis force is closely related with working conditions and working fluids, and it decreased with working temperature and increased with the heat load.(5)The influence range of the viscous force can be estimated according to the axial velocity difference of the vapor phase and the liquid phase across the interface, and the maximum viscous shear stress is located at the adiabatic section.
6.1. Outlooks
As the Coriolis force can result in the circumferential velocity, it can be inferred that the flow can take on circumferential shift, both in the evaporator section and in the condenser section. As the directions of the Coriolis force are opposite in the two sections, the circumferential velocity is also opposite. So the mutual interaction of the distribution of the circumferential velocity and the Coriolis force as well as the viscous shear stress are needed to be further clarified.
Appendix
A. Deduction of the Series Form Solution
A.1. General Solution of the Velocity Potential Solution
For the governing equation of the velocity potential function (12), it is supposed that its solution can be expressed as
Substituting it into function (12), it is obtained that
Dividing it with equation ψ = ZR, it is obtained that
If this equation is always satisfied regardless of any variables, then it can be supposed that
And for the function Z, the following function should be satisfied:
And the boundary conditions are shown below:
The solution to the above problem is a solution system, which can all be written as the form in which the eigenvalue βm is
While for function R, it satisfies the function below:
This function is a zero-order modified Bessel equation, and its boundary conditions are shown below:
According to boundary condition (A-10), it is convenient to express the solution of R as shown below: where the notation I0 stands for the first kind zero-order modified Bessel function. And thus, the solution of the problem discussed in this article can be expressed as a series: where the notation Am is a series of coefficient to be determined.
A.2. Settlement of the Coefficients Am for a Given Boundary Condition
Once the expression of the velocity potential has been obtained, the velocity field can be deducted by the partial difference of ψ to the two coordinate variables, r and z, and the expression of the radial velocity ur and the axial velocity uz can be written as follows:
And the coefficient Am can be determined based on radial velocity expression (A.14) by the orthogonal integral method, which takes into account boundary condition (A.11) and multiplies function (A.7) to the two sides of (A.14) and then makes integral at the boundary r = a, so the left side became
The symbol Jleft stands for the integral above, and by integration by parts, it can be transformed into the following form:
The specific expression is determined by the specific expression u1 and u2. For the linear distribution of boundary condition (A.11), as shown in (16), the expression of (A.17) can be obtained as shown below:
For the uniform distribution of boundary condition (A.11), as shown in (18), expression (A.17) can be written as
For the uniform and symmetric distribution of boundary condition (A.11), as shown in (20), expression (A.17) can be written as
The orthogonal integral to the right side of (A.14) can be expressed as
And its result is shown below, which is universal and regardless of the specific form of boundary condition (A.11).
The two sides of the orthogonal integral are equal, so the coefficient Am can be obtained, as shown below:
In the expression above, the derivative of the zero-order modified Bessel function I0′ is substituted by I1, which is the first-order modified Bessel function and is just the derivative of the zero-order modified Bessel function. The above equation is applicable for cases when m > 0, while for m = 0, the first term in (A.13) is
This item is a constant but cannot be settled by the method of orthogonal integral to the partial difference of the potential function, so its value is unknown. But for the problem we care, this constant has no influence on the settlement of the velocity field, because it can be eliminated by derivation. So for the sake of simplicity, its value is given as 0. So the solution of velocity potential (A.13) can be rewritten as
Nomenclature
| Denotation/unit: | Physical meaning |
| a/mm: | Radius of the vapor-liquid interface |
| Am/m2/s: | Coefficients of each item in the solution series |
| b/mm: | Inner radius of the shell |
| fc/m/s2: | Coriolis force per unit mass |
| hfg/kJ/kg: | Latent heat of evaporization |
| I0: | The zero-order modified Bessel function |
| I1: | The first-order modified Bessel function |
| k1/s−1: | Coefficient of linear velocity boundary condition at the evaporator |
| k2/s−1: | Coefficient of linear velocity boundary condition at the condenser |
| L/mm: | Length of the rotating heat pipe |
| l1/mm: | Length of the evaporator |
| l2/mm: | Length of the condenser |
| m: | Serial number of each item in the solution series |
| n/rpm: | Rotating speed of the heat pipe |
| p/Pa: | Pressure of the working fluid |
| q1/W/cm2: | Heat load at the evaporator |
| q2/W/cm2: | Heat load at the condenser |
| r/mm: | Radial coordinate |
| R: | Separated function about r |
| T/K(°C): | Working temperature of the vapor phase |
| u1/m/s: | Velocity boundary condition at the evaporator |
| U1/m/s: | Uniformly distributed velocity boundary condition at the evaporator |
| u2/m/s: | Velocity boundary condition at the condenser |
| U2/m/s: | Uniformly distributed velocity boundary condition at the evaporator |
| ur/m/s: | Radial velocity solution |
| uz/m/s: | Axial velocity solution |
| z/mm: | Axial coordinate |
| Z: | Separated function about z |
| βm: | Eigenvalues of the solution |
| ψ/m2/s: | Velocity potential |
| ω/s−1: | Angular velocity |
| ρv/kg/m3: | Density of the vapor phase. |
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgment
This project is supported by the Fundamental Research Funds for the Central Universities, China (Grant no. 2014QNB07) and Basic Research Project of Natural Science Foundation of Jiangsu Province (BK20170280).