Numerical Investigation with Experimental Validation of Heat and Mass Transfer during Evaporation in the Porous Wick within a Loop Heat Pipe

: The heat transfer performance of the evaporator signiﬁcantly affects the heat transfer capacity of the loop heat pipe (LHP). The vapor blanket can be formed once the vapor penetrates the wick especially at high heat ﬂux, resulting in an unsaturated state of the wick and deteriorating the evaporator performance. It is crucial to understand the liquid–vapor behavior for enhancing the LHP performance by investigating the fundamental heat and mass transfer in the wick with phase-change. However, previous modeling studies only considered a single-phase ﬂow or complete saturation in the wick, and the capillary effect on the ﬂuid states was rarely taken into account. The present work developed two mathematical models based on the assumptions of saturated and unsaturated wicks. The ﬂuid states were analyzed at the liquid–vapor interface under the consideration of the capillary effect, and a pore-scale evaporation model was applied to study the phase change behavior and interfacial heat and mass transfer. The relative permeability was introduced to describe the two-phase ﬂow in the porous wick, and the capillary force was modeled as a function of the local saturation in the two-phase region. The temperature results calculated by the models were compared with the experimental results, and the assumption that the vapor penetration leads to deterioration of evaporator performance at high heat ﬂux was validated. Vapor blanket thickness can be estimated through the saturation proﬁle, which provides a simple and effective method. It was also found that the capillary number ω was the key factor affecting the thickness of the vapor blanket. The greater the ω , the faster the vapor blanket thickness increases.


Introduction
A loop heat pipe (LHP) is a highly efficient two-phase heat transfer device which utilizes the evaporation and condensation of a working fluid to transfer heat [1].As illustrated in Figure 1, it consists of an evaporator coupled with a compensation chamber (CC), a capillary porous structure (wick) integrated in the evaporator, a condenser, and transport lines (connecting the evaporator and the condenser).When the evaporator is heated, the working fluid in the wick starts to evaporate, and an evaporating meniscus is formed; hence, a capillary force develops for circulating the condensate working fluid back to the liquid-vapor interface from the condenser [2].Due to the magnitude of latent heat and the capillary effect, the LHP can realize both high heat transfer capacity and long distance.
LHPs have engaged the attention of scientists and engineers since the first prototype was unveiled in 1972 [3], because they are capable of functioning at any orientation with respect to gravity and transferring heat over sufficiently large distances [4].They have been applied in multitudinous fields [5], such as spacecraft thermal control [6,7], infrared detection systems [8,9], aircraft anti-icing [10], battery thermal management [11][12][13] electronics cooling [14][15][16].In these applications of LHP, the evaporator plays an important role, because it is the part in direct contact with the heat source among the components of LHP, and its heat transfer capability significantly affects the heat transfer capacity of the LHP.The heat transfer capability of the evaporator can be characterized by an evaporative heat transfer coefficient, which can be obtained from the evaporator thermal resistance and evaporation area.The thermal resistance from the evaporator wall to the liquid-vapor interface is the major influencing factor of the evaporative heat transfer coefficient [17].In order to enhance the heat transfer capacity of LHP, it is essential to understand the liquid-vapor interface behavior in the wick.LHPs have engaged the attention of scientists and engin was unveiled in 1972 [3], because they are capable of functio respect to gravity and transferring heat over sufficiently la been applied in multitudinous fields [5], such as spacecraft t detection systems [8,9], aircraft anti-icing [10], battery therm electronics cooling [14][15][16].In these applications of LHP, t portant role, because it is the part in direct contact with the ponents of LHP, and its heat transfer capability significantly pacity of the LHP.The heat transfer capability of the evapor an evaporative heat transfer coefficient, which can be obtain mal resistance and evaporation area.The thermal resistance the liquid-vapor interface is the major influencing factor of coefficient [17].In order to enhance the heat transfer capac understand the liquid-vapor interface behavior in the wick.
Visualization is a reliable method to observe the posi phase region.Some researchers have developed visualizatio Visualization is a reliable method to observe the position and outline of the twophase region.Some researchers have developed visualization studies on the evaporator to dissect the liquid-vapor interface behavior in the porous wick.Liao and Zhao [18] conducted a visual experiment to study the two-phase zone behavior in a two-dimensional porous structure and used a high-speed video imaging system to record the photographs.The results illustrated that a co-existing zone of bubbles and liquid existed in the porous structure adjacent to the heated surface, this two-phase zone expanded with the increase of the heat flux, and a stable liquid film formed with a further increase of heat flux, which resulted in a rapid rise of the temperatures both in the heating surface and at the vapor exit.Odagiri and Nishikawara [19] conducted an observation experiment with a microscopic infrared camera and a microscope to investigate the characteristics of the thermo-fluid behavior on the surface of the porous wick in the LHP system.They reported that there are three types of thermo-fluid behaviors on the porous wick surface, corresponding to different heat fluxes.At low or moderate heat flux, evaporation occurs on the surface of the liquid bridge, but with the increase of the heat flux, the liquid bridges disappear, and evaporation occurs at the liquid-vapor interface that is formed in the porous wick.Meanwhile, a vapor blanket forms adjacent to the heating zone, which grows quickly as the heat flux increases.Other visualization tools have also been adopted to the study of phase change behavior in the porous media, such as computed tomography [20], nuclear tracer imaging [21], nuclear magnetic resonance [22], particle image velocimetry [23], and neutron radiography [24].However, all the above-mentioned methods for observing the liquid-vapor interface behavior in the porous media are limited in spatial or temporal resolution or in sample size, or require extremely expensive equipment.For that reason, a number of theoretical studies have been conducted to analyze the liquid-vapor interface behavior in porous media.Khrustalev and Faghri [25] developed a one-dimensional mathematical model to investigate the heat and mass transfer during evaporation in a porous structure by means of the film evaporation theory.They used the extended Kelvin's equation to determine the vapor saturation pressure and proved that the vapor blanket could be formed with the increasing of the heat flux, which would impair the performance of the evaporator.Lin et al. [26] presumed that there are three fluid regions in the porous wick during evaporation: a vapor blanket, a two-phase region, and a saturated liquid region.Based on this assumption, they developed a steady-state mathematical model to numerically analyze the effects of pore-size distribution on the heat transfer performance.The calculations indicated that the thermal resistance of the bidisperse wick could be held to a lower level compared to the monoporous wick, and they attributed this phenomenon to the thinner vapor blanket in the bidisperse wick.Ranjan et al. [27] proposed a numerical model for the evaporating liquid meniscus to study the steady evaporation process under the assumption of a static liquid-vapor interface shape.They used the Surface Evolver to solve the Young-Laplace equation to obtain the liquid-vapor interface morphologies and evaporative performance in four distinct wick geometries.Based on the results, it can be concluded that the evaporation heat transfer from the liquid-vapor interface decreases with increasing contact angle but increases with increasing superheat.Odagiri and Nagano [28] presented a steady-state mathematical model to study the liquid-vapor interface behavior in the evaporator, which considered the variations of the evaporator heat transfer coefficient.The quantitative relationship between the performance of the LHP and liquid-vapor interface behaviors was given based on this model.
Nevertheless, in the literatures mentioned above, the impact of the pressure jump caused by the capillary effect on the fluid states (e.g., pressure and temperature) was not taken into account while modeling the evaporator in the case of a saturated wick.On the other hand, in the case of the modeling of an unsaturated wick, the fluid flow in the porous interior was simply described by Darcy's law, without considering the two-phase flow.Although it is claimed in the references that the two-phase length is quite short, it is worth noting that the two-phase flow directly affects the momentum equation, and the capillary force can only be realized in the two-phase region.Furthermore, it was always assumed that the capillary force was equal to the maximum capillary force provided by the meniscus; the interactions between the wick saturation and the actual capillary force were neglected.
As mentioned above, the saturation state of the wick has an important impact on the heat transfer performance of the LHP.It is essential to study the heat and mass transfer during evaporation in the wick, and the main objective for the presented work is to find the relationship between the vapor blanket thickness and the heat fluxes, which can help us to set a safe boundary for the applied heat flux of the LHP to prevent the drying out.Therefore, a one-dimensional mathematical model was proposed; it provides a novel theoretical method, which can simply and effectively evaluate the vapor blanket thickness in the wick.Unlike previous models, this model adequately considered the capillary effect at the liquid-vapor interface and the two-phase flow in the porous wick, which is the novelty of the present work.

Experimental Description
This section briefly summarizes the experimental results from our previous studies using three new refrigerants which are, respectively, R245fa, R1234yf and R1234ze(E), as working fluids in a LHP with a flat-disk evaporator.Detailed experimental results have been presented in Ref. [29].Nevertheless, here we would like to introduce a notable phenomenon of the evaporator heat transfer coefficient transition with different heat fluxes.Some details of the experimental rig are given for better understanding.
The main parameters of the experimental LHP for this study are given in Table 1. Figure 2 illustrates the schematic of the experimental apparatus; the heater was provided by a regulated DC power supply, and the voltage and the current were recorded to calculate the applied heat transfer rate by digital meters that had an uncertainty value of less than ±0.5%.Twelve Pt1000 temperature sensors with an accuracy ±0.1 • C were used for monitoring the key locations temperature of the LHP in conjunction with a data acquisition system.The tee with a valve was installed near the inlet of the CC for vacuum pumping and charging, and the charging ratio for the three refrigerants was 65-70%.The condenser was connected with the cooling plate, the temperature of which was controlled by the circulating cooling system, and the cooling system temperature was set as 20 ± 0.1 • C.There were three temperature sensors evenly arranged on the condenser wall, and the mean value of these temperatures was taken as the temperature of the condenser.Similarly, four temperature sensors were symmetrically arranged on the upper and lower surfaces of the evaporator, and the mean value of these temperatures was considered as the evaporator temperature.A ceramic wick with vapor grooves on the surface was applied in the evaporator; the heat was transferred to the wick surface, vaporizing the working fluid.Then the vapor vented through the grooves into the vapor line, and the vapor temperature T v was recorded by a temperature sensor on the vapor line located 7 mm from the evaporator exit.The thermal resistance is usually the main indicator for a LHP's performance, which is defined as the ratio of the temperature difference between the evaporator and the condenser to the heat load.However, in most application scenarios, the heat sink often has sufficient cooling capacity, so one may pay more attention to the heat extraction capacity of the LHP or, we can say, to the heat transfer efficiency of the evaporator.The heat transfer coefficient can be defined to characterize the heat transfer efficiency of the evaporator as follows where q is the applied heat flux, and T e and T v are the temperatures of the evaporator and the vapor, respectively.Figure 3 illustrates the variation of the evaporator heat transfer coefficient versus the applied heat flux.It is shown that with an increase of the heat flux, the evaporator heat transfer coefficient increased, then reached a peak value, and dropped afterwards.This peak value will be referred to as the peak heat transfer coefficient h pe while the corresponding heat flux at which the peak heat transfer occurs will be referred to as the peak heat flux q pe .The h pe of the evaporator with R245fa, R1234yf and R1234ze are 10,792.56W/m 2 K, 7676.57W/m 2 K and 7753.65 W/m 2 K, and the corresponding q pe are 75 kW/m 2 , 25 kW/m 2 and 25 kW/m 2 , respectively.Here is a hypothetical explanation: when the heat flux is less than q pe , the porous wick is saturated by the liquid, and evaporation uniformly takes place at the wick-fin interface; the evaporator heat transfer coefficient increases with the imposed heat flux, which may be attributed to the gradual thinning of the evaporating liquid film near the heating surface.With a further increase of the heat flux, vapor zones form in the pores adjacent to the heating surface.As these local vapor zones keep receding into the porous structure, a vapor blanket forms along the heating surface once the heat flux exceeds q pe , which causes the evaporator heat transfer coefficient to drop and the wall temperature and vapor temperature to rise.To prove this hypothesis, two mathematical models were presented for both saturated and unsaturated wick in Section 3.

Working fluid
R245fa/R1234yf  The thermal resistance is usually the main indicator for a LHP's performance, wh is defined as the ratio of the temperature difference between the evaporator and the c denser to the heat load.However, in most application scenarios, the heat sink often sufficient cooling capacity, so one may pay more attention to the heat extraction capa of the LHP or, we can say, to the heat transfer efficiency of the evaporator.The heat tra fer coefficient can be defined to characterize the heat transfer efficiency of the evapora as follows where  is the applied heat flux, and   and   are the temperatures of the evapora and the vapor, respectively.
Figure 3 illustrates the variation of the evaporator heat transfer coefficient versus applied heat flux.It is shown that with an increase of the heat flux, the evaporator h transfer coefficient increased, then reached a peak value, and dropped afterwards.T peak value will be referred to as the peak heat transfer coefficient ℎ  while the co sponding heat flux at which the peak heat transfer occurs will be referred to as the p heat flux   .The ℎ  of the evaporator with R245fa, R1234yf and R1234ze are 10,792 W/m 2 K, 7676.57W/m 2 K and 7753.65 W/m 2 K, and the corresponding   are 75 kW/ 25 kW/m 2 and 25 kW/m 2 , respectively.Here is a hypothetical explanation: when heat flux is less than   , the porous wick is saturated by the liquid, and evaporation u formly takes place at the wick-fin interface; the evaporator heat transfer coefficient creases with the imposed heat flux, which may be attributed to the gradual thinning the evaporating liquid film near the heating surface.With a further increase of the h flux, vapor zones form in the pores adjacent to the heating surface.As these local va zones keep receding into the porous structure, a vapor blanket forms along the heat surface once the heat flux exceeds   , which causes the evaporator heat transfer coe cient to drop and the wall temperature and vapor temperature to rise.To prove this pothesis, two mathematical models were presented for both saturated and unsatura wick in Section 3.

Physical Model and Mathematical Formulation
A schematic of the computational domain for the evaporator is presented in Figure 4.The heat load applied on the evaporator wall traverses through conduction fins to wick the bottom surface, as illustrated in Figure 4, placed at y = δ along the Y coordinate.The same assumption was applied as in most of the literature, which neglected the parasitic heat flux between evaporator and ambient.Consequently, the applied heat load can be divided into two parts: one part of it is spent in heating the subcooled liquid to reach saturation temperature, and the other part is transferred through the liquid-vapor interface as evaporation occurs.A fully saturated wick was assumed in most previous analytical and numerical models; however, there are some critical circumstances that may be encountered, i.e., an extremely high heat flux, which is of more interest for electronics applications, and which will change the wick saturation status and lead to another operational regime where a vapor blanket exists inside the capillary wick.

Physical Model and Mathematical Formulation
A schematic of the computational domain for the evaporator is presented in Figure 4.The heat load applied on the evaporator wall traverses through conduction fins to wick the bottom surface, as illustrated in Figure 4, placed at y = δ along the Y coordinate.The same assumption was applied as in most of the literature, which neglected the parasitic heat flux between evaporator and ambient.Consequently, the applied heat load can be divided into two parts: one part of it is spent in heating the subcooled liquid to reach saturation temperature, and the other part is transferred through the liquid-vapor interface as evaporation occurs.A fully saturated wick was assumed in most previous analytical and numerical models; however, there are some critical circumstances that may be encountered, i.e., an extremely high heat flux, which is of more interest for electronics applications, and which will change the wick saturation status and lead to another operational regime where a vapor blanket exists inside the capillary wick.Vapor penetrates into the wick, which blocks liquid feeding and causes a partial dry-out spot in the wick at high heat flux.
The working state of a capillary wick can be different in the light of magnitude of applied heat flux.In the case of a low or moderate heat flux (Figure 4a), the wick is completely saturated with liquid, evaporation takes place at the wick-groove interface, and vapor flows through the groove into a vapor line.As a result, the liquid is drawn through the wick by capillary force to the liquid-vapor interface to preserve the mass balance.At high heat flux (Figure 4b), vapor may penetrate into wick structure and the liquid-vapor interface recede into the wick, which will cause a vapor blanket to occur and the vapor to flow through the wick thickness to the groove.In this case, the liquid feeding rate decreases and the fluid flow pressure drop in the wick increases significantly, which may lead to dry-out in the wick.Depending on the saturation state of the wick, the thermodynamic and hydrodynamic behaviors of the fluid in the porous wick can be highly distinct We shall treat the mathematical model separately in two subsections.

The Case of Low or Moderate Heat Flux
In this section, the energy equation for liquid flow in the saturated wick was solved and the temperature distribution along the thickness direction of the wick was obtained Also, the fluid states were analyzed at the liquid-vapor interface under the consideration of the capillary effect.The evaporator wall temperature was estimated based on the ther- The working state of a capillary wick can be different in the light of magnitude of applied heat flux.In the case of a low or moderate heat flux (Figure 4a), the wick is completely saturated with liquid, evaporation takes place at the wick-groove interface, and vapor flows through the groove into a vapor line.As a result, the liquid is drawn through the wick by capillary force to the liquid-vapor interface to preserve the mass balance.At high heat flux (Figure 4b), vapor may penetrate into wick structure and the liquid-vapor interface recede into the wick, which will cause a vapor blanket to occur and the vapor to flow through the wick thickness to the groove.In this case, the liquid feeding rate decreases and the fluid flow pressure drop in the wick increases significantly, which may lead to dry-out in the wick.Depending on the saturation state of the wick, the thermodynamic and hydrodynamic behaviors of the fluid in the porous wick can be highly distinct.We shall treat the mathematical model separately in two subsections.

The Case of Low or Moderate Heat Flux
In this section, the energy equation for liquid flow in the saturated wick was solved, and the temperature distribution along the thickness direction of the wick was obtained.Also, the fluid states were analyzed at the liquid-vapor interface under the consideration of the capillary effect.The evaporator wall temperature was estimated based on the thermal resistance from the wick bottom surface to the evaporator wall.

Heat Transfer Analysis
The energy equation in the saturated wick with fluid flow can be expressed as In order to simplify the energy equation, some assumptions are introduced: (i) a homogeneous, isotropic porous wick; (ii) one-dimensional, steady, incompressible flow; (iii) local thermal equilibrium between the solid and fluid phases; (iv) one-dimensional, steady heat transfer; (v) gravity and heat radiation are neglected.
Considering the computational domain as illustrated in Figure 2a, the energy equation for liquid flow in the wick for an elemental thickness ∆y is given by where u l is the liquid velocity in the wick, ρ l is the liquid density, C p,l is the specific heat of liquid, and λ e f f l,w is the effective thermal conductivity of the wick saturated with liquid.Based on assumption (i), λ e f f l,w can be represented based on volumetric averages of two phases in parallel, as in [26] where λ w is the thermal conductivity of the porous wick, λ l is the thermal conductivity of the liquid, and ε is wick porosity.The boundary conditions for the computational domain are: Solving Equation (3) with the above boundary conditions yields the temperature distribution in the wick as where q cont is the heat flux at the contact between the evaporator plate and wick, and the liquid velocity u l is related to mass flow rate of fluid .m, which can be expressed as Thus, Equation ( 5) can be deformed into To determine the mass flow-rate value, it is necessary to analyze the heat transfer process in the wick.The schematic of the cross-sectional view of the evaporator-CC structure and details of heat fluxes was presented in Figure 5.The applied heat load Q in traveled through the evaporator plate into the saturated wick; one part of it transferred to heat the liquid in the CC to reach the temperature T δ and balanced subcooling of condensate liquid.The other part transferred to the wick-groove interface, and vapor left the interface with a certain mass flow rate .m carrying the latent heat H f g .The energy balance is established in the evaporator, expressed as With the assumption that the heat exchange between evaporator case and ambie can be ignored, the contact heat flux   at the evaporator plate-wick interface co pletely enters C, and compensates the liquid sensible heat from   to   , which is a referred as the heat leak from wick to CC.The energy balance in the CC is given by    =   +   = ̇ , (  −   ) + ̇ , (  −   ) = ̇ , (  −   ) ( Combining Equations ( 9) and (10), the mass flow rate can be estimated as The wick is saturated with liquid at a low or moderate heat flux, and the liquid eva orates at the wick-groove interface.Here the vapor temperature can be expressed as where Δ  is the saturation temperature difference between the liquid and vapor at t wick-groove interface, which is caused by the thermal resistance across the evaporati meniscus, also known as interfacial thermal resistance  − , which is given by Her Knudsen equation [30] where coefficient α is the so-called accommodation coefficient that shows the portion liquid molecules going into the vapor surface and absorbed by this surface, and for t model the coefficient is equal to 0.5.However, it should be noted that the order of mag tude of α can be as high 10 3 in some particular cases. is the gas constant,  rep With the assumption that the heat exchange between evaporator case and ambient can be ignored, the contact heat flux q cont at the evaporator plate-wick interface completely enters C, and compensates the liquid sensible heat from T in to T δ , which is also referred as the heat leak from wick to CC.The energy balance in the CC is given by Combining Equations ( 9) and (10), the mass flow rate can be estimated as The wick is saturated with liquid at a low or moderate heat flux, and the liquid evaporates at the wick-groove interface.Here the vapor temperature can be expressed as where ∆T lv is the saturation temperature difference between the liquid and vapor at the wick-groove interface, which is caused by the thermal resistance across the evaporating meniscus, also known as interfacial thermal resistance R l−v , which is given by Hertz-Knudsen equation [30] where coefficient α is the so-called accommodation coefficient that shows the portion of liquid molecules going into the vapor surface and absorbed by this surface, and for this model the coefficient is equal to 0.5.However, it should be noted that the order of magnitude of α can be as high 10 3 in some particular cases.R g is the gas constant, P v represents the saturation pressure corresponding to the vapor temperature, and v v and v l are the vapor and liquid specific volume, respectively.Accordingly, applying the energy conservation equation at the evaporation interface, Equation ( 14) can be obtained .

Thermodynamics Analysis
The diagram of the LHP working cycle with respect to the saturation line of the working fluid in P-T coordinates is presented in Figure 6.The liquid feeding process from CC to wick is fulfilled due to the pressure difference between the liquid and vapor phases at the evaporation interface, which is so-called capillary force P c , and be determined by Young-Laplace equation where σ is the surface tension, θ is the contact angle, and R e f f is the average effective radius of wick pores.

Thermodynamics Analysis
The diagram of the LHP working cycle with respect to the saturation lin ing fluid in P-T coordinates is presented in Figure 6.The liquid feeding pro to wick is fulfilled due to the pressure difference between the liquid and va the evaporation interface, which is so-called capillary force   , and be d Young-Laplace equation where  is the surface tension,  is the contact angle, and   is the ave radius of wick pores.As we can see in Figure 6, the capillary force can be sustained due to difference between the vapor and liquid at the interface.The working fluid late from the evaporator to the CC, which is driven by the pressure differ these two parts.The pressure drop of the working fluid flowing from the the CC is usually defined as the external (of the wick) pressure drop, ba pressure difference between vapor and CC ∆  =   −   Similarly, the pressure drop of the working fluid flowing through the w as the internal pressure drop, which is represented by the pressure differenc and the liquid phase at the interface.
∆  =   −   As we can see in Figure 6, the capillary force can be sustained due to the pressure difference between the vapor and liquid at the interface.The working fluid starts to circulate from the evaporator to the CC, which is driven by the pressure difference between these two parts.The pressure drop of the working fluid flowing from the evaporator to the CC is usually defined as the external (of the wick) pressure drop, balanced by the pressure difference between vapor and CC ∆P EX = P v − P cc (16) Similarly, the pressure drop of the working fluid flowing through the wick is defined as the internal pressure drop, which is represented by the pressure difference between CC and the liquid phase at the interface.
Also, the internal pressure drop can be determined by Darcy's law.
where µ l is the liquid viscosity, and K w is the permeability, which is related to the pore diameter d and porosity ε of the porous wick [31].
In the steady-state operation of LHP, the sum of external and internal pressure drops is equal to actual capillary force The actual contact angle of the meniscus θ can be estimated Equation (20).The capillary force reaches maximum value, which is known as the capillary limit, when the contact angle is zero.The total pressure drops not exceeding the capillary limit is always regarded as one criterion of a functioning LHP, which is expressed as the flowing inequality Here, due to the capillary force, the pressure discontinuity must exist at the liquidvapor interface, which introduces a question as to the thermodynamic states of the liquid and vapor, which are in thermodynamic equilibrium locally.On the one hand, if we take the liquid phase as a saturated state, the vapor which is in contact with liquid at a higher pressure due to the surface tension would be a subcooled vapor.Due to the extremely unstable nature of subcooled vapor, it is hardly possible in a porous system.On the other hand, if we take the vapor phase as a saturated state, the liquid which is in contact with vapor at a lower pressure due to the surface tension would thus be a superheated liquid, which is considered to be an unconditionally stable state [32,33].
The vapor and liquid are in the thermodynamic equilibrium state at the liquid-vapor interface.Based on the equilibrium condition, the thermodynamic potential energies of two phases at the liquid-vapor interface are equal, which can be expressed by the conservation of Gibbs free energy under isothermal condition [34] where the Gibbs free energy of vapor is given by where G 0 and P 0 are the Gibbs free energy and pressure at saturation state.
Applying the ideal gas assumption to the vapor at the interface, here introducing the ideal gas state equation in Equation ( 23), the following equation can be derived where R g is the gas constant, which is equal to the general gas constant divided by the molar mass of the gas.
Similarly, the Gibbs free energy of superheated liquid can be written as Assuming the change of liquid specific volume can be neglected for the pressure difference, Equation (25) becomes Substituting Equations ( 24) and ( 26) into Equation ( 22) can obtain and vapor pressure is given by following equation, which is somewhat familiar with the Kelvin equation where P 0 is determined by the Clausius-Clapeyron equation Thus, the vapor temperature T v can be calculated based on Equations ( 12)-( 14) and Equations ( 18), ( 28) and ( 29), assuming that the vapor leaves the interface at temperature T v , which is the same as the temperature at the tip of the wick fin, as shown in Figure 7, which demonstrated the thermal resistance path from wick fin to evaporator plate.Interference fit and welding can get close contact between the wick fin and the evaporator plate, which results in good thermal connection via a contact heat transfer coefficient of 10,000 W/m 2 K, as applied in Ref. [35].The thermal conductivity of the evaporator plate (stainless steel-316) is 16.3 W/mK.Consequently, the evaporator temperature can be estimated by the heat conduction through the wick fin and the evaporator plate where R HL is the heat leak thermal resistance from evaporator shell to CC.

The Case of High Heat Flux
In the case of low or moderate heat flux, the wick is considered as completely saturated with liquid and evaporation takes place at the wick-groove interface, which is very close to the contact between the evaporator and the wick fin.With the increase of heat flux, the evaporation rate of working fluid increases, and the flow resistance of liquid in the wick increases, which will reduce the replenishment rate of CC.With a further increase of heat flux, there could be formation of a vapor blanket, which thickness is actually non-uniform and increases with the heat flux.This vapor blanket will make the liquid-vapor interface move away from the tip of the wick fin, causing an increase of the thermal resistance in the wick, which will furtherly lead to a sharp rise of the evaporator temperature.
In this section, a pore-scale evaporation model based on extended meniscus theory has been developed to investigate the heat transfer process at the liquid-vapor interface in the wick.With the postulate of a flat liquid-vapor interface and a uniform vapor blanket in the wick, the vapor blanket thickness can be estimated.Despite the fact that the vapor blanket is usually non-uniform due to the pores' distribution and surface roughness, as an approximation an equivalent uniform thickness was considered for analyzing the effect of the liquid-vapor interface location on the evaporation heat transfer coefficient.
where   is the heat leak thermal resistance from evapora

The Case of High Heat Flux
In the case of low or moderate heat flux, the wick is co rated with liquid and evaporation takes place at the wick-gr close to the contact between the evaporator and the wick f flux, the evaporation rate of working fluid increases, and th the wick increases, which will reduce the replenishment ra crease of heat flux, there could be formation of a vapor blank non-uniform and increases with the heat flux.This vapor b vapor interface move away from the tip of the wick fin, caus resistance in the wick, which will furtherly lead to a sharp r ature.
In this section, a pore-scale evaporation model based o has been developed to investigate the heat transfer process in the wick.With the postulate of a flat liquid-vapor interface in the wick, the vapor blanket thickness can be estimated.D blanket is usually non-uniform due to the pores' distributi an approximation an equivalent uniform thickness was consi of the liquid-vapor interface location on the evaporation he

Film Evaporation in Pore Scale
As shown in Figure 8a, a typical thin-film region consists of three distinct regions [2,36]: (I) the equilibrium film region, also referred to as the absorbed non-evaporating region, which is influenced by disjoining pressure (due to intermolecular forces between the liquid molecules and the solid skeleton).(II) The evaporating film region, which is influenced by both disjoining pressure and capillary pressure; the film thickness in this region increases and disjoining pressure decreases, causing evaporation and heat transfer.(III) The intrinsic meniscus region, where the disjoining pressure produces a negligible effect on evaporation, and the thick film introduces large thermal resistance.

PEER REVIEW
12 of 19

Film Evaporation in Pore Scale
As shown in Figure 8a, a typical thin-film region consists of three distinct regions [2,36]: (I) the equilibrium film region, also referred to as the absorbed non-evaporating region, which is influenced by disjoining pressure (due to intermolecular forces between the liquid molecules and the solid skeleton).(II) The evaporating film region, which is influenced by both disjoining pressure and capillary pressure; the film thickness in this region increases and disjoining pressure decreases, causing evaporation and heat transfer.(III) The intrinsic meniscus region, where the disjoining pressure produces a negligible effect on evaporation, and the thick film introduces large thermal resistance.For the extended meniscus, the Young-Laplace equation as explained in Equation ( 15) should be extended as well.An additional term, which represents the pressure difference between the liquid close to the solid skeleton   and the liquid close to the interface For the extended meniscus, the Young-Laplace equation as explained in Equation ( 15) should be extended as well.An additional term, which represents the pressure difference between the liquid close to the solid skeleton P δ and the liquid close to the interface P l , also known as the disjoining pressure P d is included, yielding the augmented Young-Laplace equation [37][38][39] P v − P δ = P c − P d Here, P c is the pressure difference between the vapor and liquid close to the interface, also known as the capillary pressure, which is related with the interfacial tension σ and the curvature of the meniscus and the meniscus curvature is the function of the film thickness given by where δ l and δ l are the first and the second derivative of the film thickness.To simplify the calculation, it is assumed that the curvature of the meniscus is 2/R men .The disjoining pressure represents the effect of intermolecular forces between the liquid and solid.For nonpolar fluids, the disjoining pressure can be represented as a function of dispersion constant A (also known as the Hamaker constant) and the film thickness, gives as [2] with A = 3.21 × 10 −21 J in the present analysis.As illustrated in Figure 8b, the heat transfer from wick skeleton to vapor was idealized to a problem of heat conduction across the liquid film with nominal thickness δ l and evaporation from the liquid-vapor interface situated in a cylindrical pore with an effective radius R e f f .The conduction heat flux across the liquid film is obtained as where T w is the temperature of wick skeleton, T lv is the interfacial temperature at liquidvapor interface which is affected by the disjoining pressure P d and capillary force P c , and also depends on the interfacial thermal resistance.The heat flux due to evaporation at the liquid-vapor interface can be obtained from [40] The interfacial pressure at the liquid-vapor interface P lv (η) in Equation ( 36) can be obtained by the extended Kelvin equation [32] P lv (η) = P s (T lv )exp P lv (η) − P s (T lv ) + P d (η) − P c ρ l R g T lv (η) (37) δ l is the nominal thickness of the liquid film.
Energies 2023, 16, 2088 14 of 21 Under steady-state conditions, the heat conduction flux across liquid film is balanced by evaporation heat flux; thus, the following expression can be derived from Equations ( 35) and ( 36) The interfacial temperature T lv (η) and pressure P lv (η) at the liquid-vapor interface can be obtained by solving Equations ( 37) and (38) simultaneously for given values of T v , T w and δ l (η).
As indicated in Equation ( 34), the disjoining pressure increases with the decrease of liquid film thickness; therefore, a liquid film will become extremely thin near the top of the pore (η → 0), and the disjoining pressure in this region will correspondingly become rather high.
The evaporation heat flux from each pore can be determined by the equivalent evaporation heat transfer coefficient h lv and temperature difference between the vapor and the interface where h lv can be derived from Equation ( 13) Thus, the evaporation heat transfer from each pore is obtained by integrating the evaporation heat fluxes of the open pore area.
where t v is the top of the meniscus on the η coordinate axis, which represents the equivalent thickness of the vapor blanket in the wick, as indicated in Figure 9. Equating the Equations ( 35) and (36), and eliminating the interface temperature T lv , the evaporation heat transfer from each pore can be expressed as with the assumption of the uniform distribution of the pores on the wick surface, and the evaporation heat transfer from each pore being equal, then total evaporation heat transfer can be estimated by where N is the number of pores on the wick surface, which can be determined by where A e f f is the across area of a pore.
where N is the number of pores on the wick surface, which can be dete where   is the across area of a pore.

Determination the Thickness of Liquid Film and Vapor Blanket
The change of local liquid film thickness   along the meniscus is a coordinate, as shown in Figure 9, which is geometrically related to the dius   and meniscus radius   given by The nominal film thickness of the heat conduction path can be exp

𝜂
As indicated in Figure 9, when the vapor blanket has a thicknes pressure can be estimated by Equations ( 31)-( 33), ( 41) Heat transfer in a pore.

Determination the Thickness of Liquid Film and Vapor Blanket
The change of local liquid film thickness δ l along the meniscus is a function of the η coordinate, as shown in Figure 9, which is geometrically related to the effective pore radius R e f f and meniscus radius R men given by The nominal film thickness of the heat conduction path can be expressed as As indicated in Figure 9, when the vapor blanket has a thickness of t v , the vapor pressure can be estimated by Equations ( 31)-( 33), ( 41) When there is two-phase flow in a porous media, the interaction between the two phases needs to be considered.The feeding of the liquid phase will be reduced by the presence of the vapor phase in pores, and the single-phase flow is not only influenced by the permeability of the porous media but also by the relative permeability of the fluid.Relative permeabilities are used as parameters in the fractional flow theory, which was introduced by Leverett [41].Whereas Darcy's law for single-phase flow has received wide validation for numerous types of porous media, in the porous media system with a two-phase flow Darcy's law should be modified by relative permeability to calculate the pressure drop of a single-phase flow in the porous media system with that of a two-phase flow in the The appropriate boundary condition at the top of the wick is s = 1 at η = 0. Therefore Equation (56) can be integrated numerically for the saturation profile given the correlations of Equations ( 50), ( 51) and (53).

Results and Discussion
The experimental temperature of the vapor and the evaporator was recorded and compared with the results from the two numerical models.The programs of the model were implemented in Matlab-R2018a, and a calling function was embedded in the program to realize the correlation between physical properties and temperature of the three working fluids through the NIST Refprop-9.1 database.
As illustrated in Figure 10, under the controlled condition that the evaporator temperature is no more than 350 K, the corresponding maximum heat flux for R245fa, R1234yf and R1234ze is 150 kW/m 2 , 100 kW/m 2 , 100 kW/m 2 , respectively.The calculated values, without considering the vapor penetration, showed good agreement with the experimental results at lower heat flux, yet the difference between the experimental and calculated temperatures keeps increasing with increase of heat flux; the obvious bifurcate points occurred once the heat flux exceed a certain specific value, which for R245fa, R1234yf and R1234ze is 75 kW/m 2 , 25 kW/m 2 and 25 kW/m 2 , respectively.The model without considering the vapor penetration took the wick as a liquid-saturated state, and the liquid-vapor interface temperature was recognized as the vapor temperature at the evaporator exit; the temperature difference between the wick skeleton and the vapor is small enough to neglect when the wick is saturated.However, with the increase of heat flux, the saturation state of the wick would change, the evaporation interface cannot be maintained at the wick-groove interface, and the liquid-vapor interface would recede into the wick for a higher replenishment rate.Thus, vapor leaving the evaporation interface will travel through the part of the dry area in the wick to the evaporator exit, and the wick skeleton will continuously heat the vapor through the path, resulting in the vapor being superheated, and causing an underestimation of the vapor temperature in the saturated model.Furthermore, the evaporator temperature will be underestimated, as inferred in Equation (30) as well.
The calculated results considering the vapor penetration coincide well with the experimental results; the model fully considered the variation of the wick saturation along with heat flux.A pore-scale heat transfer analysis was applied, and the superheat between the vapor and wick skeleton was analogous to a one-dimensional conduction heat transfer across the evaporation meniscus; the process by which vapor was heated by the dry area of the wick was considered to determine the temperature at the evaporator exit.Thus, the temperature of vapor and evaporator calculated by the second model is closer to the experimental results.From the comparison between the results of the two models and experimental results, it can be found that when the heat flux exceeds a certain value, the assumption of a saturated wick will no longer be reasonable, which proves the hypothetical explanation proposed in Section 2.
Figure 11a shows the liquid saturation profile along the wick thickness predicted by the second model (as calculated from Equation ( 56)).As a demonstration case, taking the heat flux equal to 100 kW/m 2 , the saturation changes from 0 to 1 along the wick thickness, and closer to the upper surface of the wick the saturation is closer to 1. Contrarily, the saturation is closer to 0 at the location closer to lower surface.Meanwhile, near the upper and lower surfaces the saturation changing rate appears to be accelerating.This is affected by the evaporation-replenishment process; the upper surface has a higher replenishment rate, while the lower surface has a higher evaporation rate.Figure 11b approximately characterizes the wick as three regions based on the saturation value.These regions are a liquid region (s = 1), two-phase region (s = 0~1), and vapor region (s = 0), respectively.It should be noted that the real distribution will be more complicated due to the complexity of the internal structure of porous media.Nevertheless, the aforementioned calculation results of the temperature have shown that it is a simple and effective method to estimate the vapor blanket thickness through the saturation distribution curve.As shown in Figure 11a, the equivalent vapor blanket thickness for R245fa, R1234yf, and R1234ze is 0.96 mm, 1.71 mm, and 1.57 mm at 100 kW/m 2 , respectively.The calculated results considering the vapor penetration coincide well with the experimental results; the model fully considered the variation of the wick saturation along with heat flux.A pore-scale heat transfer analysis was applied, and the superheat between the vapor and wick skeleton was analogous to a one-dimensional conduction heat transfer should be noted that the real distribution will be more complicated due to the complexity of the internal structure of porous media.Nevertheless, the aforementioned calculation results of the temperature have shown that it is a simple and effective method to estimate the vapor blanket thickness through the saturation distribution curve.As shown in Figure 11a, the equivalent vapor blanket thickness for R245fa, R1234yf, and R1234ze is 0.96 mm, 1.71 mm, and 1.57 mm at 100 kW/m 2 , respectively.The variation of the equivalent vapor blanket thickness with the heat flux for all the fluids is shown in Figure 12.The equivalent vapor thickness increases almost linearly with the increase of the heat flux; the growing rate for R245fa is the lowest, while R1234yf and R1234ze have a similar growing rate.This phenomenon is determined by the dimensionless parameter  referred to in Equation ( 53), which is also called the capillary number in fluid dynamics terminology.As the  increases, the vapor fraction enlarges; thus, a larger capillary pressure is established for a higher  in order to keep the mass balance between the replenishment and evaporation.The similar conclusion was also obtained by Wang [45].The variation of the equivalent vapor blanket thickness with the heat flux for all the fluids is shown in Figure 12.The equivalent vapor thickness increases almost linearly with the increase of the heat flux; the growing rate for R245fa is the lowest, while R1234yf and R1234ze have a similar growing rate.This phenomenon is determined by the dimensionless parameter ω referred to in Equation ( 53), which is also called the capillary number in fluid dynamics terminology.As the ω increases, the vapor fraction enlarges; thus, a larger capillary pressure is established for a higher ω in order to keep the mass balance between the replenishment and evaporation.The similar conclusion was also obtained by Wang [45].

Conclusions
In order to analyze the evaporator operational characteristics, two mathematical models have been developed based on the assumptions of saturated and unsaturated wicks.The fluid states were analyzed at the liquid-vapor interface under the consideration of the capillary effect.A pore-scale evaporation model was applied to study the phase change behavior and interfacial heat and mass transfer, the relative permeability was introduced to describe the two-phase flow in the porous wick, and the capillary force was modeled as a function of the local saturation in the two-phase region.The temperature results calculated by the models were compared to the experimental results, and the assumption that the vapor penetration leads to deterioration of evaporator performance at high heat flux was proved.Vapor blanket thickness can be estimated through the saturation profile, which provides a simple and effective method.Capillary number is the key factor affecting the thickness of the vapor blanket; the greater the capillary number, the faster the vapor blanket thickness increases.
It is worth pointing out that although the method proposed in this paper can effectively estimate the vapor blanket thickness, the model still has some limitations.It cannot get more detailed information about the liquid-vapor interface, but only the approximate position of the interface through the distribution of saturation, and it cannot predict its shape.In the future, under the theoretical framework of the method proposed in this paper, we will consider the pore distribution in the wick and develop a two-dimensional model to study the liquid-vapor interface behavior more comprehensively.

sFigure 1 .
Figure 1.Principal scheme of a loop heat pipe.

Figure 1 .
Figure 1.Principal scheme of a loop heat pipe.

Figure 2 .
Figure 2. Schematic of the LHP experimental apparatus.

Figure 3 .
Figure 3. Variation of the evaporator heat transfer coefficient versus the applied heat flux.

Figure 3 .
Figure 3. Variation of the evaporator heat transfer coefficient versus the applied heat flux.

Figure 4 .
Figure 4. Schematic of the computational domain: (a) The capillary wick is completely saturated with liquid at low heat flux, and the evaporation process occurs at the wick-groove interface.(b)Vapor penetrates into the wick, which blocks liquid feeding and causes a partial dry-out spot in the wick at high heat flux.

Figure 4 .
Figure 4. Schematic of the computational domain: (a) The capillary wick is completely saturated with liquid at low heat flux, and the evaporation process occurs at the wick-groove interface.(b) Vapor penetrates into the wick, which blocks liquid feeding and causes a partial dry-out spot in the wick at high heat flux.

) gies 2023 ,Figure 5 .
Figure 5. Schematic of (a) cross-sectional view of the evaporator-CC structure and (b) heat flux tails in the wick.

Figure 5 .
Figure 5. Schematic of (a) cross-sectional view of the evaporator-CC structure and (b) heat flux details in the wick.

Figure 6 .
Figure 6.Diagram of the LHP working cycle.

Figure 6 .
Figure 6.Diagram of the LHP working cycle.

Figure 7 .
Figure 7. Schematic diagram of thermal resistance from wic

Figure 7 .
Figure 7. Schematic diagram of thermal resistance from wick fin to evaporator plate.

Figure 8 .
Figure 8. Extended meniscus for film evaporation in pore scale.(a) Extended film region close to the liquid-vapor interface.(b) A schematic diagram of the resistance for heat transfer by conduction across the liquid film and evaporation from the liquid-vapor interface.

Figure 8 .
Figure 8. Extended meniscus for film evaporation in pore scale.(a) Extended film region close to the liquid-vapor interface.(b) A schematic diagram of the resistance for heat transfer by conduction across the liquid film and evaporation from the liquid-vapor interface.

Figure 9 .
Figure 9. Heat transfer in a pore.

Energies 2023 , 19 Figure 10 .
Figure 10.The experimental results of vapor and evaporator variations with heat flux and comparison with results from two numerical models for different working fluids.

Figure 10 .
Figure 10.The experimental results of vapor and evaporator variations with heat flux and comparison with results from two numerical models for different working fluids.

Figure 11 .
Figure 11.The saturation distribution within the wick and the subregions of fluids.

Figure 11 .
Figure 11.The saturation distribution within the wick and the subregions of fluids.

Figure 12 .
Figure 12.The variation of the equivalent vapor blanket thickness with th

Figure 12 .
Figure 12.The variation of the equivalent vapor blanket thickness with the heat flux.

Table 1 .
The main parameters of the experimental LHP.