TRANSIENT MODELING OF CPL BASED ON MESOSCALED ANALYSIS BY LATTICE BOLTZMANN METHOD

A transient model of a Capillary Pumped Loop (CPL) is proposed to predict the thermal and dynamic behavior, which includes two coupled modules: the thermal model and the hydrodynamic module. The thermal module is based on the nodal method and the dynamic module of the liquidvapor interface is based on the molecular kinetic theory. Especially, a thermal non-equilibrium model is proposed to predict the temperature of evaporator, which is unlike some traditional numerical method associate with the CPL system. In the present model, vapor is treated as a compressible phase, so an equation of state (EOS) should be used to complete the controlling equations. A mesoscaled multiphase lattice Boltzmann model (LBM) is proposed to analyze the different working fluids, which can help us to choose appropriate EOS. Combining these two models for different scale, the behavior of ammonia is analyzed by the LBM. The processes such as the start-up, the temperature oscillations in the reservoir and the heat load change in the evaporator are investigated by the transient model of a whole CPL. The numerical results of the LBM show that the PengRobinson (P-R) equation is more appropriate for the ammonia and the results of the macroscopic transient model show that why the start-up of a CPL becomes a failure easily and temperature oscillations in the reservoir are disadvantageous to the operation of CPL. With this new mixed algorithm, the mesoscaled behavior of the working fluid and the operation mechanism of a whole CPL can be investigated in detail, which can provide guidance for the design of a CPL system.


INTRODUCTION
Capillary Pumped Loop (CPL) is a two-phase heat-transport device, which rely on the surface-tension forces induced by a fine pore wick to drive a working fluid (Muraoka et al., 2001;Cao and Faghri, 1994;Kaya et al., 2008). These devices use the same basic principle of the widely known heat pipes, i.e., closed evaporation-condensation cycle being maintained by capillary pumping. It is considered to be one of the most promising thermal-control equipment in the future (Faghri, 1995;Ku, 1993). Figure 1 shows a CPL functional schematic. A CPL consists of an evaporator, a reservoir (a compensation chamber), vapor and liquid tubes and a condenser. Only the evaporator contains wicks; the rest of the loop is made of smooth wall tubing.
CPLs were invented in USA (NASA) in late 1960s (Stenger, 1966). Nowadays although the CPL technologies have reached a high level of development especially in the experimental study (Chen and Lin, 2001;Liu et al., 2008), some issues about the fundamental mechanisms remain as pints of active research (Liao and Zhao, 2000). The main subjects of investigation are the device remains permanent pressure oscillations, start-up is easy to fail, and the abrupt change of the power in the evaporator can cause a bad reaction to the whole loop, etc. Thus it is necessary to simulate the transient process of the whole CPL with numerical methods, which can provide guidance for the experimental study of a CPL system.

Fig. 1 CPL function schematic
Transient process of a CPL has been mainly studied experimentally. There are only a few works about the numerical investigations. Muraoka and his co-authors developed a mathematical model based on the nodal method and analyzed the operational limits and dry-out failure mechanisms (Muraoka et al., 2001;). Pouzet and his co-authors proposed a global model to study the fundamental response mechanisms subjected to various heat load transients (Pouzet et al., 2004).
In this work, we present a complete numerical model of the whole CPL systems. Especially, we propose a mesoscaled lattice Boltzmann model to help us to choose some equations such as EOS of the working fluid, which is used in the transient numerical model of the CPL. Our present work aims at finding the appropriate EOS of the working fluid, studying the start-up, the capillary pressure change during the whole working process, investigating the influence of the power change and analyzing the effects of the temperature fluctuation in the controlled

Frontiers in Heat Pipes
Available at www.ThermalFluidsCentral.org 2 reservoir. The results can help us to investigate the operational characteristic of the whole CPL.

MESOSCALED NUMERICAL MODEL
In this paper, the mesoscaled multiphase lattice Boltzmann model is based on the SC model, which is a typical multiphase LBM (Shan and Chen, 1993). By the modification of the EOS of the SC model, a new multiphase lattice Boltzmann model which can simulate different EOS effect is developed.

Multiphase lattice Boltzmann model
LBM was derived from the lattice gas automata (LGA). It utilizes the particle distribution function to describe the collective behaviors of fluid molecules. It is proved that the lattice BGK method is consistent with the Navier-Stokes equation for fluid flow through Chapman-Enskog expansion (Benzi et al., 1992). As a flexible and effective numerical method, it finds increasing applications in simulating complex flow, heat, and mass transfer problems (Chen and Doolen, 1998;Xuan et al., 2007;Li et al., 2005;Orazio et al., 2003). An important advantage of the lattice Boltzmann method is that it can simulate the interfacial force between the two different phases conveniently and catch the interface of the two different phases clearly. So LBM has been widely used to simulate multiphase problems (Shan and Chen, 1993;Zhang and Chen, 2003;Swift et al., 1996, Yuan andSchaefer, 2006).
where x is the location vector, Δt is the time step, i=0,1,2…N is the direction of the discrete lattice speed, e i is the discrete velocity, τ is the dimensionless relaxation time, c=Δx/Δt is the lattice speed, Δx is the grid size and f i eq is the equilibrium distribution function defined as follows: where c s is the pseudo sound speed and for D2Q9 model (see Fig. 2) we can obtain that ω 0 =4/9， ω i =1/9,， i=1,2,3,4,，ω i =1/36，i=5,6,7,8 and c s 2 =(1/3)c 2 . By the Chapman-Enskog expansion, the density and the velocity can be calculated as For the SC model, the interfacial force can be calculated as where c 0 is a constant, which is associate with the model of the lattice, for D2Q9 model c 0 =6.0, g is a parameter which describes the strength of the interfacial force, ψ(ρ(x)) is the effective density in the SC model and decides the EOS of the model.
The interracial force is introduced into the LBM by the modification of the equilibrium velocity as follows: where U is the real velocity of the two phase flow.
By the Chapman-Enskog expansion, the EOS of the SC model is According to Eq. (6), it is can be found that for any EOS we can find a corresponding effective density in the LBM, which can be calculated as follows: In the simulation we must ensure the positive value under the square root in Eq. (7) (Zeng et al., 2010).

Relationship between lattice and real unit
Without loss of generality and for convenience, in most lattice Boltzmann model, the lattice unit is widely used. Here, for different working fluids, there are different physical parameters. For convenience, the lattice unit is used in the mesoscaled simulation in this paper. The transformation between the lattice and real unit is based on the state principle. We take the Van der Waals equation (vdW) as an example to illustrate the transformation.
The format of vdW in real unit is as follows: where the superscript "real" represents the real unit, a real , b real are two equation parameters associated with the critical condition.
According to the state principle, the relations between two different units are as follows: where subscript "c" represents the critical condition and the superscript "Lu" represents the lattice unit. Substitute Eq. (9) to Eq. ( So according to Eq. (14), we prove that the format of the vdW is same under different unit system. Use this principle above, the parameters and the format for different equations of state under lattice unit can be obtained conveniently.
vdW EOS is a very simple EOS of the real gas, which can describe the behavior of CO 2 very well, but not good at other gases. It is not appropriate to use a unified EOS for simulating different working fluids in the numerical model of a CPL system. We must choose a suitable EOS for our working fluids. The lattice Boltzmann model above can help us analyze different EOS very conveniently.
The other EOSs analyzed in this paper under lattice unit are given below (for simplicity the superscript "Lu" is omitted): Redlich-Kwong (R-K ) EOS: In the above EOS, the R-K EOS, RKS EOS and P-R EOS are classified as cubic EOS, in which the R-K EOS is a two-parameter EOS, using two parameters (T c , p c ); the RKS EOS and P-R EOS are threeparameter EOS, with an additional parameter ω, which is the acentric factor. This additional parameter gives us more flexibility in modeling different fluids.
The parameters for lattice units in LBM are fixed, which is that Δx=Δt=1, c=Δx/Δt=1. For all of these cubic EOS, we set a=2/49, b=2/21 and R=1 in our simulations. The values of these parameters in real units can be calculated according to the processes mentioned in section 2.2.

Simulation results of LBM
The coexistence curves for the above three EOS are calculated by this multiphase LBM. Comparison between the simulated results and the experimental data can help us to know which EOS is closer to the experimental data. Ammonia is simulated in this section. For RKS EOS and P-R EOS, the acentric factor for ammonia is ω ammonia =0.25.
(a) (b) Fig. 3 Coexistence curves of ammonia for three EOS (a) vapor and liquid branch (b) vapor branch Figure 3 gives the comparison of the saturated density obtained from simulations for the P-R EOS, the R-K EOS and the RKS EOS with the experimental data from the steam table. From this comparison with experimental data, we can see that for the liquid density (see Fig.  3(a)) the P-R EOS gives better results and the RKS EOS gives the worst results, while for the vapor density, the P-R EOS also gives best results, but compared with the liquid branch, the results of the three EOS for vapor phase are much better. According to Fig. 3(a), we can find that for different EOS, the temperature range which can be simulated by LBM is different. In our simulation, the temperature range of the R-K EOS is the smallest. This is associated with the model of the LBM. The maximum relative error of the liquid density is 7.71% for the P-R EOS. However, we can generally consider the P-R EOS to be superior to the other two EOS for the working fluid of ammonia. Furthermore, unlike the P-R EOS, the R-K EOS is a two-parameter EOS. Therefore, no matter what fluid is being simulated, the R-K EOS gives the same coexistence curve for the reduced properties. The P-R EOS is a three-parameter EOS, and the third parameter (the acentric factor) gives us more flexibility in simulating different fluids. For example, by setting ω=0.344, we are essentially using the properties of water. Generally, the multiphase lattice Boltzmann model mentioned above can help us to choose an appropriate EOS for the working fluid which is used in the transient model of a CPL system.

NUMERICAL MODEL OF A CPL SYSTEM
Based on the work by Qian et al. (Qian et al., 2006), we developed a mathematical model for describing the transport of heat and mass inside the loop during all its operational regimes. This model comprises two different modules: the thermal and dynamical modules. The two modules are developed by the following assumptions: (1) the evaporation in the wick is uniform along the axis and perimeter; (2) the temperature gradient exists only along the radial; (3) the temperature in the reservoir is equal to the set point temperature; (4) the vapor is compressible and the liquid is uncompressible; (5) the gravity is ignored.

Thermal model
The thermal module, which computes the distribution of temperatures as a function of time, is based on the nodal method. Each node is characterized by a thermal capacitance, and is connected to other nodes by means of generalized thermal conductance. The nodes distribution is shown in Fig. 4.

) Evaporator
A thermal non-equilibrium model is proposed to computer the temperature of evaporator. The solid wick and the working fluid in the module is divided into different isothermal nodes accordingly, the heat transfer between the solid wick and the fluid can be investigated.
The nodes distribution in the evaporator is shown in Fig. 15. Each equation represents the thermal balance at the nodal level, and is given by Cover plate: where Q e is the heat load.
Porous wick: where M i , T i , H i are mass, temperature and enthalpy of node i, respectively, k ij is the conductivity between node i and j, i m  is the mass flow rate and the volumetric heat transfer coefficient h v is given by (Jiang and Ren, 2001) where ε is the porosity of the wick, d p is the particle diameter of the wick, k f is the conductivity of the working fluid.
(2) Vapor tube, liquid tube and condenser The energy equation in the tube wall is as follows: where h i is the convection coefficient between wall and the working fluid.

Hydrodynamic model
The hydrodynamic module determines the pressure and phase distributions inside the loop, as well as the flow rates of liquid and vapor. The thermal and hydrodynamic systems are coupled since, on the one hand, the transport of heat depends on the fluid flow, on the other hand, the pressure distribution, that drives the fluid along the loop, is a function of the temperature distribution.
(1) Evaporator According to the molecular kinetic theory, the evaporation rate of liquid in the evaporator is given by (Zhang and Wang, 2002)   where A eva is the vapor-liquid interface area, R g is the gas constant, T i is the interface temperature, p v is the vapor pressure above the interface, p s (T i ) is the saturated pressure under T i , r is the mean meniscus curvature radius and p d is the disjoining pressure, which is ignored in this model.
The pressure lose in the cylindrical wick is given by where l m  is the mass flow rate in the wick, r o , r i are the inner and outer radius of the cylindrical porous wick, respectively, K p is the permeability of the porous wick, M w is the mass of the working fluid in the wick and L w is the length of the porous wick.
The pressure drop in the vapor groove is given by where n is the number of the grooves and r g is the hydraulics diameter of the groove.
(2) Condenser Because of microgravity and low Bond number in the condenser, the two phase flow in the condenser is assumed to be the annular flow, which is confirmed by the Allen et al (Allen and Hallinan, 2001). The diameter of the tube we used is very small. So that the interface tension is very high and it leads to a very thin liquid film. In our model, some assumptions for the condenser should be proposed, which is as follows: (1) The thickness of the liquid film does not change with time; (2) The temperature of the liquid film equals to the temperature of the tube wall; (3) The superheat of the steam is low, so the two-phase zone happens just at the entrance of condenser; The mass conservation equation is given by out in where L con , L con,v are the length of the condenser and the vapor region respectively.
The momentum conservation equation is as follow where f v , f l are the friction factor of vapor and liquid respectively, whose form are as follow: Similar with the evaporation rate in the wick, the condensation rate in the node i is given by The length of the vapor region in the condenser can be calculated as follow (

3) Vapor and liquid tube
The vapor in the vapor tube is also considered to be compressible. The EOS and the mass and the momentum conservation equations are similar with those in the condenser. While for the liquid tube, because there is no phase change and two-phase flow in it, the equations in the liquid tube are much simpler. Thus we do not give them in detail here.
(4) Reservoir In this paper, the fluid temperature in the reservoir is at the set point temperature T r , and the mass conservation equation is as follows: where L r , L r,v are the height of the reservoir and the height of vapor in the reservoir respectively.
Other equations are similar with Eq. (34)-(37). Fig. 6 Schematic of the numerical procedure As mentioned above, the complete numerical algorithm contains two different modules: (1) mesoscaled model based on LBM (2) macroscopic transient model for the whole CPL system. The mesoscaled model in this paper is used to help us find a better equation of state for the working fluid, which is used in the latter system level simulation. In this module, the only input data is which working fluid is used in our simulation. After the mesoscaled simulation, the suitable EOS has been chosen for the transient system module. The system module includes two coupled models: the thermal and the dynamic model. The change of the temperature may cause the change of the mass flow rate and the pressure, while the change of the pressure and the mass flow rate also affect the temperature. It is a complicated process to solve the energy and momentum equations. The detailed numerical procedure is shown in Fig. 6.

RESULTS AND DISCUSSION
The geometry parameters of the CPL system for the simulation are based on our laboratory. As we can see from Fig. 1, the whole system includes a cylindrical evaporator, a condenser, a temperature control reservoir, a liquid loop and a vapor loop. The wall of the evaporator is made of stainless steel with 12 longitudinal grooves. The mean hydraulics diameter of groove is about 1.2 mm. The sintered porous wick is made of nickel powder, and the particle size is about 5~10 μm. The permeability and porosity of the porous wick are about 5×10 -13 m 2 and 0.5, respectively. The length of evaporator is about 100 mm and the outer diameter is 12 mm. The condenser tube is combined with a radiated plat. The length of this tube is 1.2 m, whose diameter is 2.6 mm. The reservoir is controlled by a temperature controlled equipment, whose accuracy is 0.5 K. the volume of the reservoir is 5.25×10 -5 m 3 . The liquid and vapor tube is made of stainless steel (No. 304), whose length is 3 m and the outer diameter is 2.6 mm. The wall thickness of the tube is 0.7 mm. The working fluid is ammonia. All the components of the device are arranged on the same horizontal level to minimize the effectshe whole system is divided into 237 nodes: evaporator: 37; vapor tube: 40; condenser: 40; liquid tube: 40 and the wall of tube: 100. of gravity. Because of the evaporation in the evaporator, the time step should be enough small, in this paper, we choose 1 ms.

Fig. 7
Evaporator temperature versus time Figure 7 gives the evaporator temperature versus time during the start-up of the CPL. The numerical results show that at the begin of the start-up, the temperature of the evaporator rises rapidly, then it slows down and keeps around a value with small fluctuation. It can be explained that at the begin of the evaporation, the condensation rate in the vapor tube is very small, so that, the liquid-vapor interface moves to the condenser quickly, thus the liquid rate in the evaporator is very high, which indicates that the capillary pressure is very high. At this time, there are much liquid goes into the evaporator from the reservoir. The vapor pressure in the reservoir is relatively higher than the saturated pressure, so that, the liquid from the reservoir causes the risen of vapor pressure and the evaporation temperature. The small fluctuation may be due to the fluctuation of the liquid-vapor interface in the evaporator. Comparison between the temperature on the heated surface and on the interface of the wick and the vapor groove shows that the temperature on the heated surface is much higher than the temperature inside the wick and the fluctuation is much smaller. It can be explained that during the evaporation, the liquid-vapor interface in the wick is unstable, which can leads to more intensive fluctuation of the inner temperature near the working fluid. Figure 8 shows the pressure difference between the two sides of the evaporator and the liquid-vapor interface position in the vapor tube during the start-up of the CPL. It show that there is an abrupt rise of the pressure difference at the beginning of the start-up, and it keeps rising to a maximum, then it goes down through a intensive fluctuation to a stable value. When the pressure difference reaches a stable value, there is some small fluctuation. The abrupt rise of the pressure difference may destroy the liquid-vapor interface and cause the start-up failure. It can be explained that at the beginning of the start-up, the evaporation in the evaporator and the condensation in the vapor tube happen simultaneously. The liquid mass flow rate is very high, which leads to a very large pressure difference. When the liquid-vapor interface in the vapor tube goes into the condenser, there is some intensive fluctuation of the pressure difference. This may be attributed to the fluctuation of the condensation rate (see Fig. 10). Fig. 10 gives the mass flow rate in the different position of the CPL system. It shows that when the CPL starts up, the evaporation rate and the mass flow rate of the working fluid in the liquid tube rise very quickly, but because the liquid-vapor interface is in the vapor tube at the beginning, the condensation rate is very small. When the liquidvapor interface goes into the condenser, the condensation rate rises quickly, and catches up with the evaporation rate generally. Finally, these three mass flow rates reach a steady-state. We find that the position of the liquid-vapor interface has great influence in the mass flow rate, such as evaporation rate, condensation rate and so on. Also it may cause the fluctuation of these flow rates. The temperature in the liquid tube and the vapor volume of the reservoir are shown in Fig. 9. According to the curve of the vapor volume in the reservoir, we can find that at the beginning of the start-up, the vapor volume decreases very quickly, the slope of the curve is very high, which indicates that the change rate of the vapor volume is very quick. It is conclude that at the beginning of the start-up, there is much liquid goes into the reservoir due to the large evaporation rate at the beginning. Fig. 9 also shows that at the beginning of the start-up, the temperature in the liquid tube decreases very slowly, but when t>60s, the temperature will decrease very quickly.  Figure 11 shows the influence on the CPL operation characteristic of the heat load. It shows that the pressure difference between the two sides and the temperature of the evaporator increase with the heat load. A larger heat load leads to a larger pressure difference, which is not good to the operation of the CPL system. The larger pressure difference is easy to destroy the liquid-vapor interface in the wick and then the dry-out may happen. Comparison between the heated surface temperature and the wick temperature, we can find that the heated surface temperature is more sensitive to the heat load. The evaporation rate, the condensation rate, the liquid-vapor interface position and the vapor volume in the reservoir under different heat load are shown in Fig. 12 and 13, respectively. According to Fig.   12, the evaporation and condensation rate, as well as the fluctuation of the flow rate increases with the heat load. Especially, the time of the abrupt rise of the condensation rate is different under the two heat loads. The time when the abrupt rise happens under Q e =40 W is shorter than the time under Q e =20 W. It can be explained that the time when the liquid-vapor interface in the vapor tube under larger heat load reaches the condenser is much shorter than the time under smaller heat load. It is can be seen in Fig. 13. Because of the large evaporation rate under large heat load, the flow rate in the tube is very large and there is much working fluid inrush into the reservoir, which leads to a very quick decrease of the vapor volume in the reservoir. This is also shown in Fig.  13.
(a) Temperature and mass flow rate (b) Liquid-vapor interface and pressure difference  Figure 14 presents the transient CPL operation for the pulsed heat input. At time is 30 s the heat load is suddenly increased from 10 W to 40 W, and then at time is 60 s, the heat load is suddenly decreased from 40 W to 30 W. When the heat load is suddenly increased, the temperature in the evaporator increases to a higher level in few seconds. The heated the heated surface temperature is more sensitive to the pulsed heat input. The evaporation and the condensation rate are under an equilibrium state before the change of the heat load. When the heat load increases suddenly, they are both increase to a higher level, and than reach a new equilibrium state. As can be seen in Fig. 14(b), the pressure difference between the two sides of the evaporator keeps around a small value before the suddenly increase of the heat load. When the heat load change to 40 W, the pressure difference increases to 1500 Pa suddenly, and then decreases to a new steady-sate through some fluctuation. When the heat load decreases to 30 W, there is also an abrupt decrease of the pressure difference. All these changes of the temperature and the pressure are both not good to the CPL operation.  Fig. 15 The influence of the pressure fluctuation in the reservoir During the operation of the CPL system, the temperature of the reservoir is fixed at the set point (313K in our model). But in fact, it has some fluctuation. A small change of the temperature in the reservoir may cause a large change of the saturated pressure, which has a great influence on the operation of the CPL system. Fig. 15 presents the temperature, the pressure difference and the liquid-vapor interface versus the fluctuation of the saturated pressure in the reservoir. The frequencies of the fluctuation in the reservoir are 0.2 Hz and 0.1 Hz respectively, and the heat load Q e =20 W. According to Fig. 15(a), the temperature of the evaporator increases with the pressure increasing. But the pressure and the liquid-vapor interface are opposite. When the pressure of the reservoir begins to increase, the pressure difference decreases first and then increases and the liquid-vapor interface falls back to the vapor tube first. It may be explained that when the pressure in the reservoir increases, the pressure at the entrance of the evaporator increases first, if the pressure exceeds the saturated pressure under the operation temperature, the evaporator will stop working, and the CPL system will reach an operation failure. The fluctuation range of the pressure difference and the liquid-vapor interface position increases with the frequency of the fluctuation in the reservoir, but the temperature is not the same as the pressure difference and the liquid-vapor interface position. So that it indicates that the high frequency of the fluctuation in the reservoir has a great influence on the whole system performance, which should be avoided in the experiment.

CONCLUSIONS
In this paper, a mesoscaled multiphase lattice Boltzmann model and a transient numerical model of a CPL is developed to investigate the thermodynamic behavior of the working fluid and the thermal as well as dynamic behavior of the whole CPL system. With the lattice Boltzmann model, the coexistence curve of ammonia is analyzed, which shows that the P-R EOS is more suitable to describe the thermodynamic behavior of ammonia. Then the P-R EOS is introduced to the transient model of the CPL. With this transient model, the startup, the temperature oscillations in the reservoir and the heat load change in the evaporator are investigated. The results show that the pressure drop is highest at the beginning of the start-up, which may be one of the main reasons cause the start-up failure. The heated surface temperature is more sensitive to the heat load. The very high heat load, the pulsed heat input and the oscillations of the reservoir temperature adverse to a CPL system, which should be avoided in this two phase loop. All the results above can give us some guidance for the experiment design of a CPL system.

ACKNOWLEDGEMENTS
This work was supported by the National Natural Science Foundation of China (Grant No. 50176018) and Program for New Century Excellent Talents in University. The authors would like to thank the helpful discussions from Dr. Jiyu Qian about the system model. A special thank is also given to Dr. Jianbang Zeng for his patient and useful discussions on multiphase lattice Boltzmann model.