A Multi-Node Lumped Parameter Model Including Gravity and Real Gas Effects for Steady and Transient Analysis of Heat Pipes

Abstract

This work presents a multi-node lumped parameter model able to predict the steady and transient behavior of capillary heat pipes, taking into account the effects of gravity (orientation angle) and the real gas effects in the vapor modeling. The model was validated against experimental results acquired by Leonardo S.p.A., which were obtained by simulating the behavior of a heat pipe embedded in a chassis cover, subject to seven cycles of transient thermal loading. After the validation, the analysis is focused on the model accuracy when using the ideal and real gas assumptions, using different working fluids (water, ammonia, acetone, HFC134a). The results showed that when using water or ammonia as working fluid, the error in modeling the vapor as an ideal instead of as real gas is negligible, both for the vapor temperatures and pressures predictions. On the contrary, when using acetone or HFC134a as working fluid, modeling the vapor as a real gas leads to a significant increase in the accuracy of the vapor pressure predictions.

1. Introduction and State of the Art

Heat pipes (HP) are very attractive components in the area of spacecraft cooling and temperature stabilization due to their low weight penalty, very low maintenance, and high reliability [1]. Nowadays, heat pipes have been studied for more than fifty years, and they are one of the most popular passive thermal control devices for many engineering fields.

The structure of a heat pipe consists of a sealed system containing a fluid, the so-called working fluid, which cyclically evaporates and condenses. Through the evaporation process, the working fluid removes heat from the component to cool down. This thermal power removed is dissipated into the external environment by means of the condensation process. Therefore, the heat pipe is a device with very high thermal conductivity that enables the transportation of heat whilst maintaining almost uniform temperature along its heated and cooled sections [2]. This means that it works in nearly isothermal conditions.

The capillary driven heat pipe (CHP) is the most commonly used and known type of heat pipe. In fact, it is used in many industrial applications, such as computer cooling [3] and spacecraft thermal control [1]. This type of heat pipe consists of a sealed container, whose inner surface is covered with a porous structure, the wick. The container is partially filled with a working fluid [4]. One end of the container, the evaporator, is heated by a heat source. Thus, the working fluid evaporates, and the vapor moves to the other end of the HP, the condenser, possibly passing through an adiabatic zone. The condenser area is cooled so that the vapor can condense. The function of the wick material on the inner surface of the heat pipe is to guide the condensate back to the evaporator [4]. This process occurs mainly because of the action of capillary forces. However, gravitational force can have a great influence on the behavior of the CHP. In fact, its presence can make the HP operation orientation-dependent, and only a proper design of the wick can allow the HP to operate in all directions and to reach the requested performance [3]. The material composing the wick must be compatible with the working fluid, in order to avoid undesired effects, such as corrosion or generation of non-condensable gas. For this reason, the most common wick/fluid couplings are: copper/water, mainly used for electronics applications; aluminum/ammonia, and aluminum/acetone, mainly used for spacecraft applications; copper/HFC134a refrigerant, mainly used for energy recovery applications.

In the past decades, there have been many attempts to model the different types of heat pipes.

A possible way to model both the steady-state and transient behavior of a loop heat pipe was described by Cullimore and Baumann [5]. This is done by using the NASA-standard thermohydraulic analyzer, SINDA/FLUINT, to design and simulate thermal/fluid systems that can be represented in networks corresponding to finite difference, finite element, and/or lumped parameter equations. The only disadvantage of this model is that it requires a good knowledge of loop heat pipe operation and of the thermophysical processes involved in those devices, at least for detailed studies.

A transient model for microgrooved heat pipes by using a macroscopic approach was presented by Suman et al. [6]. This model was developed by using an equilateral triangular heat pipe as test case, but the equations are suitable to be used with heat pipes of any polygonal shape. The flow is modeled as unsteady, incompressible, and one-dimensional along the length of the heat pipe. The coupled non-linear governing equations of heat, mass, and momentum transfer are simultaneously solved to obtain the transient and steady-state profiles of the main parameters which characterize the heat pipe behavior (temperature, pressure, flow velocity). However, this model can be used only for microgrooved heat pipes with a polygonal shape.

A lumped parameter model for capillary heat pipes was developed by Ferrandi et al. [7,8]. This model, based on the electrical analogy, divides the heat pipe in two parts, the solid network and the fluidic network, that are connected by coupling equations. The model considers the vapor as an ideal gas and it does not take into account the influence of gravity (and therefore of the orientation angle) on the heat pipe behavior.

Another lumped parameter model for loop heat pipes was developed by Bernagozzi et al. [9]. As the previous one, this model is based on the electrical analogy and it does not consider the real gas effects in the vapor modeling and the gravity effects.

A model able to predict the heat transfer performances of a closed vertical meandered pulsating heat pipe, filled with water, was proposed by Wang et al. [10], and it is based on artificial neural networks. This model can be used to predict the thermal resistance of the pulsating heat pipe, given filling ratio and related geometric parameters, at different heat fluxes. The analysis showed a good agreement between numerical and experimental results. However, in order to use artificial neural networks for predicting the pulsating heat pipe behavior, some other studies with different working fluids and more experimental data are necessary.

A simple model able to predict the thermal performance and working temperature of a multi-channel flat plate heat pipe was developed by Guichet et al. [11], and it includes the cooling manifold to estimate the characteristics of the complete system. However, this model was developed only for a very specific type of heat pipe, so it could be difficult to use in different applications.

An advanced conduction-based heat pipe model was recently developed by Zimmermann et al. [12]. An analytic expression to calculate the pressure drop inside the vapor core was developed, and the results were validated through CFD simulations. This model is quite easy to implement and it showed good results, but it is designed for steady conditions. Thus, there is no evidence that it also works in transient conditions.

The aim of this work is to build and validate a simple, multinode, lumped parameter model, which is able to predict the steady and transient behavior of capillary heat pipes, including the gravity (orientation angle) effects and the real gas effects in the vapor modeling. The model was built to be quite versatile, in fact, it can be used for heat pipes of every possible shape and geometry that can be described and/or modeled analytically, as well as for all the solid–fluid couples whose thermophysical properties are known. The model was validated by using a real dataset provided by the company Leonardo S.p.A., referring to an experiment performed by using a flat heat pipe made of copper and filled with water as working fluid. After the validation, the model was used for predicting the heat pipe behavior, when four different solid–fluid couples are used (copper–water, aluminum–ammonia, aluminum–acetone, copper–HFC134a refrigerant). A specific focus is placed on assessing for which fluids and working conditions the ideal gas assumption grants a good accuracy, and for which ones, on the contrary, its use is not suitable, and real gas properties must be included in the model.

2. Heat Pipe Model Description

Figure 1 shows the conceptual model of the heat pipe used in this work. To use the lumped parameter approach, the heat pipe scheme can be conceptually divided into two major parts, the solid network and the fluid network, as in [7,8,9]. Each of these networks can be further divided into different minor parts, called nodes or lumps. Each node represents a section of the complete heat pipe, and all its thermophysical properties are considered uniform across the whole section, and equal to their average value in the section, which is assumed to be the one at the center of the volume of the node which represents the considered section (midpoint rule).

The solid network is modeled through the use of the electrical analogy, in which the currents are represented by the heat transfer rates, and the voltages are represented by the temperature differences between two nodes. All the nodes are connected to the adjacent ones, to the fluid network, and to the external environment through thermal equivalent resistors and capacitors.

Since the fluid network must also take into account the presence of gravitational terms, it cannot be modeled by the electrical analogy, because it is not a passive network as it would have been without gravitational forces. Therefore, this part of the system is modeled through the one-dimensional Navier–Stokes equations and the equations of real gases.

Then, the coupling between the solid and fluid networks, and the one between the liquid and the vapor must be considered in order to obtain the complete set of equations describing the behavior of the entire heat pipe.

2.1. Solid Network

The scheme of the solid network is represented in Figure 1, where dashed lines separate the vapor zone, wick, solid wall, and external regions (evaporator and condenser).

The solid network is modeled using a total of six nodes, three in the solid wall ($T_{PE}$, $T_{PA}$, $T_{PC}$), and three in the wick ($T_{WE}$, $T_{WA}$, $T_{WC}$). For both the wall and the wick, there is a node representing the evaporator section (subscript E), a node representing the adiabatic section (subscript A), and a node representing the condenser section (subscript C). These nodes are connected to three nodes in the vapor zone ($T_{VE}$, $T_{VA}$, $T_{VC}$), which belongs to the fluid network, and to two nodes on the external walls, namely, $T_E$ at the evaporator, $T_C$ at the condenser. The latter is in its turn connected to a node representing the external environment ($T_F$).

In Figure 1, the symbol T represents the temperature of a node, the symbol $\dot{Q}$ represents a heat power, which can be the one entering the system from the evaporator wall, ${\dot{Q}}_{IN}$, or the one exiting the system from the condenser wall, ${\dot{Q}}_{OUT}$. The symbol R represents the resistance of a thermal equivalent resistor, and the symbol C represents the capacitance of a thermal equivalent condenser. The resistances in the solid wall and in the wick are conductive, the ones connecting the wick to the vapor zone are convective, and the one connecting the condenser wall to the external environment is mixed, which means that it can be conductive, convective, radiative, or a combination, depending on the boundary conditions of the specific case.

The resistances $R_{1PA}$, $R_{2PA}$, $R_{1WA}$, and $R_{2WA}$ are axial, in fact, they model what happens along the longitudinal dimension of the heat pipe, namely, from the evaporator to the condenser. Conversely, all the other resistances are radial, describing the thermal behavior along the radial direction of the heat pipe, namely, from the center of the vapor zone to the solid wall.

The resistance and capacitance of each branch of the network depend on the geometry of the heat pipe section. Therefore, the equations used to describe them will be analyzed in Section 3, where the geometry of the case study is described.

At this point, writing an energy balance on each solid node, it is possible to draw six equilibrium equations.

For the evaporator wall node, the equilibrium equation is:

$$C_{PE}·\frac{d{T}_{PE}}{d\tau }={\dot{Q}}_{IN}+\frac{T_{PA}-T_{PE}}{R_{1PA}}+\frac{T_{WE}-T_{PE}}{R_{2PE}+R_{1WE}}$$

(1)

where $\tau$ is time. In Equation (1), either the input thermal power, ${\dot{Q}}_{IN}$, or the evaporator wall temperature, $T_E$, has to be known. In this second case, the input thermal power to be used in Equation (1) is:

$${\dot{Q}}_{IN}=\frac{T_E-T_{PE}}{R_{1PE}}$$

(2)

For the condenser wall node, the equilibrium equation is:

$$C_{PC}·\frac{d{T}_{PC}}{d\tau }=-{\dot{Q}}_{OUT}+\frac{T_{PA}-T_{PC}}{R_{2PA}}+\frac{T_{WC}-T_{PC}}{R_{2PC}+R_{1WC}}$$

(3)

In Equation (3), either the output thermal power, ${\dot{Q}}_{OUT}$, or the external environment temperature, $T_F$, has to be known. As before, in this second case, the output thermal power to be used in Equation (3) is:

$${\dot{Q}}_{OUT}=\frac{T_{PC}-T_F}{R_{1PC}+R_M}$$

(4)

For the adiabatic wall node, the equilibrium equation is:

$$C_{PA}·\frac{d{T}_{PA}}{d\tau }=\frac{T_{PE}-T_{PA}}{R_{1PA}}+\frac{T_{PC}-T_{PA}}{R_{2PA}}+\frac{T_{WA}-T_{PA}}{R_{PA}+R_{WA}}$$

(5)

For the three wick nodes (evaporator, condenser, and adiabatic zone), the equilibrium equations are respectively:

$$C_{WE}·\frac{d{T}_{WE}}{d\tau }=\frac{T_{PE}-T_{WE}}{R_{1WE}+R_{2PE}}+\frac{T_{WA}-T_{WE}}{R_{1WA}}+\frac{T_{VE}-T_{WE}}{R_{2WE}+R_{EF}}$$

(6)

$$C_{WC}·\frac{d{T}_{WC}}{d\tau }=\frac{T_{PC}-T_{WC}}{R_{1WC}+R_{2PC}}+\frac{T_{WA}-T_{WC}}{R_{2WA}}+\frac{T_{VC}-T_{WC}}{R_{2WC}+R_{CF}}$$

(7)

$$C_{WA}·\frac{d{T}_{WA}}{d\tau }=\frac{T_{PA}-T_{WA}}{R_{WA}+R_{PA}}+\frac{T_{WE}-T_{WA}}{R_{1WA}}+\frac{T_{WC}-T_{WA}}{R_{2WA}}+\frac{T_{VA}-T_{WA}}{R_{3WA}+R_{AF}}$$

(8)

The vapor temperature in the evaporator zone, $T_{VE}$, in Equation (6), and the vapor temperature in the condenser zone, $T_{VC}$, in Equation (7), are taken from the fluid network. Moreover, in Equation (8), the vapor temperature in the adiabatic zone, $T_{VA}$, has to be known. It is computed as a weighted average of temperatures $T_{VE}$ and $T_{VC}$, as shown in System (9):

$$\left\{\begin{array}{c}{T}_{VA}=w_E·T_{VE}+w_C·T_{VC}\hfill \\ w_E=\frac{\frac{L_C}{2}+\frac{L_A}{2}}{L_{eff}}\hfill \\ w_C=\frac{\frac{L_E}{2}+\frac{L_A}{2}}{L_{eff}}\hfill \end{array}\right.$$

(9)

where $L_E$, $L_A$, and $L_C$ are the evaporator, adiabatic zone, and condenser lengths, respectively; $L_{eff}$ is the effective duct length ($L_{eff}=L_E/2+L_A+L_C/2$), $w_E$ and $w_C$ are the weights that take into account the fact that the center of the adiabatic zone can be closer to the condenser than to the evaporator, or vice versa.

Therefore, Equations (1), (3) and (5)–(8) are the six equilibrium equations which describe the temperatures evolution in the solid network.

2.2. Fluid Network

In order to write the equations which describe the fluid network, six assumptions are necessary. It is assumed that:

• the fluid flow is one-dimensional and laminar for both phases;
• vapor transformations in the evaporator/condenser tanks are isoentropic;
• evaporation and condensation only occur in the evaporator and condenser zones, therefore, the mass flow rates of the phases do not change in the adiabatic region;
• the liquid temperature is exactly equal to the wick temperature for each point along the wick length;
• the thermophysical properties of the liquid, the surface tension, and the latent heat of vaporization are evaluated in saturation conditions;
• All the condensing fluid reaches the wick, with no “accumulation” before it.

With reference to Figure 1, heat enters the fluid network arriving from the external heat source into the evaporator liquid tank, so that the working fluid in liquid state evaporates and conceptually flows into the vapor tank, thanks to the capillary pressure jump. Then, the vapor formed flows through the adiabatic vapor duct and arrives to the condenser vapor tank, from which a certain quantity of heat is removed. This phenomenon causes the vapor to condense and flow back to the evaporator, where the cycle restarts.

At this point, it is possible to analyze each piece of the fluid network, and to write the equations which describe its behavior. In such equations, the mass flow rate of the vapor is indicated with ${\dot{m}}_V$, the one of the liquid is indicated with ${\dot{m}}_L$, the one that evaporates is indicated with ${\dot{m}}_E$, and the one that condenses is indicated with ${\dot{m}}_C$.

2.2.1. Condensed Liquid

For the last assumption in Section 2.2, it is possible to state that:

$${\dot{m}}_L={\dot{m}}_C$$

(10)

Equation (10) highlights the fact that there is no accumulation of liquid in the condenser.

It represents the first of the nine equations for the fluid network.

2.2.2. Vapor Tanks

In order to write the equations that link the vapor pressures and mass flow rates, it is necessary to start from the continuity equation. For the vapor, it can be written as:

$$V_V·\frac{d{\rho }_V}{d\tau }=\sum _{i=1}^N{\dot{m}}_{Vi}$$

(11)

where $V_V$ is the volume of the vapor tank considered, ${\rho }_V$ is the vapor density, ${\dot{m}}_{Vi}$ are the N vapor mass flow rates crossing the tank boundaries (taken as positive when entering). However, it is possible to write the vapor density derivative with respect to time as a function of the vapor pressure $p_V$, as:

$$\frac{d{\rho }_V}{d\tau }=\frac{\partial {\rho }_V}{\partial p_V}\frac{d{p}_V}{d\tau }$$

(12)

For the second assumption in Section 2.2, the vapor in the tanks undergoes isoentropic expansion or compression. In isoentropic conditions, it is possible to state that:

$$\frac{\partial p_V}{\partial {\rho }_V}=a_V^2$$

(13)

where $a_V$ is the speed of sound in the vapor at the current temperature and pressure. Therefore, substituting Equation (13) in Equation (12), we obtain the two equations which link the vapor pressures and mass flow rates, from the point of view of the evaporator and condenser tanks, respectively:

$$\left\{\begin{array}{c}\frac{V_{VE}}{a_{VE}^2}·\frac{d{p}_{VE}}{d\tau }={\dot{m}}_E-{\dot{m}}_V\hfill \\ \frac{V_{VC}}{a_{VC}^2}·\frac{d{p}_{VC}}{d\tau }={\dot{m}}_V-{\dot{m}}_C\hfill \end{array}\right.$$

(14)

where $p_{VE}$ and $p_{VC}$ are the vapor pressures in the evaporator and condenser zones, respectively; $V_{VE}$ and $V_{VC}$ are the volumes of the evaporator vapor tank and condenser vapor tank, respectively; $a_{VE}$ and $a_{VC}$ are the sound speeds in the vapor in the evaporator and condenser zones, respectively.

The two equations in System (14) must be added to the complete set of equations for the fluid network.

2.2.3. Vapor Pressure–Temperature Relation

The general relation between the vapor temperature and pressure is not trivial. However, since the considered vapor transformations are isoentropic, it is possible to write the entropy differential, $ds$, as a function of T and p, and to impose this differential equal to zero. In particular [13]:

$$ds={\left(\frac{\partial s}{\partial T}\right)}_{p}dT+{\left(\frac{\partial s}{\partial p}\right)}_{T}dp=\frac{c_p}{T}dT-\frac{k_p}{\rho }dp$$

(15)

where $c_p$ is the specific heat at constant pressure of the considered gas, and $k_p$ is the isobaric expansion coefficient. Imposing the entropy differential equal to zero, dividing Equation (15) by the time differential, $d\tau$, and referring the terms of Equation (15) to the vapor in the evaporator and condenser, one obtains:

$$\left\{\begin{array}{c}\frac{d{T}_{VE}}{d\tau }=T_{VE}\left(\frac{k_{pVE}}{{\rho }_{VE}·c_{pVE}}\right)\frac{d{p}_{VE}}{d\tau }\hfill \\ \frac{d{T}_{VC}}{d\tau }=T_{VC}\left(\frac{k_{pVC}}{{\rho }_{VC}·c_{pVC}}\right)\frac{d{p}_{VC}}{d\tau }\hfill \end{array}\right.$$

(16)

where $k_{pVE}$ and $k_{pVC}$ are the vapor isobaric expansion coefficients in the evaporator and condenser zones, respectively; ${\rho }_{VE}$ and ${\rho }_{VC}$ are the vapor densities in the evaporator and condenser zones, respectively; $c_{pVE}$ and $c_{pVC}$ are the vapor specific heats at constant pressure in the evaporator and condenser zones, respectively.

The two equations of System (16) are two of the nine equations that must be added to the complete set for the fluid network.

2.2.4. Evaporation and Condensation

Since the evaporation and condensation phenomena are lumped in the evaporator and condenser zones, for the third assumption in Section 2.2, the mass flow rates which evaporate and condense can be computed as:

$$\left\{\begin{array}{c}{\dot{m}}_E=\frac{{\dot{Q}}_E}{h_{LVE}}\hfill \\ {\dot{m}}_C=\frac{{\dot{Q}}_C}{h_{LVC}}\hfill \end{array}\right.$$

(17)

where ${\dot{Q}}_E$ is the thermal power going from node $T_{WE}$ to node $T_{VE}$, namely, the effective thermal power used for evaporation; ${\dot{Q}}_C$ is the thermal power going from node $T_{VC}$ to node $T_{WC}$, namely, the effective thermal power extracted for condensation; $h_{LVE}$ is the latent heat of vaporization computed at an intermediate temperature between $T_{WE}$ and $T_{VE}$; and $h_{LVC}$ is the latent heat of vaporization computed at an intermediate temperature between $T_{WC}$ and $T_{VC}$.

The equations of System (17) have to be added to the complete set of equations for the fluid network.

2.2.5. Vapor Line

Considering the vapor line, it is possible to write the one-dimensional and laminar form of the Navier–Stokes equations describing the conservation of momentum, as:

$$\rho ·\frac{\partial u}{\partial \tau }-\mu ·\frac{{\partial }^{2}u}{\partial x^2}+\frac{\partial p}{\partial x}+\rho gsin\left(\alpha \right)=0$$

(18)

where u is the vapor velocity, g is the gravity acceleration, $\mu$ is the vapor dynamic viscosity, and $\alpha$ is the heat pipe orientation angle. In particular, $\alpha$ can range between −90° and 90°: it is equal to 0° when the heat pipe is oriented horizontally; it is negative when the evaporator is above the condenser; it is positive when the evaporator is under the condenser. In most cases, this last condition is the most favorable for heat pipe operation because gravity helps the liquid to flow from the condenser to the evaporator. Conversely, when $\alpha$ is negative, the liquid pressure in the condenser must be high enough to compensate the gravity pushing the liquid back to the condenser zone.

Equation (18) can be integrated over the effective duct length, $L_{eff}$, as:

$${\int }_{L_{eff}}\left(\rho ·\frac{\partial u}{\partial \tau }\right)dl-{\int }_{L_{eff}}\left(\mu ·\frac{{\partial }^{2}u}{\partial x^2}\right)dl+{\int }_{L_{eff}}\left(\frac{\partial p}{\partial x}\right)dl+{\int }_{L_{eff}}\left(\rho gsin\left(\alpha \right)\right)dl=0$$

(19)

Starting with the time derivative term, it is possible to write (using the relation ${\dot{m}}_V=\rho u{A}_V$ in the last step):

$${\int }_{L_{eff}}\left(\rho ·\frac{\partial u}{\partial \tau }\right)dl={\rho }_{VA}·L_{eff}·\frac{d{u}_V}{dt}=\frac{L_{eff}}{A_V}·\frac{d{\dot{m}}_V}{dt}$$

(20)

where $u_V$ is the vapor velocity, ${\rho }_{VA}$ is the vapor density in the adiabatic section, and $A_V$ is the cross-sectional area of the vapor duct.

As the flow is assumed as laminar (first assumption), the second term of Equation (19) can be rewritten as:

$$-{\int }_{L_{eff}}\left(\mu ·\frac{{\partial }^{2}u}{\partial x^2}\right)dl={\rho }_{VA}·\frac{64}{Re}·\frac{L_{eff}}{d_h}·\frac{u_V^2}{2}=\frac{32{\mu }_{VA}{L}_{eff}}{{\rho }_{VA}{d}_h^2{A}_V}·{\dot{m}}_V$$

(21)

where $d_h$ is the hydraulic diameter of the vapor duct, ${\mu }_{VA}$ is the vapor dynamic viscosity in the adiabatic section, and $Re$ is the Reynolds number ($Re={\rho }_{VA}{d}_h{u}_V/{\mu }_{VA}$).

The third term of Equation (19) can be rewritten as:

$${\int }_{L_{eff}}\left(\frac{\partial p}{\partial x}\right)dl=p_{VC}-p_{VE}$$

(22)

Finally, the last term of Equation (19) can be rewritten as:

$${\int }_{L_{eff}}\left(\rho gsin\left(\alpha \right)\right)dl={\rho }_{VA}g{L}_{eff}sin\left(\alpha \right)$$

(23)

Therefore, the equation for the vapor line to be added in the complete set of equations for the fluid network is:

$$\frac{L_{eff}}{A_V}·\frac{d{\dot{m}}_V}{d\tau }+\frac{32{\mu }_{VA}{L}_{eff}}{{\rho }_{VA}{d}_h^2{A}_V}·{\dot{m}}_V=p_{VE}-p_{VC}-{\rho }_{VA}g{L}_{eff}sin\left(\alpha \right)$$

(24)

Equation (24) expresses the relation between the vapor mass flow rate, the pressures in the vapor tanks, and the gravitational contribution.

2.2.6. Liquid Line

Analyzing the liquid line, first of all, it is possible to write the liquid mass flow rate ${\dot{m}}_L$ as a function of the wick porosity $\epsilon$, that is, the liquid fraction capability of the wick, namely, the ratio between the volume of voids in the wick and its total volume:

$${\dot{m}}_L={\rho }_{LA}\epsilon A_W{u}_L$$

(25)

where $u_L$ is the liquid velocity; ${\rho }_{LA}$ is the liquid density in the adiabatic zone; and $A_W$ is the cross-sectional area of the wick.

Next, the pressure drop in a porous medium can be expressed through Darcy’s law [14]:

$$\frac{\mu }{K}{u}_x=-\frac{\partial p}{\partial x}+f_x$$

(26)

where $u_x$ is the velocity in the x direction, $\mu$ is the dynamic viscosity, $f_x$ is the body force in the x direction, and K is the wick permeability, which can be estimated by using the Blake–Kozeny equation [15,16]:

$$K=\frac{4r_p^2{\epsilon }^3}{150{\left(1-\epsilon \right)}^2}$$

(27)

where $r_p$ is the pore radius. Consequently, the equation describing the liquid line can be written as

$$\frac{{\mu }_{LA}{L}_{eff}}{K{\rho }_{LA}\epsilon A_W}·{\dot{m}}_L=p_{LC}-p_{LE}+{\rho }_{LA}g{L}_{eff}sin\left(\alpha \right)$$

(28)

Therefore, the complete system of equations which describe the behavior of the fluid temperatures, pressures, and mass flow rates in the fluid network, is composed of Equations (10), (24), (28) and Systems (14), (16), (17).

2.3. Solid–Fluid Coupling

The solid–fluid coupling is described by four equations. The first two equations describe the coupling between the wick temperature and the liquid temperature, which, given the forth assumption in Section 2.2 are always equal, so that the equations read:

$$\left\{\begin{array}{c}{T}_{LE}=T_{WE}\hfill \\ T_{LC}=T_{WC}\hfill \end{array}\right.$$

(29)

where $T_{LE}$ and $T_{LC}$ are the liquid temperatures in the evaporator and condenser zones, respectively.

The other two equations for the solid–fluid coupling are those that define the heat transfer rates ${\dot{Q}}_E$ and ${\dot{Q}}_C$, as functions of some temperatures and resistances, which are taken from both the solid and fluid networks. In particular, by looking at Figure 1, it is possible to state that:

$$\left\{\begin{array}{c}{\dot{Q}}_E=\frac{T_{WE}-T_{VE}}{R_{2WE}+R_{EF}}\hfill \\ {\dot{Q}}_C=\frac{T_{VC}-T_{WC}}{R_{2WC}+R_{CF}}\hfill \end{array}\right.$$

(30)

The two equations of System (30) describe the mathematical expression of the heat transfer rates that is necessary to give (${\dot{Q}}_E$) or extract (${\dot{Q}}_C$) from the fluid, in order to make it evaporate or condense.

Therefore, the four equations of Systems (29) and (30) constitute the complete set of equations which describe the coupling between the fluid and the solid part of the heat pipe model.

2.4. Liquid–Vapor Coupling

The liquid–vapor coupling is described by only one equation, namely, the Laplace–Young equation [17]. It describes the capillary pressure jump in the evaporator as a function of the liquid surface tension at the evaporator temperature, $\sigma$, the contact angle, $\theta$, and the capillary radius, $r_c$, as:

$$p_{VE}-p_{LE}=\frac{2\sigma cos\left(\theta \right)}{r_c}$$

(31)

In Equation (31), the contact angle $\theta$ is assumed to be constant during the operation (neglecting contact angle hysteresis), but dependent on the fluid–solid couple, while the effective capillary radius, $r_c$, is computed as a function of the pore radius, $r_p$, as [8]:

$$r_c=0.41r_p$$

(32)

Equation (31) is the last one to be taken into account for the description of the complete model.

2.5. Implementation of the Complete Set

From the previous sections, it is possible to conclude that the heat pipe behavior is described by twenty equations, namely, eleven ordinary differential equations and nine algebraic equations, which contain twenty unknowns. Therefore, the system is fully determined.

The equations are:

$$\left\{\begin{array}{c}{C}_{PE}·\frac{d{T}_{PE}}{d\tau }={\dot{Q}}_{IN}+\frac{T_{PA}-T_{PE}}{R_{1PA}}+\frac{T_{WE}-T_{PE}}{R_{2PE}+R_{1WE}}\hfill \\ C_{PC}·\frac{d{T}_{PC}}{d\tau }=-{\dot{Q}}_{OUT}+\frac{T_{PA}-T_{PC}}{R_{2PA}}+\frac{T_{WC}-T_{PC}}{R_{2PC}+R_{1WC}}\hfill \\ C_{WE}·\frac{d{T}_{WE}}{d\tau }=\frac{T_{PE}-T_{WE}}{R_{1WE}+R_{2PE}}+\frac{T_{WA}-T_{WE}}{R_{1WA}}+\frac{T_{VE}-T_{WE}}{R_{2WE}+R_{EF}}\hfill \\ C_{WC}·\frac{d{T}_{WC}}{d\tau }=\frac{T_{PC}-T_{WC}}{R_{1WC}+R_{2PC}}+\frac{T_{WA}-T_{WC}}{R_{2WA}}+\frac{T_{VC}-T_{WC}}{R_{2WC}+R_{CF}}\hfill \\ C_{PA}·\frac{d{T}_{PA}}{d\tau }=\frac{T_{PE}-T_{PA}}{R_{1PA}}+\frac{T_{PC}-T_{PA}}{R_{2PA}}+\frac{T_{WA}-T_{PA}}{R_{PA}+R_{WA}}\hfill \\ C_{WA}·\frac{d{T}_{WA}}{d\tau }=\frac{T_{PA}-T_{WA}}{R_{WA}+R_{PA}}+\frac{T_{WE}-T_{WA}}{R_{1WA}}+\frac{T_{WC}-T_{WA}}{R_{2WA}}+\frac{T_{VA}-T_{WA}}{R_{3WA}+R_{AF}}\hfill \\ {\dot{m}}_L={\dot{m}}_C\hfill \\ \frac{V_{VE}}{a_{VE}^2}·\frac{d{p}_{VE}}{d\tau }={\dot{m}}_E-{\dot{m}}_V\hfill \\ \frac{V_{VC}}{a_{VC}^2}·\frac{d{p}_{VC}}{d\tau }={\dot{m}}_V-{\dot{m}}_C\hfill \\ \frac{d{T}_{VE}}{d\tau }=T_{VE}\left(\frac{k_{pVE}}{{\rho }_{VE}·c_{pVE}}\right)\frac{d{p}_{VE}}{d\tau }\hfill \\ \frac{d{T}_{VC}}{d\tau }=T_{VC}\left(\frac{k_{pVC}}{{\rho }_{VC}·c_{pVC}}\right)\frac{d{p}_{VC}}{d\tau }\hfill \\ {\dot{m}}_E=\frac{{\dot{Q}}_E}{h_{LVE}}\hfill \\ {\dot{m}}_C=\frac{{\dot{Q}}_C}{h_{LVC}}\hfill \\ \frac{L_{eff}}{A_V}·\frac{d{\dot{m}}_V}{d\tau }+\frac{32{\mu }_{VA}{L}_{eff}}{{\rho }_{VA}{d}_h^2{A}_V}·{\dot{m}}_V=p_{VE}-p_{VC}-{\rho }_{VA}g{L}_{eff}sin\left(\alpha \right)\hfill \\ \frac{{\mu }_{LA}{L}_{eff}}{K{\rho }_{LA}\epsilon A_W}·{\dot{m}}_L=p_{LC}-p_{LE}+{\rho }_{LA}g{L}_{eff}sin\left(\alpha \right)\hfill \\ T_{LE}=T_{WE}\hfill \\ T_{LC}=T_{WC}\hfill \\ {\dot{Q}}_E=\frac{T_{WE}-T_{VE}}{R_{2WE}+R_{EF}}\hfill \\ {\dot{Q}}_C=\frac{T_{VC}-T_{WC}}{R_{2WC}+R_{CF}}\hfill \\ p_{VE}-p_{LE}=\frac{2\sigma cos\left(\theta \right)}{r_c}\hfill \end{array}\right.$$

(33)

The unknowns can be classified as:

• Ten temperatures: $T_{PE},T_{PC},T_{WE},T_{WC},T_{PA},T_{WA},T_{VE},T_{VC},T_{LE},T_{LC}$;
• Four pressures: $p_{VE},p_{VC},p_{LE},p_{LC}$;
• Four mass flow rates: ${\dot{m}}_E,{\dot{m}}_C,{\dot{m}}_V,{\dot{m}}_L$;
• Two heat transfer rates: ${\dot{Q}}_E,{\dot{Q}}_C$.

In order to solve System (33), an explicit method is used. In particular, the eleven differential equations are discretized by the forward finite differences approach and solved for each time step. Then, the remaining nine algebraic equations are sequentially solved.

Since the forward finite differences approach is an explicit method, and since the number of equations to be solved is quite large, the stability of the solution requires a very low time step. However, this does not greatly affect the computational time, which remains absolutely manageable.

In the solution process, the thermophysical properties of the liquid (saturation pressure, density, viscosity, thermal conductivity), the surface tension, and the latent heat of vaporization are evaluated by using a polynomial interpolation of the saturation table of the working fluid (fifth assumption in Section 2.2). Conversely, the thermophysical properties of the vapor (density, viscosity, specific heat at constant pressure, isobaric expansion coefficient, sound speed) are evaluated at the current temperature and pressure, with the aid of some look-up tables. Such tables were created using the REFPROP software (version 9.1), from the National Institute of Standards and Technology [18]. The data contained in these tables are extracted with a nearest neighbor approach, in order to obtain meaningful results also near saturation conditions.

Further details about the numerical implementation method can be found in Appendix A.

3. Case Study

As already said, the model was first of all validated against experimental results acquired by Leonardo S.p.A. The heat pipe used for the experimental test is shown in Figure 2a. It is a copper heat pipe using water as the working fluid. In order to be easily brazed to the aluminum chassis, the HP was covered by a very thin electrolytic nickel plate (thickness $5÷8\mu \mathrm{m}$).

The external dimensions of the HP are shown in Figure 2b. It is worth mentioning that since the two ends of the heat pipe are not straight, but sharpened, the total length that will be given as an input to the model is not $120\mathrm{m}\mathrm{m}$, but $115\mathrm{m}\mathrm{m}$, in order to compensate for the loss of thermal exchange area due to the sharpening.

From Figure 2b, it is also possible to notice that the section of the HP is not cylindrical. In fact, the thickness is much lower than the width. Concerning the internal dimensions:

• The thickness of the wall is equal to $0.296\mathrm{m}\mathrm{m}$;
• The thickness of the wick is equal to $0.500\mathrm{m}\mathrm{m}$;
• The wick porosity, $\epsilon$, is approximately equal to $0.5$;
• The average grain radius is approximately equal to 100 μm.

3.1. Description of the Experiment

In the experimental test, three heat pipes like the one shown in Figure 2 were embedded in a chassis cover made of the aluminum alloy 6082, as shown in Figure 3. In particular, the empty chassis with the three cavities can be seen in Figure 3a, the heat pipe ready to be embedded in a cavity can be seen in Figure 3b.

After the heat pipe embedding, a $TGHLV1R00JE$ (Ohmite Manufacturing Company) resistor was placed on the evaporator zones of the three heat pipes, to be used as thermal source. The resistor was fixed with a bracket to ensure mechanical and thermal contact of the power component.

Then, the front face of the chassis was entirely blackened, in order to perform a proper measurement by using an infrared thermocamera, and two thermocouples were placed on the two opposite sides of the central heat pipe.

Given the placing of the setup, it is possible to state that the heat pipes have a vertical bottom-heating orientation, with the condenser over the evaporator. Therefore, the orientation angle $\alpha$ to be used in the model is equal to 90°.

Furthermore, the heat pipe condenser zone is able to exchange thermal power by conduction with the chassis cover, by natural convection with the air, and by radiation with the room walls. Thus, the resistance connecting the condenser external wall and the external environment, $R_M$, is a mixed resistance.

The experiment was performed by turning on the resistor for $60\mathrm{s}$, then turning it off for $5\min$, and repeating the same cycle of thermal charge and discharge six more times, for a total of seven cycles, and a total time of $2520\mathrm{s}$.

The total thermal power given by the resistor in the thermal charge phases is $115\mathrm{W}$, but only a part of it is taken by the heat pipes. In fact, a great part is directly dissipated by the chassis cover to the external environment.

The results of this experiment, namely, the temperatures acquired by the thermocouples, are the output to be compared with that from the model object of this study. However, before this step, it is necessary to adapt the latter to the geometry of the real heat pipe, and to extract the data of ${\dot{Q}}_{IN}$ and $R_M$ to be given to the model as input parameters.

3.2. Model Adaptation

For the aim of this work, the shape of the heat pipe cross-section can be modeled as a rectangle with two semicircles at the two ends, as can be seen in Figure 4. Therefore, the formulas used for computing the resistances and capacitances of the solid network, and the formulas for the quantities multiplying $\frac{d{\dot{m}}_V}{d\tau }$, ${\dot{m}}_V$, and ${\dot{m}}_L$ in the fluid network equations must be adapted to this geometry.

First, it is useful to define the cross-sectional areas of the wall, $A_P$, wick, $A_W$, and vapor zone, $A_V$, in the new domain:

$$\left\{\begin{array}{c}{A}_P=\pi \left(r_{ep}^2-r_{ip}^2\right)+2w_{eff}{t}_p\hfill \\ A_W=\pi \left(r_{ip}^2-r_{iw}^2\right)+2w_{eff}{t}_w\hfill \\ A_V=\pi r_{iw}^2+2w_{eff}{r}_{iw}\hfill \end{array}\right.$$

(34)

where:

• $r_{ep}$, $r_{ip}$, and $r_{iw}$ are the external wall, internal wall, and internal wick radii of the circular regions at the ends of the section, respectively;
• $w_{eff}$ is the effective width of the section, which corresponds to the total width minus twice the external wall radius: $w_{eff}=w_{tot}-2r_{ep}$;
• $t_p$ is the wall thickness: $t_p=r_{ep}-r_{ip}$;
• $t_w$ is the wick thickness: $t_w=r_{ip}-r_{iw}$.

At this point, it is possible to define all the aforementioned quantities in this specific domain.

3.2.1. Solid Network

For the capacitances of the solid network, the adaptation is straightforward.

The thermal capacitance, $C_n$, of a node n is defined as the product of the mass of the node, $m_n$, and its specific heat, $c_n$. In its turn, the mass of the node can be expressed as the product of its density, ${\rho }_n$, and its volume $V_n$. Thus, the thermal capacitance $C_n$ can be generically expressed as:

$$C_n={\rho }_n{c}_n{V}_n$$

(35)

Starting from Equation (35), two cases have to be considered: the node n is in the solid wall ($C_{PE}$, $C_{PA}$, $C_{PC}$), and the node n is in the wick ($C_{WE}$, $C_{WA}$, $C_{WC}$).

In case n is in the solid wall, the evaluation of the thermal capacitance is trivial. In fact, ${\rho }_n$ is the density of the material that constitutes the solid wall, ${\rho }_W$. The specific heat $c_n$ is, in the same way, that of the material that constitutes the solid wall, $c_W$. The volume $V_n$ can be computed as the product of the longitudinal dimension of the node, $L_n$, and the area of the node section, $A_n$, which is $A_P$ for the considered case. Therefore, the thermal capacitance of a node in the solid wall can be computed as:

$$C_P={\rho }_W{c}_W{L}_n{A}_P={\rho }_W{c}_W{L}_n\left(\pi \left(r_{ep}^2-r_{ip}^2\right)+2w_{eff}{t}_p\right)$$

(36)

In case n is in the wick, since the wick is composed of a solid and a liquid part, it is necessary to compute the effective value of the product ${\rho }_n{c}_n$, which is an average of the one of the solid, ${\rho }_W{c}_W$, and that of the liquid, ${\rho }_L{c}_L$, whose properties vary depending on the temperature. This average is weighted by the wick porosity, $\epsilon$. Thus, in this case:

$${\rho }_n{c}_n={\left(\rho c\right)}_{eff}=\epsilon ·{\rho }_L{c}_L+\left(1-\epsilon \right)·{\rho }_W{c}_W$$

(37)

The volume $V_n$ is computed exactly as in the previous case, but the area of the node, $A_n$, is equal to $A_W$. Thus, the thermal capacitance of a node in the wick can be written as:

$$C_W={\left(\rho c\right)}_{eff}{L}_n{A}_W=\left(\epsilon ·{\rho }_L{c}_L+\left(1-\epsilon \right)·{\rho }_W{c}_W\right)L_n\left(\pi \left(r_{ip}^2-r_{iw}^2\right)+2w_{eff}{t}_w\right)$$

(38)

At this point, knowing the heat pipe dimensions and the thermophysical properties of its materials, all the six thermal capacitances shown in Figure 1 can be easily computed.

As previously mentioned, the resistances can be axial or radial; conductive or convective (radiation inside the heat pipe has been neglected); and, as with the capacitances, in the solid wall or in the wick.

The convective resistances are all radial, and they are evaluated starting from the convective coefficient, $h_F$, and the area of the thermal exchange section, $A_F$, namely:

$$R_F=\frac{1}{h_F{A}_F}=\frac{1}{h_F\left(2\pi r_{iw}+2w_{eff}\right)L_n}$$

(39)

These convective resistances are assumed to be constant and known. In fact, it is assumed that the variations of $h_F$ during the heat pipe operation are negligible. Convective coefficients during phase change and in the adiabatic section were estimated by literature correlations [19].

The mixed resistance from the external wall to the external environment, $R_M$, can be either constant or variable during the heat pipe operation. This resistance can be evaluated using Equation (39), by estimating a possibly time-varying equivalent convective coefficient, or through empirical correlations, or using Equation (40), if the thermal power of output, ${\dot{Q}}_{OUT}$, is known for each time step:

$$R_M\left(\tau \right)=\frac{T_C\left(\tau \right)-T_F}{{\dot{Q}}_{OUT}\left(\tau \right)}$$

(40)

The conductive resistances must be studied separately in case they are axial or radial.

Starting from the axial case, for each type of geometry of the section, the axial resistance, $R_{A12}$, connecting two nodes, 1 and 2, can be evaluated as:

$$R_{A12}=\frac{L_{12}}{\lambda A}$$

(41)

where $L_{12}$ is the distance between node 1 and node 2, $\lambda$ is the thermal conductivity of the material crossed by the resistance, and A is the cross-sectional area of the considered section, which can be either the wall area, $A_P$, if the considered resistance is in the wall region, or the wick area, $A_W$, if the considered resistance is in the wick region.

Considering the radial case and starting from the effective thermal conductivity, there is a difference between those for the wall and for the wick. In case of a wall resistance, $\lambda$ is the thermal conductivity of the material that constitutes the solid wall, namely, $\lambda ={\lambda }_W$. Conversely, in case of a wick resistance, $\lambda$ can be evaluated as a weighted average of the wall thermal conductivity, ${\lambda }_W$, and the liquid thermal conductivity, ${\lambda }_L\left(\tau \right)$, which varies at each time step depending on the liquid temperature. For a sintered wick, the effective thermal conductivity of the wick ${\lambda }_{eff}$ can be computed as [8]:

$${\lambda }_{eff}\left(\tau \right)=\frac{{\lambda }_W\left(2+\frac{{\lambda }_L\left(\tau \right)}{{\lambda }_W}-2\epsilon \left(1-\frac{{\lambda }_L\left(\tau \right)}{{\lambda }_W}\right)\right)}{2+\frac{{\lambda }_L\left(\tau \right)}{{\lambda }_W}+\epsilon \left(1-\frac{{\lambda }_L\left(\tau \right)}{{\lambda }_W}\right)}$$

(42)

Then, the expression for the radial resistances of the wall and of the wick are computed by dividing the heat pipe section into four parts: two of them are concentric semicircles, the other two are flat plates. Each of these parts has a radial resistance, and all these resistances are in parallel configuration. However, the resistances of the two concentric semicircles are equal, therefore, their equivalent resistance, $R_{eqC}$, is equal to the value of each single resistance divided by two, which corresponds to the resistance of two concentric circles. Moreover, the resistances of the two flat plates are equal too, therefore, their equivalent resistance, $R_{eqL}$, is equal to the value of each single resistance divided by two, which corresponds to the resistance of a flat plate with double section area. At this point, only two resistances, $R_{eqC}$ and $R_{eqL}$, in parallel configuration, are left, so it is straightforward to compute the final equivalent resistance, namely, the radial resistance of the entire section, $R_R$.

In particular, the resistance of each of the concentric semi-circles, $R_C$, is:

$$R_C=\frac{log\left(\frac{r_e}{r_i}\right)}{\pi \lambda L}$$

(43)

where $r_e$ and $r_i$ are the external and internal radii delimiting the region under consideration, respectively, L is the length of such region, $\lambda$ is equal to ${\lambda }_W$ for the wall, and ${\lambda }_{eff}$ for the wick.

Therefore, the equivalent radial resistance of the two concentric semi-circles is:

$$R_{eqC}=\frac{R_C}{2}=\frac{log\left(\frac{r_e}{r_i}\right)}{2\pi \lambda L}$$

(44)

The resistance across a flat plate of width $w_{eff}$ and thickness $t=r_e-r_i$, is:

$$R_L=\frac{t}{\lambda w_{eff}L}=\frac{r_e-r_i}{\lambda w_{eff}L}$$

(45)

Therefore, the equivalent radial resistance of the two flat plates is:

$$R_{eqL}=\frac{R_L}{2}=\frac{r_e-r_i}{2\lambda w_{eff}L}$$

(46)

At this point, it is possible to compute the equivalent radial resistance of the complete section as:

$$R_R=R_{eqC}//R_{eqL}=\frac{R_{eqC}·R_{eqL}}{R_{e{q}_C}+R_{eqL}}=\frac{log\left(\frac{r_e}{r_i}\right)·\left(r_e-r_i\right)}{2\lambda L\left(\pi \left(r_e-r_i\right)+w_{eff}log\left(\frac{r_e}{r_i}\right)\right)}$$

(47)

3.2.2. Fluid Network

The expression of the coefficient that multiplies $\frac{d{\dot{m}}_V}{dt}$ in Equation (24) is quite easy to adapt. In fact, looking at Equation (20), it is possible to notice that it is the ratio between the effective duct length and the cross-sectional area of the vapor duct. Therefore, in the new domain:

$$\frac{L_{eff}}{A_V}=\frac{L_{eff}}{\pi r_{iw}^2+2w_{eff}{r}_{iw}}$$

(48)

The expression of the coefficient multiplying ${\dot{m}}_V$ in the equation for the vapor line, namely, again, Equation (24), involves the computation of the hydraulic diameter of the vapor section, $d_h$. In general, the hydraulic diameter is defined as the ratio between four times the cross-sectional area of the vapor duct, $A_V$, and the duct perimeter, $P_V$. Therefore, for the considered geometry, it can be computed as:

$$d_h=\frac{4A_V}{P_V}=\frac{4\left(\pi r_{iw}^2+2r_{iw}{w}_{eff}\right)}{2\left(\pi r_{iw}+w_{eff}\right)}=\frac{2\left(\pi r_{iw}^2+2r_{iw}{w}_{eff}\right)}{\pi r_{iw}+w_{eff}}$$

(49)

In this way, starting from Equation (21), the coefficient becomes:

$$\frac{32{\mu }_{VA}{L}_{eff}}{{\rho }_{VA}{d}_h^2{A}_V}=\frac{8{\mu }_{VA}{L}_{eff}{\left(\pi r_{iw}+w_{eff}\right)}^2}{{\rho }_{VA}{\left(\pi r_{iw}^2+2r_{iw}{w}_{eff}\right)}^3}$$

(50)

Finally, from Equation (28), it is possible to notice that the coefficient multiplying ${\dot{m}}_L$ in the equation for the liquid line can be expressed as:

$$\frac{{\mu }_{LA}}{K{\rho }_{LA}\epsilon A_W}=\frac{{\mu }_{LA}}{K{\rho }_{LA}\epsilon \left(\pi \left(r_{ip}^2-r_{iw}^2\right)+2w_{eff}{t}_w\right)}$$

(51)

3.3. Extraction of Input Data

After the adaptation of the model to the real geometry, it is necessary to find the input data for the model, namely, the initial conditions and the boundary conditions.

The only initial condition to define is the initial temperature, $T_0$. In fact, the initial mass flow rate is zero by definition and the initial pressure, $p_0$, is computed by inverting the Kelvin equation [20], as:

$$p_0=\frac{p_{sat0}}{e^{\left(\frac{2{\sigma }_0{V}_m}{R{T}_0{r}_m}\right)}}$$

(52)

where:

• $p_{sat0}$ is the saturation pressure evaluated at temperature $T_0$;
• ${\sigma }_0$ is the liquid surface tension evaluated at temperature $T_0$;
• $V_m$ is the liquid molar volume, which is computed dividing the molar mass of the working fluid by the liquid density evaluated at temperature $T_0$;
• R is the universal gas constant;
• $r_m$ is the mean radius of curvature of the liquid–gas interface, and for this application, it is equal to the capillary radius, $r_c$.

In this case, the initial temperature is taken from the experimental data, and it is equal to the ambient temperature, 19.5 °C.

Concerning the boundary conditions, also the temperature of the external environment, $T_F$, is provided by the experimental data for all the time steps of the thermal test. Even if this temperature varies during the experiment, its variations are so small that they can be neglected. Therefore, $T_F$ is assumed to be constant and equal to its measured value at the first time step.

At this point, the only missing data are the input thermal power, ${\dot{Q}}_{IN}$, and the resistance between the condenser wall and the external environment, $R_M$. In order to obtain their values, a preliminary computational fluid dynamics simulation was performed using the finite element approach by means of Comsol Multiphysics [21]. The idea is to set a very simple model of the chassis cover, heat pipes, and thermal source, to run a basic simulation in order to find the temperature evolution of the heat pipe in the middle, and to process the temperature data obtained in order to find the values of ${\dot{Q}}_{IN}$ and $R_M$.

In particular, at first, the geometry is set. The heat pipes are modeled as conductive bars with the following features:

• An equivalent mass, computed as the sum of the masses of the solid, liquid, and vapor parts;
• An equivalent thermal conductivity, estimated by a preliminary run of the heat pipe model in a basic case;
• An equivalent thermal capacitance, computed as weighted average of the capacitances of the solid, liquid, and vapor parts.

The heat source is modeled as a rectangular box, in contact the lower part of the three heat pipes, and whose external dimensions are the ones of the $TGHLV1R00JE$ resistor. The chassis is modeled as a nearly rectangular box with raised edges, made of aluminum.

After defining the geometry, the domain is meshed with a standard tetrahedral mesh. Both the geometry and the mesh are shown in Figure 5.

The heat pipes can exchange heat in different ways:

• Their sides and rear face can exchange by conduction with the chassis cover;
• Their front face can exchange by natural convection with the air in the room;
• Their front face can exchange also by radiation with the walls of the room.

Inside the chassis cover, heat is obviously transported by conduction. Then, all the external surfaces of the chassis cover can exchange heat by convection and radiation with the surroundings. Radiation is modeled with the simplifying assumption of “small body exchanging with a large cavity”. Emissivity of the surfaces is assumed as $0.1$ for the nontreated aluminum parts, while it is assumed as $0.93$ for the blackened faces. The convective coefficient between the chassis cover and the surrounding air is roughly assumed as uniform on all faces, and estimated by correlations for natural convection over a vertical flat plate (resulting in a value of $4\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$).

The input thermal power given by the resistor is modeled as a periodic rectangular wave with $1/6$ duty cycle ($60\mathrm{s}$ pulse width, $360\mathrm{s}$ total cycle duration) and maximum power $115\mathrm{W}$. The results of the simulation provide the temperature evolution of the whole chassis with the embedded heat pipes.

The temperature distribution on the chassis at four significant instants of time is shown in Figure 6 (at the end of the two phases of the first cycle) and Figure 7 (at the end of the two phases in a stabilized cycle). In particular, Figure 6a shows the temperature at the end of the first thermal charge phase, Figure 6b shows the temperature at the beginning of the second thermal charge phase, Figure 7a shows the temperature at the end of the last thermal charge phase, and Figure 7b shows the temperature at the end of the last thermal discharge phase.

The temperature evolution of the central heat pipe is extracted and given as an input to a MATLAB code. In this code, the heat pipe is modeled by using a single-node lumped parameter approach, with the following assumptions:

• The input thermal power, ${\dot{Q}}_{IN}$, is constant and equal for all the thermal charge phases, and it is null for all the thermal discharge phases;
• The value of the resistance $R_M$ is constant during the phases of charge and discharge of the same cycle, while it changes between the cycles to take into account the equivalent convective coefficient variation with temperature.

Therefore, the two equations describing the temperature evolution of the heat pipe as a function of ${\dot{Q}}_{IN}$ and $R_M$ are:

• For the thermal charge phases: $C·\frac{dT}{d\tau }={\dot{Q}}_{IN}-\frac{T-T_F}{R_M}$;
• For the thermal discharge phases: $C·\frac{dT}{d\tau }=-\frac{T-T_F}{R_M}$.

In the above, C is the equivalent heat capacity of the heat pipe.

In the previous equations, the equivalent heat capacity, the temperature, and the temperature derivative are assumed as known for each time step—they are obtained by sampling the temperature profile in time calculated by the transient CFD simulation in two points, corresponding to the positions of the thermocouples in the real experiment, and averaging them.

Consequently, for each thermal cycle, the MATLAB code computes $R_M$ from the second equation and then ${\dot{Q}}_{IN}$ from the first equation, thus returning the two needed pieces of input data.

4. Results

In this section, the results obtained using the described model will be presented. The section is divided in two parts: the first one describes the comparison between numerical and experimental data for the case study, while the second one focuses on the analysis of the model accuracy when using the ideal and real gas assumptions, for a representative simulation with constant input power, using different fluids.

4.1. Model Validation

By using the input parameters described in Section 3.3, the simulation is run. As previously mentioned, the total time for the seven cycles is $2520\mathrm{s}$, which corresponds to $42\min$. Since the ODEs of the model are approximated through the forward finite differences method, which is an explicit method, a small time step is needed to keep the simulations stable. After some preliminary tests, a fixed time step equal to 5 μs was selected.

In the experiment, two thermocouples were placed on the two opposite extremities of the central heat pipe. In particular, the one closer to the resistor acquires a temperature that will be indicated with the symbol $T_a$ in the following, while the other one acquires a temperature that will be indicated with the symbol $T_b$ in the following.

Since the thermocouples are placed on the external wall of the heat pipe, at the two opposite extremities, it is reasonable to compare $T_a$ with the evaporator wall temperature, $T_E$, and $T_b$ with the condenser wall temperature, $T_C$. It is worth noting that the results are partially biased by the fact that the real evaporator wall temperature is the one of the heat pipe face below the resistor, so it is reasonably higher than $T_a$. Moreover, the condenser wall temperature computed by the model is the temperature at the center of the condenser wall. Thus, since more than the $70\%$ of the heat pipe length is taken up by the condenser, its center is not close to its end. Consequently, in practice, it must be taken into account that $T_C$ cannot be exactly equal to $T_b$.

Figure 8 shows the comparison between numerical data computed through the model and experimental data acquired by the thermocouples for the seven thermal cycles. In particular, Figure 8a shows the comparison between $T_a$ and $T_E$, while Figure 8b shows the comparison between $T_b$ and $T_C$.

From Figure 8, it is possible to see that the model provides results very close to the experimental ones in the last three cycles, namely, the stabilized ones. In fact, in these cycles, the error on the peaks of the evaporator temperature goes from $2.26$ to $3.10\%$, with an average value of $2.68\%$, while the error on the peaks of the condenser temperature covers the range between $4.66$ and $6.27\%$, with an average value of $5.36\%$. Conversely, in the first four cycles, the maximum temperature of each cycle estimated by the model is much higher than the measured one, even if the minimum temperature for each cycle is very close to the experimental one. In these cycles, the error on the peaks of the evaporator temperature goes from $21.9$ to $25.0\%$, with an average value of $23.4\%$, while the error on the peaks of the condenser temperature covers the range between $13.7$ and $24.0\%$, with an average value of $17.0\%$.

The temperature evolution computed by the model is more regular than the measured one. This can be explained by the fact that when the heat pipe operation starts, the activation process takes a certain time, which varies depending on the considered heat pipe. In its present version, the activation process in the model is idealized, no activation resistances are considered.

In the experimental test analyzed, the heat pipe activation process reasonably takes a significant fraction of the total time for the first thermal charge phases, so, at the end, the heat pipe has not reached a temperature as high as it would have with a lower or null activation time. However, after four cycles, the temperature reached is high enough to stabilize the following cycles, so the model becomes much more reliable.

From Figure 8b, it is possible to notice that in the last three cycles, the measured temperature is higher than the one predicted by the model. This discrepancy is probably due to how the input data are extracted. In fact, after a preliminary CFD simulation, the boundary conditions are obtained by approximating the heat pipe through a single-node lumped parameter approach. When using this method, it is assumed that each region of the heat pipe exchanges the same amount of heat. However, in the real case, it is reasonable to assume that the heat exchanged by the front face of the condenser, by natural convection and radiation, is lower than the heat exchanged by the other faces of the condenser, by conduction with the chassis cover. Therefore, the temperature of the condenser external wall, which can be computed as an average of the temperatures of the whole external wall of the condenser, may be lower than $T_b$.

Therefore, it is possible to state that, in view of the results of the comparison between the numerical data computed through the model and experimental data acquired by the thermocouples for the seven thermal cycles, the validation of the model can be considered satisfactory, despite some predictable limitations.

4.2. Comparison between Ideal and Real Gas

In this subsection, four fictitious cases with constant evaporator heat power are considered. The heat pipe considered is the one used for the case study in terms of shape and dimensions, but the analysis is performed by considering four different solid–fluid couples, namely, copper–water, aluminum–ammonia, aluminum–acetone, and copper–HFC134a refrigerant. For each case, the vapor temperatures and pressures obtained by using the model proposed in this work are compared with the ones obtained by modeling the vapor as an ideal gas, as in the work of Ferrandi et al. [8].

In order to perform such analysis, it is assumed that the considered heat pipe is subject to a constant heat power input equal to $20\mathrm{W}$, and the resistance $R_M$ is modeled as an equivalent convective resistance. The initial temperature and the external temperature are both set to be equal to 25 °C. Furthermore, it is assumed, without loss of generality, that the heat pipe is oriented horizontally.

4.2.1. Copper-Water

Figure 9 shows the comparison between the vapor temperatures in case of using ideal or real gas equations for the vapor modeling, in a heat pipe made of copper and filled with water as working fluid. In particular, Figure 9a shows the comparison between the vapor temperatures in the evaporator, $T_{VE}$, while Figure 9b shows the comparison between the vapor temperatures in the condenser, $T_{VC}$.

It is evident how the vapor temperatures in the two cases are practically equal. Looking at the vapor pressures in the evaporator, $p_{VE}$ in Figure 10a, and in the condenser, $p_{VC}$ in Figure 10b, it is possible to notice a small difference in the results when the steady state is reached.

Nevertheless, this difference leads to an error lower than $1\%$ in the pressure estimation. Therefore, for the given conditions, modeling the vapor as a real gas is not necessary when the considered fluid is water. This is reasonable, as the compressibility factor of water in the considered case is very close to the one in [22], confirming the fact that water vapor behaves almost as an ideal gas.

4.2.2. Aluminum–Ammonia

As a second case, the same simulation was repeated using a heat pipe made of aluminum and filled with ammonia as working fluid.

The difference between the vapor temperatures computed by approximating the vapor as ideal or real gas are shown in Figure 11.

Again, it is possible to notice that the temperature profiles are similar to the point that it is hard to distinguish the respective curves.

The comparison between the vapor pressure curves in case of ideal and real gas shows a difference of about $1.7\%$, as can be seen in Figure 12. This difference can be relevant in some applications, but, in general, such an error will be much lower than those caused by the use of the lumped parameter approach, so for this case also it can be concluded that there is no relevant difference in modeling the vapor as ideal or real gas.

4.2.3. Aluminum–Acetone

At this stage, the same case as before is analyzed, but using a heat pipe made of aluminum and filled with acetone as working fluid.

The comparison between the vapor temperatures obtained in case of modeling the vapor as ideal or real gas is shown in Figure 13 and, as in the previous cases, the difference between the curves is negligible.

However, comparing the vapor pressures in the two cases (Figure 14), it is possible to notice that the error in modeling the vapor as an ideal gas rather than a real gas is close to $18\%$, and this is certainly not negligible.

Therefore, when acetone is used as working fluid, it is necessary to model the vapor by using the real gas equations.

4.2.4. Copper–HFC134a Refrigerant

Finally, again the same case is analyzed, but using a heat pipe made of copper and filled with HFC134a as working fluid.

The results are the same of the acetone–aluminum case: the temperature profiles (Figure 15) show a practically null difference, while the pressure profiles (Figure 16) evidence a deviation of more than $13\%$.

Therefore, also when using HFC134a as working fluid of the heat pipe, modeling the vapor as a real gas offers a significant increase in the accuracy.

In summary, the results confirmed that when the working fluid and the operating conditions are such that the intermolecular forces play a significant role, real gas equations and properties should be used. The amount of deviation from the ideal gas behavior at a certain temperature and pressure depends on the specific substance. To better evidence this aspect, Figure 17 shows the paths followed by the four fluids (water, ammonia, acetone, HFC134a) in the evaporator, during the transition from start-up to steady operation, in a $p-T$ chart (using reduced thermodynamic properties, $p_r$ and $T_r$).

As it can be seen, the fluids work in different regions. The average ratios between real gas and ideal gas values for compressibility, speed of sound, isobaric specific heat, and isobaric expansion coefficient are also reported. The latter evidence how for some fluids the values are all close to unity, and this is reflected in a negligible difference between real and ideal gas predictions, while for others the difference is large, resulting in the need for real gas modeling.

5. Conclusions

In this study, a multinode lumped parameter model for steady and transient analysis of heat pipes was developed, starting from other models already existing in literature, but including the effects of gravity (orientation angle) and modeling the vapor as a real gas.

The model was validated by using an experimental case study, in which three parallel heat pipes were embedded in a chassis cover and heated by a resistor, during seven cycles of $1\min$ of thermal charge and $5\min$ of thermal discharge.

The evaporator and condenser wall temperatures of the heat pipe in the middle, measured by two thermocouples, were compared to the ones obtained through the model. The results showed good agreement in predicting data in the stabilized thermal cycles, while the model provides temperatures much higher than the measured ones in the first cycles. This is probably due to the heat pipe activation process, the method through which the boundary conditions were calculated, and also to the position of the thermocouples which acquired the experimental results.

Additionally, it is analyzed whether, for some working fluids, it is important to model the vapor as a real gas, rather than as an ideal gas. This analysis showed that the vapor temperatures obtained by modeling the vapor as an ideal or real gas are almost the same, for all the considered cases. However, the vapor pressures can be significantly different.

In particular, when using water or ammonia as working fluids, the difference between the vapor pressures is negligible, while when using acetone or HFC134a as working fluids, modeling the vapor as an ideal gas rather than a real gas leads to an error in the vapor pressure prediction larger than $10\%$.

Therefore, it can be concluded that when the operating conditions are such that the intermolecular forces play a significant role, it becomes important to model the working fluid as a real gas, both in terms of relations and of values of the thermophysical properties included in the model. Otherwise, the latter may be significantly under- or overestimated, resulting in a bad prediction, particularly of the vapor pressures.