LOW-ORDER MODEL OF THE DYNAMICS AND START-UP OF A PULSATING HEAT PIPE

A simple open-loop pulsating heat pipe model is proposed, which allows to analytically determine the start-up behavior by a linear stability analysis. Two distinct types of instability can occur in such a pulsating heat pipe: oscillatory and non-oscillatory. This paper demonstrates that for bubbles consisting of non-condensible gas, large temperature gradients along the wall are required to achieve start-up, whereas start-up is fairly easy to achieve when there is only a single working medium that forms bubbles from its vapor. The study also finds that surface tension as such only influences start-up indirectly, while contact angle hysteresis dampens out any instabilities.


INTRODUCTION
With recent progress in miniaturization of integrated circuits come increased cooling requirements. This has raised the interest in pulsating heat pipes (PHPs). PHPs have also been investigated for use for battery cooling in electric vehicles (Singh et al., 2021). PHPs have the potential to remove high heat fluxes from small surfaces at low temperature gradients. In addition, PHPs are easy to manufacture, and they do not require any moving parts, so they are not prone to mechanical failure. Unfortunately the physics of PHPs is still poorly understood, preventing reliable prediction of performance and thus commercial use.
A pulsating heat pipe is a tube with multiple bends (Fig. 1). Alternatingly, the bends lie in evaporator or condenser zones. The tube is filled with a working fluid, which in general is in two-phase state, i.e. * Corresponding author. Email: schily@tfd.mw.tum.de liquid and gas coexist. The diameter of the tube is sufficiently small that the two phases cannot pass one another, but form distinct bubbles and slugs. When heat is added in the evaporator zones, bubbles may expand due to increased temperature, but mostly due to partial evaporation of the adjacent liquid slugs. This moves the slugs and other bubbles towards the condenser zones, where heat is removed, and bubbles contract. When continuously heating the evaporator and cooling the condenser, one would naïvely expect that the motion of bubbles and slugs goes on until there is no more liquid in the evaporator sections and no more vapor in the condenser sections, and then dies down. Instead, it is often observed that the PHP enters a state of continuous, self-sustained, and irregular oscillatory motion.
In the literature on PHPs, two foci have developed: On the one hand, there are many papers that investigate the performance of PHPs experimentally (Quan and Jia, 2009;Ma et al., 2008;Song and Xu, 2009;Qu et al., 2009;Khandekar et al., 2003;Singh et al., 2021), often with the goal to correlate performance with design and operational parameters, such as tube diameter, filling ratio, and orientation with respect to the field of gravity. Also note recent experimental investigations on PHP operation in micro-gravity by Taft and Irick (2019). More recently, performance improvements were achieved by deviating from the classical structure of a PHP and introducing e.g. multiple branches or cross-talk between sections (Fairley et al., 2015;Shi et al., 2011;Stevens et al., 2019). Also related are wire-bonded micro heat pipes (Sobhan and Peterson, 2019).
On the other hand, the dynamics (i.e. the movement of slugs and bubbles) of PHPs were often modeled by a set of ordinary differential equations treating each slug and bubble as a separate, homogeneous entity (so-called lumped model). Subsequently these models were numeri-cally simulated in the time domain (Shafii et al., 2001;Kim et al., 2005;Dilawar and Pattamatta, 2013). The authors are not aware of a model that yielded satisfactory quantitative results for PHP dynamics in comparison to experiments. Bridging this gap between numerical and experimental results proved to be a challenge. Partially, this is due to the now well-established fact that PHPs display chaotic behavior (Maezawa et al., 1996;Dobson, 2004;Song and Xu, 2009;Qu et al., 2009). But also the multitude of physical effects interacting with one another in PHPs, and the dynamics of many bubbles and slugs moving through several evaporators, condensers, and adiabatic sections impede the development of a sufficient understanding of the working principle of a PHP.
More recently there have been promising efforts to investigate PHP dynamics with CFD (Pouryoussefi and Zhang, 2016;Bhagat and Deshmukh, 2021).
To avoid loosing overview of the problem, this study is confined to the start-up mechanism of a very reduced example of a PHP, with only one liquid slug and two bubbles -i.e. there are no large amplitudes, no non-linearities, no chaotic behavior. The start-up mechanism is a key part to understanding the dynamics of a PHP, since it allows insight into the origin of the kinetic and potential energy stored in the continued oscillation of a PHP. The spontaneous onset of motion indicates that a linear instability might be involved 1 .
Studies on similarly minimalist setups have been performed, e.g. by Dilawar and Pattamatta (2013), and also by Dobson (2004), who investigated an open-ended single-bubble PHP experimentally and numerically, remarkably observing chaos in both the experiment and the model. Also the groundbreaking categorization of flow regimes in a PHP by Khandekar and Groll (2004) was performed using an extremely reduced setup.
On the start-up behavior, there have been experimental investigations in the past e.g. by Quan and Jia (2009), and also some numerical ones like Yang and Luan (2012). This work adds a theoretical study of the physics behind PHP start-up.
The present paper will first investigate the start-up behavior of a simplistic PHP with bubbles consisting of a non-condensible gas by analyzing its linear stability in the equilibrium state (Section 2). Then successively more physics and complexity will be added to the model (vapor bubbles with phase change in Section 3, surface tension effects in Section 4), while repeating the linear stability analysis to assess the effects of the added phenomena on start-up behavior.
This study finds that there are two types of instability that can in principle cause the onset of motion in PHP model under investigationan oscillatory instability, and a non-oscillatory one. Depending on the temperature gradient along the wall, one of them may exist, but never both. Another result is that it takes a prohibitively strong wall temperature gradient to support any of said instabilities, if the bubbles consist of non-condensible gas. When allowing phase change, instabilities become much more likely. Surface tension does not directly influence the conditions for start-up by linear instability in case the wall consists of a high-energy material, i.e. if there is no contact angle hysteresis. If there is contact angle hysteresis, no linear instabilities occur at all.

MINIMUM WORKING EXAMPLE OF A PULSATING HEAT PIPE
This section serves to introduce a simplistic PHP setup (a so-called "minimum working example"), and describe it mathematically with a set of simple ordinary differential equations. To describe the start-up behavior, 1 Note that investigating linear stability disregards non-linear triggering (Juniper and Sujith, 2018), i.e. the onset of motion by non-linear instability. In the absence of linear instabilities, non-linear instabilities require perturbations of finite magnitude. It is safe to assume that a sufficiently large external perturbation, e.g. shaking of the entire PHP, is not always present (as opposed to infinitesimal perturbations required by linear instabilities). Therefore this mechanism cannot explain the spontaneous start-up of PHPs as reported in the literature. However, a sufficient perturbation for non-linear triggering might result from the formation of a bubble within a liquid slug. The investigation of bubble formation and the relevance of non-linear triggering in PHPs is, however, beyond the scope of this paper. the model is linearized around the equilibrium point to obtain a characteristic equation. The zeros of the characteristic equation correspond to eigenvalues. An instability, i.e. onset of motion, occurs if there is at least one eigenvalue with a positive real part.
The geometry of the minimum working example is shown in Fig. 2. The object of the study is a U-shaped, open-loop PHP with closed ends. It has two bubbles of non-condensible ideal gas, and one incompressible liquid slug between the bubbles. There is an arbitrary wall temperature distribution. The U-bend does not occur explicitly in the model. This allows to consider a straightened model with a proper axial coordinate (Fig. 2).
Such a setup can show oscillatory behavior if the bubbles act as heat engines, i.e. heat is added after compression and rejected after expansion. This mechanism was first described by Rayleigh (1878) as the founding statement of thermo-acoustics. The work output generated by the oscillating bubble is then used to further increase the oscillation amplitude or to overcome viscous friction. The latter is, however, neglected in the following. The analysis in the present section serves the purpose of clarifying the circumstances required to add and reject heat with the right timing, so the Rayleigh criterion is satisfied.
The nomenclature is presented in Fig. 3. The indices 1 and 2 refer to left and right bubble, respectively. Each bubble is characterized by respective uniform values of pressure p, volume V , temperature T , extensive internal energy U and mass m. For the case of non-condensible bubbles, the mass in both bubbles is the same mg and constant. This constraint will be dropped in the following sections.Q denotes heat flow rate (again, extensive not mass or area specific). For the present section, only heat transfer from the wall into the bubbles is considered.
x is the axial displacement of the center of mass of the slug from its neutral position, w the velocity of the slug, and m l its mass. In this section the mass of the slug remains constant, but also this constraint will be relaxed in the following sections.
With x, not only the slug position, but also the bubble volumes can be expressed by the equilibrium length L of each bubble and the cross sectional area A of the PHP: where the index j refers to the left (j = 1) and right bubble (j = 2), respectively. The ± expresses that for the left bubble (j = 1), x is added, whereas it is subtracted in the case of the right bubble (j = 2).
The following governing equations describe the setup: where t is time. In this simplistic model, friction, capillary effects, and heat transfer in the axial direction are neglected. Heat transfer between bubbles and wall only depends on temperature difference. Changes in the surface area available for the heat transfer and changes in the heat transfer coefficient are disregarded, such that the constant K can completely describe wall-to-bubble heat transfer: K will later be varied as a part of a non-dimensional parameter. Tj is the bubble temperature, andTw,j is the average temperature of the wall section in contact with the respective bubble:  Fig. 3 Nomenclature: Indices 1 and 2 for left and right bubble, respectively, slug mass m l , displacement x and slug velocity w. Mean wall temperatures Tw1 andTw2 of the bubbles are obtained by integration and thus also the mean temperature depends on the integration limits.
As presented in the lower part of Fig. 3, not only the normalization factors, but also the integration domains ∂Vi depend on the slug displacement x. Therefore the mean wall temperatureTw,j also depends on the slug position (and thus on time t), although the wall temperature profile itself is constant Tw = Tw(t).
Heat transfer between the liquid slug and the wall may be neglected, since the slug is assumed incompressible and cannot evaporate. Consequently the slug temperature cannot influence the dynamics of the PHP. The assumptions that there is no phase change and that there are no capillary effects will be dropped in Sections 3 and 4, respectively.
After inserting the definition of heat flow rates -Eqs.
(2) are linearized around the equilibrium state. The equilibrium state must be chosen as the initial state for the linearization, because otherwise the setup would not describe the start-up of a PHP, but some transient state. Since geometry and boundary conditions are symmetric, the equilibrium should be, too: x = 0, peq ≡ p1 = p2. Since mass and volume of both bubbles are the same, their temperatures must also be equal (ideal gas law): T1 = T2. Equilibrium means that there is no change of state, i.e. in particular the internal energy of the bubbles may not change. Therefore there may not be any heat transfer in the equilibrium state, requiring the mean wall temperature in each bubble to be the same as the bubble temperature. Since all temperatures occurring are the same, the equilibrium temperature is defined as Teq ≡ T1 = T2 =Tw1 =Tw2. As the system regarded is a U-shaped tube, the wall temperature distribution in the straightened model must be symmetric with respect to the center of the tube (Fig. 2). Therefore the wall temperature gradients of the two bubbles in equilibrium relate to one another by Ultimately, the linearization results in the following set of equations: where the denotes small, first-order perturbations of the respective variable. All other variables are not a function of time anymore. Now one closes the system with the caloric equation of state U = mgcvT = mgRT / (γ − 1), the ideal gas law pV = pA (L ± x ) = mgRT , and the kinematic constraint V 1 = −V 2 = Ax : Non-dimensionalizing with x + ≡ x /L, x * ≡ x/L, t + ≡ t/τosc, w + ≡ w τosc/L, U + ≡ U / (mgcvTeq) = U /Ueq, and T * ≡ T /Teq allows not only a more general statement, but also to reduce the number of variables: τosc ≡ m l L 2 /(mgRTeq) defines the characteristic time scale of the oscillator in order to avoid non-dimensional parameters in the momentum balance. In physical terms this corresponds to the natural period of the (adiabatic) gas-spring oscillator of the same dimensions times a constant factor: τosc = τgas spring √ 2γ/(2π). In addition, the characteristic time scale (the relaxation time) of the heat transfer is defined as τq ≡ mgcv/K. The ratio of the two time scales yields the nondimensional time lag τ + ≡ τq/τosc. These measures have reduced the complexity of the system to three non-dimensional constants, one of which is the gas property γ (ratio of specific heats).
The general solution of this system of linear ordinary differential equations is the superposition of all non-dimensional eigenvectors z + : Of particular interest are unstable solutions, as the existence of an unstable solution indicates that the solution will grow to finite values from an infinitesimal perturbation, which would physically be observable as the start-up of the PHP. For an unstable system, there must be at least one eigenvalue λ + with a positive real part λ + = σ + > 0, where real and imaginary part of the eigenvalue λ + = σ + + iω + correspond to non-dimensional growth rate σ + and frequency ω + of the solution, respectively (both real-valued). To find the eigenvalues one rewrites the system of equations as: with the coefficient matrix M . Then one solves the characteristic equation: with the unity matrix I. It is immediately clear, that there is always a stable real-valued eigenvalue at λ + = −1/τ + < 0, corresponding to a stable, non-oscillatory eigenmode. Looking for instability, one can focus on the other factor, a polynomial of third order. For now, this study confines itself to merely finding the stability limit, i.e. the minimum requirements for PHP start-up. Hence, one looks for combinations of γ, τ + , and ∂T * w /∂x * that allow σ + ≡ λ + = 0 and therefore λ + = iω + . Insertion yields: This equation can be split up into real and imaginary part and yields two conditions that a marginally stable eigenvalue must fulfill: Strictly speaking, it is not yet known, whether these two conditions are actually stability limits. Instead, one or both of them could be a minimum or maximum of the polynomial in λ + given as the second factor of Eq. (11).
To check whether the two conditions found really are stability limits, one can investigate the sensitivity, i.e. the derivative of the nondimensional growth rate with respect to the non-dimensional wall temperature gradient dσ + /d ∂T * w /∂x * at the temperature gradient where σ + = 0. Only for zero derivative, both larger and smaller wall temperature gradient have the same stability, i.e. there is no stability limit. If the derivative is positive, increasing the wall temperature gradient from its original value at σ + = 0 results in a transition from stable oscillation into instability. For negative derivative, decreasing the wall temperature gradient corresponds to transition into instability. The analysis in Appendix A shows that: Thus, for an oscillatory eigenvalue, the criterion ∂T * w /∂x * = 1 − γ is always a stability limit.
Remarkably, with ω + = ± √ 2γ, the ratio of the characteristic time scales τ + does neither influence stability itself, nor the non-dimensional frequency of oscillation at marginal stability.
From a physics point of view, this result can be interpreted as follows: If the wall temperature gradient is negative and sufficiently strong 2 , heat is added to the bubble during and after compression 3 . The reason is that the compressed bubble is only in contact with hot zones of the wall. The gas temperature of the bubble increases due to compression, but whenever Eq. (16) is satisfied, the increase of the mean wall temperature due to displacement to hotter section of the wall is stronger than the increase of the gas temperature due to compression. Eventually, the pressure increase in the bubble due to compression and heat addition builds up a restoring force that stops the slug and then pushes it back. Due to its kinetic energy, the slug overshoots the neutral position and the bubble expands, such that the corresponding opposite effect happens: Although the gas temperature decreases, the decrease of the mean wall temperature dominates and results in cooling after expansion.
To illustrate the effect, Fig. 4 shows the non-dimensional perturbations of volume V + 1 = x + , pressure p + 1 = U + 1 − x + and heat fluẋ Q + 1 = ∂T * /∂x * x + − U + 1 /τ + for the left bubble (index j = 1). Since internal energy increase and decrease from heat transfer has a delay (relaxation time), heat is added when the pressure is high and rejected at low pressure. Similar to a heat engine, this thermodynamic cycle generates work, which is then stored in the kinetic energy of the slug. In thermo-acoustics, this is expressed by the Rayleigh criterion (Rayleigh, 1878), which is also shown at the bottom of Fig. 4. The integral of p + 1Q + 1 over one cycle is positive, therefore the process generates work, i.e. increases its amplitude.
For the non-oscillatory eigenvalue, the sensitivity of the growth rate with respect to the temperature gradient can similarly be determined: i.e. if the wall temperature gradient is positive and sufficiently large 4 , one bubble grows more and more at the cost of the other 5 . To understand the physics behind this effect, consider the quasi-steady caseTw ≈ T , i.e. heat transfer is sufficiently fast that the bubble temperature is always equal to the mean wall temperature. Upon compression and the corresponding displacement, the bubble is exposed to a section of the wall with lower temperature such that its temperature decreases. If the wall temperature gradient is sufficiently strong a pressure decrease may result despite the volume reduction of the same bubble. The other bubble suffers the inverse effect: pressure increases despite volume expansion. The pressure difference moves the slug such that compression and expansion of the bubbles are accelerated once they started. Both instabilities, oscillatory and non-oscillatory, are in principle in favor of the start-up of the PHP, since they create motion.
The above results can be put into context by numerically solving the characteristic equation and varying the three parameters. Figure 5 displays what happens to the eigenvalues when varying ∂T * /∂x * at γ = 4/3 and τ + = 1 (colored points). The constant eigenvalue λ + = −1/τ + is represented as a black ×. The red entries in the diagram denote eigenvalues computed from a configuration with ∂T * w /∂x * ≤ 1 − γ (recall the type: oscillatory unstable eigenmode). There are always two stable eigenmodes with real-valued eigenvalue, and a pair of unstable, oscilla- tory eigenmodes, (i.e. complex valued eigenvalues). When increasing the temperature gradient to lie within the interval 1 − γ ≤ ∂T * w /∂x * ≤ 1, there are only stable eigenmodes (black entries). Finally, for ∂T * w /∂x * ≥ 1 (blue) there exist the known stable eigenmode with the real eigenvalue λ + = −1/τ + (black ×), a pair of stable oscillatory eigenmodes, and one unstable eigenmode with real eigenvalue.
In gray, Fig. 5 presents corresponding results for four additional values of τ + : 2, 1/2, 1/3, and 1/4. The pattern evolving confirms that τ + does neither influence whether the system is stable or not, nor does it affect the oscillation frequency at the stability limit for oscillatory instability -the curves for constant τ + intersect on the imaginary axis. Physically this means that stability of the simplistic model itself does not depend on all those quantities that only occur in τ + , e.g. filling ratio and heat transfer coefficients. This stands in contradiction to experimental findings, which will be discussed in Section 5.3.
τ + however has an effect on growth rate and frequency, if the system configuration does not perfectly match the stability limit. Figure 6b shows that for smaller τ + , the eigenvalue pattern displays two exceptional points (black circles). Exceptional points are points where two or more eigenvalues and their corresponding eigenvectors coincide (Heiss, 2004(Heiss, , 2012. Whereas it is apparent from Fig. 6b that the eigenvalues coincide, it is mathematically tedious to check whether also the eigenvectors coincide. Therefore in Appendix B, the method is sketched, but not all equations are explicitly shown. When increasing the wall temperature gradient the system initially has a pair of unstable, oscillatory eigenmodes and two non-oscillatory, stable ones (the evolution of one real-valued eigenvalue is displayed in red, the other is represented as black × and always at λ + = −1/τ + ). Then the pair of oscillatory eigenvalues becomes stable (black), and subsequently encounters the first exceptional point turning it into two real eigenvalues. One of them then becomes more and more unstable and moves to the right along the real axis until it reaches stability limit and becomes unstable again (blue). The other becomes more stable and unites with the moving one of the two originally stable real-valued eigenvalues in the second exceptional point, to form again a pair of stable oscillatory eigenvalues (blue).
Following the evolution of eigenvalues for larger τ + (Fig. 6a) is less complicated. There are no exceptional points, therefore the initially unstable pair of oscillatory eigenmodes grows stable, whereas one of the two stable eigenmodes (real eigenvalues) grows unstable without interacting with the oscillatory eigenvalue pair. All patterns have in common that the real eigenvalue turns unstable after the oscillatory eigenvalue pair has become stable, leaving a range of wall temperature gradients, where the system is stable. Figure 7 summarizes the behavior of the PHP model in a stability map. The stability limits describe above are marked as red lines. The gray zone in the top left corner is the parameter range, where there are no complex eigenvalues at all -neither with stable eigenmodes, nor unstable ones.

PULSATING HEAT PIPE WITH EVAPORATION
Among the assumptions made and the constraints imposed in Section 2, the first to be relaxed is that the working medium cannot change its phase. This is the most important step, since practically there are no PHPs with non-condensible gas bubbles. The simplistic example in Section 2 merely is the simplest way to reproduce the two characteristic instabilities of the PHP and served top introduce our way of analysis.
The main difference between the PHP with non-condensible gas bubbles and the one that allows evaporation and condensation is the way heat transfer influences the dynamics. In case of non-condensible bubbles, the pressure in a bubble can only be changed by compression or heat addition to the gas. Additional heat can be transferred from the wall into the slug, but the slug temperature does not influence the PHP dynamics.
When allowing phase change at the interface between bubble and slug, heat transfer from wall to slug can be used for phase change (Fig. 8).
Since the density of the vapor is orders of magnitude lower than the density of the liquid, it can be expected that this mechanism -when present -is the dominant one. This justifies an important assumption that we can make at this point: that the heat transfer between bubble and wall is negligible compared to the influence of phase change.
As a linear stability analysis, it regards infinitesimal perturbations around an equilibrium state. In the equilibrium state, no net heat transfer into or out of the slug may take place since that would result in a longterm change in temperature and eventually phase change, i.e. change of mass. Ultimately relevant for the analysis is therefore merely the perturbation of the heat exchange between slug and wall. Because the perturbation of the heat flux must first propagate to the phase interface to become relevant for the PHP dynamics, and because typical slugs are clearly longer than the tube diameter, the following analysis considers only the heat flux between wall and slug that takes place in the vicinity of a phase interface.
In particular once the slug is moving, the hydrodynamic boundary layer close to the phase interface is very thin (for both advancing and receding phase interface, see (Srinivasan and Khandekar, 2017;Thulasidas et al., 1997) for details on the flow pattern), which favors the effect of phase change even more over direct heat transfer from the wall into the bubble (Srinivasan and Khandekar, 2017;Thulasidas et al., 1997). Retrospectively, it becomes clear that in Section 2, the direct heat transfer between bubble and wall was only used, because it is the only relevant wall heat transfer term under the conditions given.
Assuming that the heat transfer between wall and phase interface is the only relevant one, brings another advantage: The fact that all heat transfer takes place at the same axial position eliminates the integration over the wall temperature, i.e. one uses a local value Tw instead of a spatial averageTw.
The governing equations are quite similar to the ones from Section 2: Note that the momentum balance for the liquid slug does not change despite the fact that it is now an open system: Convective momentum transfer cancels with the time derivative of mass. Instead of the -now negligible -heat transfer term, there is a convection term that transports enthalpy of saturated vapor (hv,j with j = 1, 2) into the bubbles at the respective mass flow ratesṁj. The mass flow rate of saturated vapor added to a bubble can be computed from the heat flow rate at the phase interface:ṁ where ∆h lv is the latent heat of the working medium under equilibrium conditions. As in Section 2, the heat transfer is modeled with a constant K, which incorporates heat transfer coefficient as well as reference area 6 . The temperature at the phase interface must be saturation temperature Tsat (see Fig. 8 Note that distinct from e.g. Dobson (2004), this study regards a PHP with dry walls along the bubbles (no liquid film). Therefore the heat transfer takes place at the phase interface only and the relevant temperature difference is the one between the (local, not mean) wall temperature at the (a) Colored eigenvalues for τ + = 1.
phase interface (Tw(x) instead ofTw(x)) 7 and the saturation temperature (not the vapor temperature in the bubble T1 or T2, which are in general different from Tsat). For an investigation of start-up, the assumption of dry walls makes sense, because initially the PHP is in equilibrium state, and heated walls will dry out eventually, if there is no slug motion. In addition this avoids the trouble of properly modeling the film thickness. Clausius-Clapeyron's relation (for ideal gases) allows to express temperature perturbations T sat from the equilibrium state, where Tsat = Teq: where the ∓ refers to left (j = 1) and right (j = 2) bubble, respectively. 7 see Eq. 4 for definition of the average and its dependence from the slug displacement x.
Also the above equation defines the modified Stefan number: with the specific enthalpy heq and isobaric heat capacity cp in a bubble at equilibrium state. Now one proceeds with the same strategy as in Section 2. First insert and linearize in equilibrium state, where x = 0, peq ≡ p1 = p2, Teq ≡ T1 = T2 = Tw1 = Tw2, and ∂Tw/∂x ≡ ∂Tw1/∂x = −∂Tw2/∂x. Then eliminate p and V as for non-condensible gas case: Finally, non-dimensionalize the same way as the non-condensible case: with the characteristic equation: which is similar to the one for non-condensible gas, but has different coefficients. Again, one stable eigenvalue is found immediately, so that ∂T Physically, the criterion for an oscillatory instability reflects the same mechanism as in the non-condensing case. The scaling factor Ste < 1 reflects that the saturation temperature (which is now the relevant temperature for heat transfer) is less pressure-sensitive than the bubble temperature. This makes the system more prone to instability. Compared to the case with non-condensible gas bubbles (Section 2), the non-oscillatory instability mechanism is modified as follows: The condition for the existence of the instability is now: i.e., when there is a very small, but positive temperature gradient along the wall of the left bubble (and a negative one on the wall of the right bubble, for symmetry reasons) and the bubble position is perturbed slightly to the right, the wall temperature at the left phase interface increases a little, whereas the left saturation temperature decreases due to the pressure reduction from expansion. Now, the liquid at the left phase interface will start evaporating until pressure is high enough that the saturation temperature in the bubble is the same as the wall temperature again. The opposite mechanism takes place in the right bubble. In the end the left bubble has a higher pressure (due to higher wall temperature at the phase interface) than the right bubble, which furthers the movement of the bubble to the right.
Other than for non-condensible gas bubbles, this effect can be sustained until the complete collapse of one bubble, and it does not require extraordinary circumstances, merely that the slug is placed in a hot zone of the wall. Although this type of instability can cause a slug to accelerate significantly (and contribute to the increase of kinetic energy in the PHP), it is questionable whether it can be stated that this instability favors PHP start-up, since it ultimately converges to a more stable state than the original one.

SURFACE TENSION
It is well known that effects related to surface tension play a role in PHP operation, at least by preventing stratified flow and enforcing the bubbleslug pattern, which is vital for the working principle of the PHP. There have been experimental (Srinivasan and Khandekar, 2017) and numerical (Dilawar and Pattamatta, 2013;Srinivasan and Khandekar, 2017) studies on the role of surface tension in PHP operation. Therefore this section shall discuss such effects, based on the minimalist PHP model with phase change as introduced in Section 3.
The fundamental equations are the same, except that now there is a distinction between the pressures on either side of the phase interface. p1 and p2 refer to the pressure in the bubbles, i.e. on the vapor side of the interface. Newly introduced are the pressures on the liquid side p l,1 and p l,2 : The only equation affected is the momentum balance. One can express the pressure in the liquid adjacent to bubble j by the Young-Laplace equation: where σ lv,j is the surface tension of the liquid-vapor interface, and Θj is the contact angle at the contact line of the three phases, see Fig. 9. Note that giving an index j to both variables is required to describe their variation, in particular a temperature dependence and contact angle hysteresis. For a PHP made of a round tube, d is just its diameter, and therefore constant (omit the index j), however, for other cross-section geometries the situation can be more complicated. Since it will not influence the core message of this paper, the following analysis is restricted to a cylindrical geometry.

Contact Angle Hysteresis
The term contact angle hysteresis (CAH) refers to the phenomenon that the contact angle Θa of an advancing contact line (i.e. the liquid covers more and more of the solid substrate) is greater than the contact angle Θr of a receding contact line. For a static contact line, the contact angle may assume any value between the two limiting cases Θr ≤ Θ ≤ Θa. The contact angles of moving contact lines also depend on the velocity of the contact line, however, the velocity influence is minor (Eral et al., 2013) and will therefore be neglected in the following analysis.
There is a finite pressure difference that has to be overcome in order to make a static slug move. Since this mimics the physical behavior of kinetic dry friction, it will be referred to as frictional pressure difference in the following.
Now rewrite the momentum balance: When defining ∆pfric, the distinction between the two values of surface tension was kept up, but they were assigned to the direction of contact line motion σ lv,r and σ lv,a , instead of σ lv,1 and σ lv,2 , see Fig. 10. This allows to take into account the temperature difference influencing the surface tension, while still exploiting the symmetry of the problem. Θr, Θa, and σ lv are assumed to be temperature dependent. The temperature of relevance here is the temperature of the phase interface, which is saturation temperature Tsat, which again depends on pressure. This leads to the concept of boiling delay, i.e. superheating the liquid phase beyond the regular point of evaporation. For at least partially wetting surfaces (contact angle Θ < π/2) the pressure in the bubble is greater than in the adjacent liquid phase. Therefore the equilibrium saturation temperature in the bubble is also greater and defining the reference temperature for contact angles and surface tension proves to be cumbersome. The following discussion will neglect this effect and instead assume that contact angles and surface tension depend on the bubble saturation temperature, i.e. the bubble pressure. This decision does not affect the result, which is of qualitative nature anyway. Concluding one can state: and similarly with the index k = r, a for receding and advancing contact line, respectively. Again, one can non-dimensionalize the momentum equation as in Sections 2 and 3, and obtain a factor for the frictional pressure ∆pfric, which defines the non-dimensional frictional pressure ∆p * fric : introducing a non-dimensional the surface tension σ * r/a = σ lv,r/a /σ lv,eq with the surface tension σ lv,eq at equilibrium temperature, and a nondimensional constant: One can linearize the frictional pressure by distinguishing between the fluctuating (pressure dependent) component ∆p + fric and the one that is constant ∆p * fric .
Note that the value of the non-dimensional surface tension is unity in the central position and that T + sat = (γ − 1) /γ Ste p + expresses the non-dimensional saturation temperature fluctuation (see Section 3). The pressure fluctuation is again defined as: with the case w + = 0 undefined. It is impossible to linearize around the static central position in this case, since the right hand side of the momentum balance does not have a unique differential for w = 0. It is however possible to linearize by regarding the different cases individually. Subsequently, introduce the same non-dimensional variables as before.
The following two abbreviations allow convenient presentation: Although previous studies present some approaches (Carey, 1992;Pomeau and Vannimenus, 1985;Joanny and de Gennes, 1984), it is very hard to make a general statement on the temperature dependence of CAH ∂Θr/∂T * , ∂Θa/∂T * . For small perturbations around the static central position, ∆p * fric is always the dominating term in the momentum balance. One can hence integrate it immediately, which then also allows to integrate the kinematic constraint.
where x + 0 = x + t + = 0 and w + 0 = w + t + = 0 are the initial conditions. As expected the movement dies down after a short interval of time t + end = w + 0 /∆p * fric . Regarding the behavior of the system after the initial perturbation has died down t + > t + end , the system has been separated into two independent systems of first order: Both ODEs are damped, thus one concludes that CAH makes the PHP model completely stable against infinitesimal perturbations. Theoretically the most asymmetric stable positions can be determined by p1 − p2 = ∆pfric. However the temperature dependence of CAH is too unclear to do this in practice. This study will therefore refrain from investigating stability in such a position. The result that PHP start-up is impossible in the presence of CAH stands in contradiction to experimental results of PHPs actually working (Quan and Jia, 2009;Song and Xu, 2009;Qu et al., 2009;Khandekar and Groll, 2004). This discrepancy of course requires further discussion, which can be found in Section 5.4.

High Surface Energy Materials and Surfactants
The above analysis suggests that CAH is unfavorable for PHP operation. It is therefore a worthy approach to design a PHP with the goal to minimize or even eliminate CAH. Smith et al. (2014) find experimentally that the refrigerants R-134a and HFO-1234yf are a lot less prone to contact angle hysteresis than classical fluids (water and acetone were investigated). Whereas this is a highly valuable insight, it would be very tedious to systematically investigate a broader range of material combinations experimentally. It follows a small literature survey with the goal to find a more general lever to reduce the influence of CAH.
In textbooks (Carey, 1992) on the matter, one finds that the contact angle itself is the result of the system minimizing its Helmholtz free energy by making a compromise between the three types of interfaces, which each have a surface tension (solid-liquid σ sl , solid-vapor σsv, liquid-vapor σ lv ). Young's law summarizes this compromise (see Fig. 11 This relation is valid as long as σsv ≤ σ sl + σ lv . Otherwise the contact angle is zero Θ = 0, corresponding to full wetting. The qualitative phenomenon of CAH is explained in Carey (1992) by inhomogeneities of the surface, like stains, heterogeneous materials (e.g. multiphase alloys and imperfections in the crystal structure), and wall roughness. All these imperfections represent perturbations of the surface free energy of the solid, and consequently the theoretical contact angle. Eral et al. (2013) provide an excellent review on the topic. Among other sources, they cite attempts of a mathematical description of aforementioned qualitative explanation (most notably Pomeau and Vannimenus, 1985;Joanny and de Gennes, 1984). Also models for the velocity dependence of CAH based on molecular dynamics as well as continuum assumption can be found here (Eral et al., 2013), none of which however works without fitting parameters.
In combination these studies allow the conclusion that CAH can be completely eliminated by finding a material with a sufficiently high surface free energy that perturbations of the latter do not result in the solid-vapor surface tension decreasing below the sum of solid-liquid and liquid-vapor surface tension: i.e. what is needed is a solid substrate that is fully wetted by the working fluid of the PHP. The substrate must be smooth and very clean, the latter of which is difficult to achieve, since a stain will decrease the free energy of the system substantially and is therefore a highly favorable state from a thermodynamic point of view. Note that what is actually required is a substrate-gas combination with high surface tension σsv, not necessarily a solid material with a high surface free energy σs. However the nature of surface tension -resulting from Lennard-Jones potential, i.e. short-range repelling and long-range attracting forces between molecules -suggests that the surface tension between two partners with very different densities (like solid-gas and liquid-gas) is dominated by the high-density partner. The comparison of experimental values for the surface tension of water in its vapor (Vargaftik et al., 1983) and water in air at 100% relative humidity (Pérez-Díaz et al., 2012) shows, that in this case the difference for a given temperature is lower than measurement precision, which confirms the initial idea. Finally, the determination of surface free energies of solids by the Zisman method is based on the approximation that the surface free energy of a solid is the same as the surface tension between solid and vapor (Carey, 1992;Shafrin and Zisman, 1972). High surface energy materials (metals, glass) typically have surface energies σsv > 500 mN/m -as opposed to low surface energy materials, like polymers (also acrylic), which often range between 15 mN/m < σsv < 40 mN/m (Carey, 1992). To compare, the surface tension of non-metallic liquids is mostly between 15 mN/m < σsv < 75 mN/m (Carey, 1992). Material data can be found in (Carey, 1992;Shafrin and Zisman, 1972;Eustathopoulos et al., 1999;Weirauch and Ownby, 1999).
Also according to Young's law, σsv σ lv is not a sufficient condition for full wetting, but only an indicator. It is also necessary that the surface tension between liquid and solid is sufficiently low, which according to the data given by Carey (1992) should often be the case with high surface energy materials. As examples, Carey mentions low-viscosity silicone oils on metal, or exotic combinations like liquid helium on glass.
An alternative approach would be to add a surfactant to the liquid, as has been experimentally demonstrated by Srinivasan et al. (2015). A surfactant accumulates at the phase interfaces on which it has an energetically favorable effect, i.e. it reduces σ lv and possibly σ sl , whereas it does not affect σsv. The selection of a surfactant is not trivial because of its evaporation properties. Too low saturation temperature (compared to the bulk of the working medium) may have an adverse effect on the concentration accumulated at the respective phase interfaces. A surfactant with a high saturation temperature could however leave a stain on the solid surface after the liquid film dries out, effectively reducing σ sl .
Assume there is a material combination that completely avoids CAH. The contact angle will be Θ = 0. In this case the three constants Cr, Ca, and ∆p * fric are: Therefore the fundamental equations of minimalist PHP model will have the following non-dimensional, linearized form: Since there are no constant terms and linearization is again possible, one can interpret this result by similar means as in Sections 2 and 3. The characteristic polynomial is the same as in Section 3, except the factor Cr = Ca in front of the first order and constant terms: Again, there is the same stable non-oscillatory eigenvalue as in Section 3. Looking for the stability margin, i.e. λ + = iω + , one finds that the temperature dependence of the surface tension merely influences the oscillation frequency of the oscillatory eigenvalue, but leaves the stability limits untouched: In hindsight this was to be expected, as surface tension, as well as its temperature dependence, merely poses a modification to the spring constant of the problem, but not to the actual mechanism of instability. Note that it is not necessary to investigate the case of Cr ≤ 0, since ∂σ * /∂T * < 0.

DISCUSSION
This section reviews the results presented in the previous sections with respect to experimental and numerical findings from literature.
5.1. Non-condensible gas PHPs require prohibitively high wall temperature gradients As already stated in Section 2, the example with bubbles from a noncondensible gas merely served to present the layout and working principle of our PHP model, as well as to introduce linear stability analysis and the two types of instability that can occur. The authors are not aware of an experimental PHP with this feature. There are however studies investigating the presence of non-condensible gases in the PHP working medium (Quan and Jia, 2009;Jia and Yin, 2007), i.e. bubbles can never be fully condensed. These studies come to the conclusion that the presence of non-condensible gases is unfavorable for PHP operation. The present study can reproduce this trend by comparing the results for the non-condensible gas (Section 2) and the evaporator / condenser case (Section 3). One difference between the cases is the relevant wall temperature. For the non-condensible gas case, this is the mean value over the entire wall in contact with the bubble: whereas for the evaporator / condenser case it is the wall temperature T * w at the phase interface in position x * . Note that whereas the slug position is denoted by x * , the axial coordinate for the integration of the wall temperature isx * . The exact relation between the two gradients ∂T * w /∂x * and ∂T * w /∂x * therefore depends on the wall temperature profile T * w (x * ). For a linear profile of T * w (x * ), one finds the simple relation: The other difference is the start-up criterion itself, involving the modified Stefan number for a phase-change working fluid. For water at Teq = 300 K, the modified Stefan number is Ste = 0.236. Therefore, a temperature difference of −23.6 K is required over the bubble length to initiate the first oscillatory instability in an evaporator / condenser PHP filled with water, whereas it takes a much larger temperature difference of −240 K for the same instability to occur in a PHP with air bubbles. Since this would result in a freezing water slug (T w,slug ≤ 300 K + −240 K/2 = 180 K), one would have to increase the equilibrium temperature to Teq = 460 K, resulting in a required temperature difference of −368 K in order to avoid freezing the slug. To achieve a non-oscillatory instability with a non-condensing gas PHP, the model presented requires a temperature difference of +600 K at an equilibrium temperature of Teq = 300 K, whereas for the evaporator / condenser PHP the non-oscillatory instability is present for all non-negative temperature gradients. The unrealistic numbers for the noncondensing gas case lead to the conclusion that the model -although simplistic -confirms that phase change of the working medium is crucial for PHP operation. Shafii et al. (2001) found that their computational model always converges to a number of bubbles equal to the number of evaporator sections. This can be explained by the non-oscillatory instability found in this paper: It is impossible to conceive a distribution of slugs and bubbles in a PHP that has more bubbles than evaporator sections and at the same time does not involve a non-oscillatory instability. One can rephrase the argument in physical terms: In its essence the mechanism of the non-oscillatory instability (given phase change is allowed) is that as long as the wall temperatures at the two ends of a slug differ, the bubble at the warmer end grows at cost of the other one. Therefore, a slug in the condenser section is stable: If moved to one side, the wall temperatures at its ends induce a restoring force, see Fig. 1, assuming a continuous wall temperature distribution. Contrarily, a slug in the evaporator is unstable, i.e. it will eventually move towards one of the two adjacent condensers. A slug between evaporator and condenser will continue its motion towards the condenser, thereby shrinking the bubble on the condenser side till collapse.

Number of bubbles and non-oscillatory instability
In experimental PHPs, convergence towards one bubble per evaporator section is not observed in general (Khandekar and Groll, 2004), which Shafii et al. (2001) explain by bubble formation in the evaporator -a physical process that was not taken into account in their model. In order to model the dynamics of a PHP properly, it is necessary to improve the understanding and models of bubble formation.

Influence of the filling ratio
The influence of the filling ratio on PHP performance, and in particular on start-up reliability has been investigated repeatedly in the past (Han et al., 2016). Results are contradictory. Most studies find the best performance i.e. highest heat conductivity, of the PHP at a filling ratio around 50%, there are however strong deviations, e.g. Khandekar et al. (2003) find the best performance at 10%.
In the models investigated, the filling ratio occurs implicitly in the constant τ + , since it relates the mass of the slug m l with the bubble length L. In none of the cases regarded has τ + an effect on whether startup occurs or not. It does, however, affect the growth rate as well as the frequency of instabilities with significant growth rate.
This shows the limitations of the extremely reduced setup used for this study. In theory the finding that filling ratio does not influence startup behavior may be correct, however in a practical PHP, the filling ratio influences the number of phase interfaces. Since capillary forces define a minimum size of bubbles and slugs, a filling ratio of about 50% allows the highest amount of phase interfaces. Since start-up by linear instability requires a phase interface in a zone with a wall temperature gradient, the start-up probability is higher with more phase interfaces. This argument stands in agreement with Stevens et al. (2019), who performed an extensive investigation on reproducibility of PHP start-up and operation and found that there may be significant differences in start-up reliability between two identically built PHPs.

Suppression of instabilities by contact angle hysteresis
Contact angle hysteresis makes the PHP model under investigation completely stable against infinitesimal perturbations (see also Srinivasan et al., 2015;Srinivasan and Khandekar, 2017). This result stands in contradiction to experimental results of PHPs actually working (Quan and Jia, 2009;Song and Xu, 2009;Qu et al., 2009;Stevens et al., 2019;Khandekar and Groll, 2004;Han et al., 2016). Nevertheless, literature shows quite a lot of evidence for the adverse effect CAH has on PHP operation.
In a real PHP there are many more bubbles and slugs than in the models investigated. Also, typically at least on the hot part of the wall there is no boundary condition of fixed temperature, but rather of fixed heat flux, since the main prospected application of the PHP is cooling of micro-electronics. As long as the PHP does not contribute to heat transfer the evaporator wall temperature will increase. If the initial bubble distribution is asymmetric, a pressure difference sufficient to overcome frictional pressure is eventually built up causing a first movement. Such a situation should not be confused with a linear instability 8 . Contact angle hysteresis creates a finite range of stable bubble positions, which also significantly increases the probability of so-called stop-over, i.e. the prob-8 It could however pave the way to a linear instability by laying a film on the bubble walls (see Section 5.5). A film on the walls would eliminate CAH and its stabilizing effect. ability of the PHP coincidentally finding a state that allows to cease all motion.
Most experimental setups of PHPs use high surface energy materials (find a list in Han et al. (2016)), such as copper, aluminum or glass. Such materials reduce the relevance of CAH, and thus the visibility of its adverse effect. Using a PHP made of two acrylic plates -a material of low surface free energy - Khandekar et al. (2003) obtained poor performance in horizontal operation, i.e. when the inherent instability resulting from bottom-heat mode is missing. They also concluded that part of the reason must be CAH, and they had a validation experiment with glass tubes, but did not report quantitative results for the latter, since such results would not be comparable due to a change in geometry.
A mechanism that would circumvent the barrier posed by CAH is bubble formation within a slug in the evaporator zone. A newly formed bubble initially forms a liquid film between gaseous and solid phase. The film eliminates CAH and hence enables the instability mechanisms presented in this paper. Due to its non-linear nature, bubble formation could not be taken into account in the models of this paper. As discussed before, bubble formation is a mechanism that is not sufficiently understood and requires more research.
Contradicting all the evidence collected above, Dilawar and Pattamatta (2013) found in their numerical study that CAH does not have a significant influence on PHP performance in steady-state operation. They argue that the pressure differences caused are small compared to the pressure differences resulting from evaporation and condensation at different (wall-)temperatures. The present study has not tried to reproduce this result, but from Srinivasan et al. (2015), as an exemplary value it is known that in case of a water slug in a 1.5 mm glass tube the pressure difference from CAH was 240 Pa, which is clearly more than what would be considered an infinitesimal perturbation.
In conclusion, CAH has an adverse effect on PHP start-up, which can be mitigated by use of surfactants (Srinivasan et al., 2015) and appropriate choice of materials. It is apparent from the discussion that CAH in PHPs is not fully understood and requires further research.

Film models and dry-out
In particular because it avoids the adverse effect of CAH on PHP startup, but also because of the high area available for evaporation, a liquid film between bubble and wall is favorable for PHP operation and startup. In contrast, this paper assumes a dry wall throughout. Of course, partially this is to not further complicate the model, and avoid the need to arbitrarily assume a film thickness (Dobson, 2004).
The authors also think that a dry wall is a good assumption for a start-up scenario, since Khandekar and Groll (2004) found that agglomeration of all bubbles in the evaporator is an important precursor to stopover, i.e. the complete cease of motion in the PHP. This situation requires the dry-out of the evaporator wall, since a film would still constitute evaporation, i.e. motion. Stop-over then is the initial condition for start-up.
As a criterion for the existence of a liquid film in partially wetting conditions, researchers have often resorted to the capillary number Ca = η l w/σ lv > 10 −3 , with the dynamic viscosity η l of the liquid (Eral et al., 2013;Srinivasan et al., 2015;Srinivasan and Khandekar, 2017). Since for start-up, the initial condition is zero velocity w = 0, the criterion for a liquid film to cover the wall is never fulfilled.

CONCLUSIONS
From the theoretical analysis of linear stability of a minimum working example of a PHP we conclude the following: Firstly, evaporation and condensation are crucial for the start-up of a PHP. Whereas this has been understood for a while, the authors are not aware of another publication deriving this fact from first principles. This paper has demonstrated that PHP operation is actually not impossible with a non-condensible gas, but that the required temperature differences are prohibitively high.
Secondly, there are two types of instability: an oscillatory one and a non-oscillatory one. The oscillatory instability follows the pattern of a thermodynamic cycle. Given a small perturbation that reduces the size of a bubble, its pressure increases. If by the same perturbation the bubble moves to a higher relevant wall temperature, heat is added, increasing the pressure even more. Eventually the pressure increase reverses the slug movement, resulting in bubble expansion, pressure decrease, and -due to lower wall temperature -heat rejection at low pressure. The work resulting from this cycle is invested in kinetic energy of the slug, i.e. higher and higher oscillation amplitudes. This mechanism is the same as in a thermo-acoustic instability.
The non-oscillatory instability occurs in case of the opposite temperature gradient. When the bubble volume is slightly reduced, the relevant wall temperature drops at the same time. A sufficient amount of heat is removed from the bubble to decrease its pressure despite volume reduction. Remarkably, in case of evaporation / condensation the non-oscillatory instability has an onset at zero temperature gradient.
Finally, contact angle hysteresis suppresses both types of instability. For PHP design it is hence favorable to choose high surface energy materials and make sure they are smooth and very clean. With contact angle hysteresis present, the start-up mechanism of a PHP cannot be explained by linear instability. Instead the PHP must rely on a non-linear mechanism, like bubble formation -crucial for start-up -or film evaporation. If the influence of contact angle hysteresis can be eliminated, the surface tension itself and its temperature dependence do not affect start-up reliability.