Interfacial instability as a trigger for dryout inception in two-phase CO2 flow

Interf a cial inst ability as a trigger f or dr yout inception in tw o-phase CO 2 flo w A Preprint G. Cantini Departmen t of Mec hanical Engineering Univ ersity of Bath Cla verton Do wn, Bath BA2 7A Y United Kingdom gc803@bath.ac.uk G. Arnone Istituto Nazionale di Alta Matematica F. Sev eri INdAM, Piazzale Aldo Moro 5, Roma, 00185 Italy arnone@altamatematica.it F. Cap one Dipartimen to di Matematica e Applicazioni R. Cacciopp oli Univ ersity of Naples F ederico I I Nap oli, 80126 Italy fcapone@unina.it J.A. Gianfrani Dipartimen to di Matematica e Applicazioni R. Cacciopp oli Univ ersity of Naples F ederico I I Nap oli, 80126 Italy jacopoalfonso.gianfrani@unina.it M. Carnev ale Departmen t of Engineering and Applied Sciences Univ ersity of Bergamo Via Marconi 4, 24044 Dalmine (BG) Italy mauro.carnevale@unibg.it Marc h 25, 2026 Abstract Progress in particle physic leads to increasing in detector luminosity and a consequen t in- creasing o verheating induced b y Joule effect. An effectiv e co oling strategy is the exploitation of CO 2 heat latency in phase-c hange. An additional c hallenge, relev ant to detectors for High Energy Particles, is the consequent geometrical constrain due to the limited space avialable for the co oling system within the detector arrangement, leading to the implemen tation of co oling system by means of millichannels. In this context, at relative high v ap our quality the liquid phase exhibits annular flo w, an ticipating the dry out. Dryout is a critical condition where the heat transfer co efficien t dramatically drops and dangerous temp erature lev els can b e reached, p otentially leading to catastrophic consequences. Exp erimental evidences reveal that its b eha vior in tw o-phase annular flows differs from conv en tional refrigeran ts and the fundamen tal inception-mechanism is not yet understo od. This study aims at inv estigating the k ey new idea whereby dryout inception is triggered b y instability of the liquid-v ap our in terface. A mathematical mo del for t wo-phase ann ular flow is presen ted and the stability of the in terface b etw een the tw o fluids is studied through the linear theory . The stability anal- ysis reduces to solving a coupled forth-order differential eigenv alue problem that is treated n umerically with an in-house code based on the Cheb yshev- τ metho d. Numerical inv esti- gations identify a critical v alue for the v ap our quality , named x dry , that leads to interface instabilit y . The resulting predictions on x dry are confirmed by exp erimental data collected from t wo indep endent experimental campaigns, v alidating the hypothesis that dryout in- ception is gov erned by in terfacial instabilities. A Preprint 1 In tro duction The dryout phenomenon is known in heat transfer pro cesses as a condition to b e a voided in order to keep a co oling system working at its b est [ 1 ]. It is defined as the condition where the Heat T ransfer Co efficien t (HTC) drops drastically within the heat exchanger, reac hing v alues comparable to those of gaseous flows [ 2 ]. In heat managemen t processes, during the ev ap oration the flo w is bi-phase and the heat transfer is more efficien t than in flows without phase c hange [ 3 ]. Nonetheless, as the v apour quality increases along the channel, the progressive depletion of the wall film may trigger dryout. In suc h a condition, the liquid film detaches from the pip e wall that is then in direct contact with the v ap our phase, leading to an abrupt degradation of the local HTC and p otentially compromising component integrit y . In applications where there is a constant imp osed heat flux at the wall, like in n uclear reactor or in electronics co oling, the dryout can raise the wall temp erature enough to melt the exc hanger bringing the devices to failure or dangerous consequences [ 4 ]. Therefore, a more accurate prediction of the onset of dryout can b e useful to supp ort less conserv ative safet y margin in co oling system design, with a more efficient result. Among all refrigerants, b oiling carbon dio xide offers numerous adv antages ov er synthetic refrigerants. Its cost, reliability and minimal Global W arming P oten tial mak e CO 2 c heap, efficient and en vironmentally friendly . The high w orking pressure and low liquid-to-v ap our density ratio of CO 2 mak e it w ell-suited for efficien t heat exc hange in milli-scale diameter pipes and for saturation temp eratures T sat in the interv al ( − 40 ◦ C, − 10 ◦ C ) [ 5 ]. This is b ecause the ev ap oration tak es place at higher pressure than standard refrig- eran ts and the pressure drop has a limited relativ e effect on b oiling pressure compared to standard fluids [ 6 ]. Additionally , CO 2 is neither toxic nor flammable, and its non-corrosive prop erties are piv otal for safe op eration of the electronics and other equipmen t. Its resistance to radiation makes it an ideal candidate for co oling applications in particle physics detectors. F urthermore, for small-diameter channels the thermal conductivit y of CO 2 is higher than standard fluids as can b e seen in Figure 2 in [ 6 ]. The ab o ve features, how ev er, make CO 2 t wo-phase flows qualitatively different from those of conv en tional refrigeran ts. As a result, phenomenological mo dels struggle to accurately predict the v ap our quality at the onset of dryout, from now on denoted as x dry . Accurate estimation of x dry is therefore pivotal in co oling- system design. In [ 7 ], tw o trends ha ve been identified for the CO 2 dry out v ap our qualit y as a function of the mass flux: • δ + regime: x dry increases with rising mass flux G . Although less frequently observed in the literature, some instances of this b ehaviour are do cumen ted in [ 8 ] and [9]. • δ − regime: x dry decreases with increasing mass flux G . This regime is widely reported in the literature as in [10] and [11]. The flow pattern map prop osed by [11] is considered the most robust for CO 2 but represents solely δ − b e- ha viour. Correlations informed b y the δ − regime will result in an underestimation of x dry when δ + b eha viour is present. An attempt to describ e the transition from δ + to δ − is prop osed in [ 7 ] where a phenomenological mo del based on a system of differential equations for liquid and v ap our phases is introduced. Ho w ever, such a mo del disregards the w all heat flux, which has b een shown to b e imp ortant in the δ + regime through an exp erimen tal study conducted in [ 9 ]. Imp ortan tly , all theoretical approaches presented in the open literature fail to mo del the comp ound effects of all contributing physical parameters, as stressed by [12]. A recen t paper [13] provides an initial step tow ard filling this gap. This research is based on an experimental campaign conducted at CERN, and is in tended to in vestigate the dryout inception in CO 2 in milliscale pip es under representativ e op erating conditions for high-heat-flux applications in the δ + regime. In [13] the authors show the inability of existing mo dels to predict dryout inception across all conditions, particularly the less understoo d δ + regime where dry out v ap our qualit y increases with mass flux. Moreo ver, as a nov elt y , the authors connect exp erimen tal observ ations to stability theory through dimensional analysis, deriving a dryout instability factor , denoted as I δ + , that naturally emerges from momentum balance at the interface. The empirical correlation b etw een I δ + and measured x dry bridges first-principles theory and data-driv en prediction, v alidated against independent data from [ 9 ]. Remarkably , this framework is indep endent of the saturation temp erature and the heat flux, namely it reduces the dep endence on exp erimen t-sp ecific conditions, making design and extrap olation easier. In this con text, motiv ated by [13], the present work aims to supp ort the h yp othesis that the dryout phe- nomenon is triggered b y instability of the interface b et ween liquid and gaseous CO 2 , pro viding a mathematical theoretical coun terpart of the exp erimen tal empirical campaign. T o this end, w e dev elop a mathematical 2 A Preprint mo del for tw o-phase annular flow, with liquid and v ap our separated by a sharp interface with appropriately defined interfacial conditions. The present study represents the first application of rigorous stability theory to the δ + regime in boiling CO 2 ann ular flow, pro viding a framew ork that sim ultaneously accounts for all relev ant mechanisms: viscous effects, surface tension, heat and mass transfer, and comp ound effects of geometry and op erating conditions. The resulting equations form a differential eigenv alue problem of a couple of fourth-order differential equations requiring eight b oundary conditions: four at the interface, t wo at the pip e centerline, and tw o at the w all. W e show that, for a prescrib ed set of parameters, solving the eigenv alue problem bridges stability c harac- teristics with x dry , b y identifying dry out inception when the real part of the leading eigenv alue b ecomes p ositiv e. Numerical results show go od agreement with exp erimen tal data from tw o indep endent campaigns [ 9 , 13]. At this point, a natural question arises: is CO 2 the only fluid suitable to the present study , or does it realize, under practical op erating conditions, the thermodynamic state required for the onset of the δ + regime? The latter is the answ er. Sp ecifically , compared to most refrigerants, CO 2 exhibits muc h smaller liquid-to-v ap our density ratio (i.e. more similar liquid and v ap our densities) at typical refrigeration op er- ating temp eratures. A relatively small densit y ratio implies a reduced slip ratio b etw een v ap our core and liquid film, hence weak er interfacial shear and smo other interface in annular flow. These are precisely the conditions under whic h the interfacial-stabilit y mechanism considered here becomes ph ysically meaningful and predictive for dryout inception. Conv ersely , for standard refrigerants at the same temp eratures, the m uch larger density ratio promotes δ − b eha viour. Actually , only outside the standard op erating conditions relev ant for ev ap orators (e.g. saturation temp eratures far abov e practical refrigeration op eration), the δ + regime dominates for standard refrigerants. The here prop osed mo del implicitly contains a hypothesis that the δ + regime might exist for other fluids at some conditions, but those conditions are not practically rele- v ant. Only CO 2 can achiev e sp ecific thermo dynamic state at practical op erating conditions required for δ + b eha viour. The reader is referred to Section 5 for details. In literature, in terfacial instability b et ween tw o sup erp osed fluids has b een widely analyzed [14, 15, 16, 17, 18]. In particular, [17] prop osed a model for interfacial instability betw een t wo sup erp osed homogeneous compressible fluids at rest that are inv ested by thermal gradien ts. These fluids occupy a horizontal lay er and the interface b et ween them (identified by a horizontal profile) divides the region in t wo parts and it is treated as a sharp in terface. Therefore, in terfacial conditions are appended to the mo del and the effect of ev ap oration and interfacial wa ves on the onset of Ra yleigh-T a ylor instability is studied. Later on, the same author prop osed in [18] a simplified mathematical mo del of the same physical setting where thermal effects mo dify the o ccurrence of interfacial instability . The main assumption that simplifies the mo del regards the replacemen t of the energy balance equation with an equation of pro duction of mass at the interface. This is relying on the fact that sinc e the tr ansfer of mass acr oss the interfac e r epr esents a tr ansformation of the fluid fr om one phase to another, ther e is invariably a latent he at asso ciate d with the phase change . The author then remarks that it is essential ly thr ough this interfacial c oupling b etwe en the mass tr ansfer and the r ele ase of latent he at that the motion of fluids is influenc e d by the thermal effe cts [18]. In other words, the net heat flux across the in terface is mo deled as a function of interface displacement. In the present study , w e follow this idea. Unlik e the ab o ve-men tioned works, the physical framew ork examined here in v olves heat flux at the pipe wall that provides external heat to the system and cylindrical geometry constraint. A fundamental assumption of this researc h is to mo del the presence of the external heat flux via a non-constant profile of the interface o ver the axial length of the pip e. F or this reason, the liquid film thic kness is assumed to be reducing along the axis of the pip e, so as to mo del the heat and mass transfer from the liquid phase to the v apour one. Moreo v er, despite the mo del b eing defined up on tw o-dimensional approximation, it keeps track of the cylindrical geometry within the interfacial conditions, via the dryout instability factor I δ + , defined up on relations that hold only in cylindrical geometries, see Eq. ( 27). In this regard, let us remark that the δ + regime in CO 2 t wo-phase flow can b e achiev ed only in millichannels where the annular flow o ccurs as larger geometries w ould promote gra vitational stratification of the flow, p oten tially shifting behavior to w ard δ − ev en when density ratios remain fa v orable. Studies using larger c hannels ma y not access the δ + regime despite using CO 2 . The plan of the paper is the following. In Section 2, the physical problem is introduced via a schematic framew ork with the aim of clarifying the geometry of the problem and in tro ducing the relev ant v ariables app earing in the mathematical mo del. In Section 3 , the deriv ation of the mathematical model and interfacial conditions is pro vided, including the main assumptions that lead to a mathematically treatable problem. The dimensional analysis is shown in Section 4 , where the nondimensional dryout instability factor is defined. 3 A Preprint Section 5 deals with the rationale for considering CO 2 in the present study instead of an y other standard refrigeran ts. Section 6 is devoted to the rigorous presentation of the linear stability analysis that leads to a differential eigenv alue problem to b e treated via the Chebyshev- τ method. In the last section, numerical results are discussed and compared with existing exp erimental data confirming that the mo del captures the main features of the physical phenomenon and frames the dryout as a consequence of in terfacial instability . 2 Problem framing This section is intended to show a schematic representation of the ph ysical setup discussed in the introduction for mathematical mo delling purp oses. Hence, let us consider a horizontal pip e of diameter D where a tw o- phase annular flow develops. The flow inv olv es tw o immiscible fluids: F l , F v , where the former is in contact with the cylinder w alls and acts as refrigerant, the latter occupies the bulk and is generated b y phase c hange (ev ap oration) of F l . Moreo ver, the tw o fluids are separated by an interface, which is assumed to b e smo oth and whose profile changes as phase transition from F l to F v tak es place. T o introduce the framework, let us notice that the physical setup can b e view ed from tw o p ersp ectiv es: stream wise and spanwise (see Figure 1 ). In the presen t inv estigation, for mo delling purp oses, a streamwise p oin t of view is adopted. Moreov er, given the cylindrical shap e of the domain, it is assumed that the flow exhibits radial symmetry with resp ect to the cen tral axis of the cylinder and a longitudinal section of the pip e is considered. In this wa y , the domain reduces to a bidimensional horizontal lay er of depth D / 2 , where the upp er b oundary is identified by the cylinder axis, while the low er b oundary is the cylinder sidewall, sub ject to the presence of an external heat flux q . In this framework, as depicted in Figure 2 , F l o ccupies the region closer to the low er b oundary and the tw o fluids are separated b y a smo oth interface with decreasing profile. This profile iden tifies the phase transition that tak es place during the cooling process, as the thickness of the liquid film decreases. Let us now introduce a bidimensional Cartesian reference frame O z y with unit vectors { k , j } , resp ectiv ely , with j pointing upw ard. The generic fluid F α , ( α = l , v where l stands for “liquid" and v for “v ap our" ) is c haracterized b y dynamic viscosity µ α , density ρ α , stress tensor σ α and the unknown fields are velocity and pressure v α ( z , y , t ) = ( u α , w α ) and p α = p α ( z , y , t ) , resp ectively . Moreo ver, the interface separating the t wo sup erposed fluids is iden tified by the implicit e quation S ( x ) = 0 , where x = ( z , y ) . The implicit function theorem guarantees the existence of y = δ ( z ) such that S ( z , δ ( z )) = y − δ ( z ) = 0 (1) Let us remark that y = δ ( z ) is the equation that iden tifies the sloping interface separating the tw o fluids, see Figure 2. Figure 1: Sk etc h of distinct viewp oin ts for a cylindrical configuration. Stream wise p oint of view is adopted in this study . 4 A Preprint Figure 2: Sketc h of physical configuration, where the liquid film is attached to the pip e wall, while the bulk is o ccupied by the v ap our phase. The sloping interface b etw een the t wo fluids is identified b y δ ( z ) . External heat is provided to the system via the heat flux q . 3 In terfacial conditions Here, the equations that guaran tee conserv ation of mass, momen tum and energy across the in terface are deriv ed from first-principles theory . Indeed, the mathematical mo del arises from general conserv ation laws, whic h hold in each fluid phase and at the in terface separating the tw o phases. Let us define the quantit y ψ α = [ ψ α ] k , ( k = 1 , . . . , 5) within each fluid F α , ( α = l , v ) . The general form of conserv ation la ws for ψ α is: ∂ ψ α ∂ t + ∇ · f α = 0 , (2) where f α is the flux of ψ α . In contin uum mechanics, by emplo ying the eulerian point of view, Eq. ( 2 ) is sp ecialized b y defining, for i, j = 1 , 2 , 3 , [ f α ] ij = [ ψ α ] i [ v α ] j + [ G α ] ij , (3) and ψ α =    ρ α ρ α [ v α ] i 1 2 ρ α | v α | 2 + ρ α ε α    , G α =    0 j − [ σ α ] ij − [ σ α ] j k [ v α ] k + [ Q α ] j    (4) where ρ α is the fluid density , v α is the fluid velocity , σ α is the Cauch y stress tensor, ε α is the internal energy and Q α is the heat flux [19]. The resulting balance equations, expressing the conserv ation of mass, momen tum and energy , are:              ∂ ρ α ∂ t + ∇ · ( ρ α v α ) = 0 ∂ ( ρ α v α ) ∂ t + ∇ · ( ρ α v α ⊗ v α − σ α ) = 0 , ∂  ρ α ε α + 1 2 ρ ( v α ) 2  ∂ t + ∇ · h  ρ α ε α + 1 2 ρ ( v α ) 2  v α − σ α v α + Q α i = 0 (5) and hold within each fluid F α , ( α = l , v ) . No w, let us assume that densities ψ α are conserved at the in terface, namely Eqs. ( 2 )-( 4 ) hold within a domain Ω astride the interface. If S ( x , t ) = 0 is the equation that identifies the interface separating the tw o fluids, the following general conserv ation law at the interface holds [17]. 5 A Preprint f l ∇ S + ψ l ∂ S ∂ t = f v ∇ S + ψ v ∂ S ∂ t . (6) By substituting Eqs. ( 3 )-( 4 ) into Eq. ( 6 ), and by considering the first comp onent ( k = 1 ) of the resulting equation, the conserv ation of mass across the interface is obtained: ρ l  ∂ S ∂ t + v l · ∇ S  = ρ v  ∂ S ∂ t + v v · ∇ S  . (7) Index k = 2 , 3 , 4 leads to the conserv ation of momentum across the interface: ρ l v l  ∂ S ∂ t + v l · ∇ S  − σ l ∇ S = ρ v v v  ∂ S ∂ t + v v · ∇ S  − σ v ∇ S. (8) A t this stage, we should write the energy balance equation obtained at k = 5 . Ho w ever, the resulting problem w ould b e so complicated that a simplification step is needed. As anticipated in the introduction, we follow the idea of Hsieh [18], which consists of replacing the conserv ation of energy with a densit y balance equation where the heat transfer is mo delled as a pro duction term of densit y of mass. Indeed, according to Hsieh [18], transfer of mass across the interface represents a transformation from one phase to the other and therefore it is linked to latent heat asso ciated to phase change. Consequently , thermal effects are mo deled via mass transfer across the interface. Moreo ver, the p osition and modification of the interface corresp onds to phase c hange asso ciated with a certain latent heat. So, thermal effects are strictly connected to the instantaneous p osition of the interface. Hence, the interfacial condition for energy balance is replaced b y the following equation ρ l  ∂ S ∂ t + v l · ∇ S  = F ( y ) L , (9) where L is the latent heat required to conv ert liquid into v ap our at constant temp erature and F ( y ) represents the net heat flux across the interface when phase c hange takes place [18, 16]. Let us note that if L → ∞ , mass flux at the interface is absent as an infinite amount of heat would b e required for a phase change to happ en. Let us remark no w that, since the actual physical setup in volv es micro c hannels, Stevin’s contribution to pressure can b e dropp ed compared to surface tension effects, therefore we assume p − ρg y ≃ p , where p is the absolute pressure and g the gra vity . Therefore, the Cauch y stress tensor takes the following form: [ σ α ] hk = − p α δ hk + µ α  ∂ [ v α ] h ∂ x k + ∂ [ v α ] k ∂ x h  = − p α δ hk + [ ¯ σ α ] hk . (10) µ α b eing the dynamic viscosity of F α . Therefore, Eq. ( 8 ) b ecomes: ρ l v l  ∂ S ∂ t + v l · ∇ S  − ¯ σ l ∇ S − ( p v − p l ) ∇ S = ρ v v v  ∂ S ∂ t + v v · ∇ S  − ¯ σ v ∇ S. (11) Finally , let us assume that normal stresses due to pressure result in mo dification, near equilibrium, of the in terface p osition and shap e according to Y oung-Laplace equation: ∆ p = τ  1 R 1 + 1 R 2  (12) where ∆ p = p v − p l is the pressure difference at the in terface, τ is the surface tension and R 1 , R 2 are the principal radii of curv ature. In conclusion, the interfacial conditions are the follo wing:                      ρ l  ∂ S ∂ t + v l · ∇ S  = ρ v  ∂ S ∂ t + v v · ∇ S  , ρ l v l  ∂ S ∂ t + v l · ∇ S  − ¯ σ l ∇ S − τ  1 R 1 + 1 R 2  ∇ S = ρ v v v  ∂ S ∂ t + v v · ∇ S  − ¯ σ v ∇ S, ρ l  ∂ S ∂ t + v l · ∇ S  = F ( y ) L , u l = u v (13) where Eq. (13) 4 represen ts a regularity condition at the interface. 6 A Preprint No w, the go verning equations that hold within each fluid F α ( α = l , v ) are the incompressible Navier-Stok es equations:    ∇ · v α = 0 ∂ v α ∂ t + v α · ∇ v α = − 1 ρ α ∇ p α + µ α ρ α ∆ v α (14) that can b e derived from Eq. (5) 1 − 2 , by virtue of the stress tensor in Eq. (10). This system of gov erning equations is coupled with the follo wing b oundary conditions: v l = 0 on y = 0 v v = V k on y = D 2 (15) whic h mo dels no-slip condition at the upp er and low er b oundaries, as if the velocity of the fictitious upp er b oundary is V = V k . 4 Equations of motion and dimensionless form ulation The nondimensional formulation of Navier-Stok es equations and interfacial conditions is obtained by setting the diameter D of the tub e as characteristic length, the ratio G/ρ l of mass flux ov er liquid densit y as the c haracteristic velocity , while the characteristic time is ( D ρ l ) /G and the c haracteristic pressure is P = ( G ) 2 /ρ l . Hence, denoting by the asterisks the nondimensional quantities x = D x ∗ , v α = G ρ l ( v α ) ∗ , t = D ρ l G t ∗ , p α = P ( p α ) ∗ , (16) the Navier-Stok es equations for b oth fluids b ecome, omitting the asterisks for notational conv enience,      ∂ v l ∂ t + v l · ∇ v l = −∇ p l + 1 Re l ∆ v l ∂ v v ∂ t + v v · ∇ v v = − b ρ ∇ p v + 1 Re v b ρ ∆ v v (17) Here b ρ = ρ l ρ v and Re α = GD µ α (18) whic h are the Reynolds num b ers relative to the tw o fluids F α , α = l , v . They mo del the interpla y b etw een inertial and viscous forces, which act in different wa ys in the tw o fluids. Let us notice that the difference with a more standard definition of the Reynolds n um b er is consequence of the choice of the c haracteristic v elo cit y as G ρ l , working with the intensit y of the mass flux rather than the mean velocity of the fluid. Therefore, by defining ℓ α =  1 if α = l , b ρ if α = v, (19) the nondimensional Navier-Stok es equations can b e written in short form:    ∇ · v α = 0 ∂ v α ∂ t + v α · ∇ v α = − ℓ α ∇ p α + ℓ α Re α ∆ v α (20) with b oundary conditions Eq. (15) assigned on y = 0 , y = 1 2 . The interfacial conditions, instead, omitting the asterisks, b ecome                            ∂ S ∂ t + v l · ∇ S = ¯ ρ  ∂ S ∂ t + v v · ∇ S  , v l  ∂ S ∂ t + v l · ∇ S  − 1 Re l ¯ σ l ∇ S = ¯ ρ v v  ∂ S ∂ t + v v · ∇ S  − 1 Re v ¯ σ v ∇ S + 1 W e  1 R 1 + 1 R 2  ∇ S, ∂ S ∂ t + v l · ∇ S = F ( y ) L , u l = u v (21) 7 A Preprint where ¯ ρ = ρ v ρ l and the ratio µ l G ρ l D has b een chosen as characteristic quantit y for the stress tensors ¯ σ α , while W e = ( G ) 2 D τ ρ l (22) is the W eb er num ber. This nondimensional num b er plays a pivotal role in analysing interfacial problems. Indeed, it is defined as the ratio b etw een inertial forces and forces due to surface tension. The Na vier-Stok es equations (20) admit a steady solution denoted by m α b = ( v α b , p α b ) , where v α b ( z , y ) = ( u α b , w α b ) and the interface at the steady state is identified by Eq. ( 1 ). Let us remark that the velocity profile of liquid and v ap our phases satisfy Eq. (21) 4 . Therefore, at the interface y = δ ( z ) , v α b ( z , δ ( z )) = V b ( z , δ ( z )) where V b = ( U b , W b ) . Nevertheless, their deriv atives are different in general, therefore deriv atives of v α b at the interface are denoted by e.g. ( v α b ) y = ( V α b ) y . A t this stage of the study , we do not specify the basic v elo cit y profile of the t wo fluids. This choice is motiv ated by the desire of k eeping the analysis as general and broadly applicable as p ossible. Therefore, we pro ceed to p erform a linear stability analysis of this generic profile in order to determine a general form of the perturbation equations. How ever, numerical in v estigation is performed by selecting base velocity profiles v α b consisten t with annular flow. 4.1 The dryout instability factor In this paragraph, the main steps to determine the dimensionless dryout instabilit y factor defined in [13] are outlined. The reader is referred to this pap er for a detailed deriv ation. Given Eq. ( 1 ), it follows that ∂ S ∂ z = − d δ d z (23) F or the sak e of con venience, the function δ ( z ) is expressed as a comp osite function δ = δ ( α ( x ( z ))) so that the chain rule can b e applied: d δ d z = d δ d α d α d x d x d z (24) where α is the v oid fraction and x the v ap our quality . The deriv ativ es in Eq. (24) can b e calculated as in [13] and are: d δ d α = − D 4 √ α d α d x = α (1 − α ) x (1 − x ) d x d z = 4 q h gl GD (25) where h gl is the relative en thalp y b et ween liquid and gaseous phases. Substituting Eqs. (23)-(25) in the last term of Eq. (21) 2 pro jected on the z axis and b eing R 1 = D 2 √ α , R 2 = ∞ , the following expression is obtained: 1 W e  1 R 1 + 1 R 2  ∂ S ∂ z = 2 (1 − α ) x (1 − x ) Bo W e (26) The Boiling num b er Bo = q h gl G is defined as the ratio b etw een the rate of phase change ov er the total mass flow rate. while the ratio Bo W e = ρ l τ q G 3 h gl D exhibits an in verse prop ortionalit y to the mass flux. This mathematical form is consistent with the exp erimen tally–deriv ed δ + b eha viour [13]. In particular, it is shown that the dryout v ap our quality x dry is correlated with the right-hand-side of Eq. ( 26), defined as dryout instability factor I δ + := 2 (1 − α ) x (1 − x ) Bo W e (27) for the δ + regime. Let us remark that the importance of this num ber lies also in the geometrical configuration that has b een considered. In millic hannels, the surface tension plays a pivotal role, whereas if a larger diameter tub e were considered, the surface tension would b e small and I δ + w ould b e negligible to o. 8 A Preprint 5 CO 2 as the reference fluid for the δ + regime This section aims to give a detailed discussion on the reason why the ob ject of the present study fo cuses on t wo-phase CO 2 ann ular flo w in the δ + regime. As already underlined in the introduction, carb on dioxide uniquely accesses, under practically relev ant ev ap orator conditions, the thermo dynamic state required for the onset of the δ + regime. It is now essential to clarify why CO 2 can exhibit a different trend compared to other refrigerants. A comparison with conv entional refrigerants (e.g. R-12, R-134a, R-290) shows that, at the same ev ap orator conditions, unlike other parameters, liquid-to-v ap our density ratio and reduced pressure v ary significantly . In particular, CO 2 at typical refrigeration temp eratures op erates at m uch higher reduced pressure and exhibits a muc h smaller density ratio (i.e. more similar liquid and v ap our densities) than conv entional refrigeran ts. This is the key discriminating feature. T able 1 summarizes the comparison at identical standard conditions. The key p oint is that W e, Re, Bd and Bo remain essentially comparable in magnitude across fluids, while ρ l /ρ v and p r differ significantly . Fluid W e Re Bd Bo ρ l /ρ v p r ( − ) ( − ) ( − ) ( − ) ( − ) ( − ) CO 2 194,189 9,311 1.263 0.0154 16.6 0.359 R-12 71,969 4,157 1.013 0.0261 132.5 0.056 R-134a 79,032 3,697 0.965 0.0199 162.0 0.050 R-290 217,579 8,175 0.440 0.0106 84.3 0.056 T able 1: Dimensionless parameters for all refrigeran ts at standard conditions ( T sat = − 15 ◦ C, G = 1200 kg m − 2 s − 1 , D = 1 mm, q = 5 k W m − 2 ), data collected from [20]. CO 2 is a clear outlier in terms of liquid-to- v ap our density ratio and reduced pressure, indicating a near-critical op erating state with muc h more similar phase densities. This shows that CO 2 is the only fluid among standard refrigerants that realizes, at useful ev ap orator tem- p eratures, the lo w density-ratio condition required for δ + b eha viour. The small density ratio ac hieved by CO 2 reduces phase slip and w eakens in terfacial shear. This suppresses droplet entrainmen t, helps preserve a coheren t annular film, and allows increasing mass flux to enhance conv ectiv e ev aporation in the film without immediately disrupting it. Under these conditions, dryout can b e delay ed as G increases, which is precisely the signature of the δ + regime. The same conclusion can b e dra wn from T able 2 , where the p ercen tage deviation relative to CO 2 is reported. In fact, v ariations in W e, Re, Bd and Bo are significant but remain within factors of order unity . By con trast, the density ratio differs significantly , whic h yields different qualitative ann ular-flow dynamics. Fluid W e Re Bd Bo ρ l /ρ g p r (%) (%) (%) (%) (%) (%) CO 2 – – – – – – R-12 − 62.9 − 55.4 − 19.8 + 69.5 + 698 − 84.4 R-134a − 59.3 − 60.3 − 23.6 + 29.2 + 876 − 86.1 R-290 +12 . 0 − 12 . 2 − 65 . 2 − 31 . 2 +408 − 84 . 4 T able 2: Percen tage differences from CO 2 at the same standard conditions. The key discrepancy is the liquid-to-v ap our density ratio, which differs b y + 408% to + 876%, far exceeding the v ariation of the classical dimensionless groups. W e no w wan t to underline that the observed difference is not due to op erating conditions. In fact, let us fo cus on R-290 so as to repro duce the CO 2 reference case by adjusting the operating conditions (c hannel diameter, mass flux, and wall heat flux) so as to match the standard dimensionless groups as closely as p ossible. This is a stringent test of classical similitude arguments: if W e, Re, Bd, Bo num b ers were sufficient to repro duce similar v alues for the densit y ratio, then the R-290 case should repro duce the same δ + b eha viour as CO 2 . T able 3 shows that this is not the case. Ev en when W e, Re, Bd, and Bo are matched with excellen t accuracy , the densit y ratio remains muc h larger than in CO 2 . This remains true both at the same saturation temp erature and at elev ated saturation temp erature. 9 A Preprint Fluid D G q W e Re Bd Bo ρ l /ρ g (mm) (kg/m 2 s) (k W/m 2 ) ( − ) ( − ) ( − ) ( − ) ( − ) CO 2 ( − 15 ◦ C) 1.000 1200 5.0 194 9 1.263 0.0154 16.6 R-290 ( − 15 ◦ C) 1.694 871 5.3 194 10 1.263 0.0154 84.3 R-290 ( +15 ◦ C) 1.469 744 4.0 194 10 1.263 0.0154 32.1 T able 3: Optimization attempts to match CO 2 dimensionless parameters using R-290 by adjusting D , G , and q . Although W e, Re, Bd, and Bo can be matched closely , the liquid-to-v ap our density ratio remains substan tially higher than for CO 2 . This demonstrates tw o key p oin ts: 1. matc hing the classical dimensionless groups is p ossible by tuning op erating and geometric parame- ters; 2. matc hing the density ratio is not p ossible at practical ev aporator temp eratures. This shows that the density ratio is a key parameter for the δ + regime. In conclusion, to achiev e density ratios comparable to CO 2 , con ven tional refrigeran ts would need to operate m uc h closer to their critical p oin ts, i.e. at substantially higher saturation temp eratures. F or fluids such as R-290, this w ould require ev ap orator temp eratures far ab ov e the practical range used in refrigeration and air-conditioning applications. Suc h conditions are thermodynamically irrelev ant for co oling systems: an ev ap orator op erating at very high saturation temp erature provides little or no useful co oling relative to typical heat-sink temperatures. Therefore, even if a δ + b eha viour were theoretically p ossible for another refrigerant, it would o ccur outside the range where the ev ap orator is used as a refrigeration device. This is the k ey reason why CO 2 is the natural and practically relev ant fluid for the present study . Its lo w critical temp erature allows near-critical op eration (and hence low density ratio) within the temp erature range where ev ap orators are actually employ ed for co oling. The rationale for using CO 2 in this work can therefore b e summarized as follo ws: 1. The δ + regime is primarily controlled by thermo dynamic state, esp ecially phase-density similarity , rather than by fluid chemistry . 2. CO 2 uniquely provides low liquid-to-v ap our density ratio at practical ev aporator temp eratures b e- cause of its low critical temp erature. 3. Con ven tional refrigerants can match standard dimensionless groups but not the density ratio under useful refrigeration conditions. 4. The present in terfacial-stability-based model is physically consisten t precisely in the lo w-densit y- ratio annular-flo w conditions that CO 2 can realize. F or these reasons, CO 2 is not simply one refrigerant among many in the presen t study: it is the ph ysically and tec hnologically relev ant working fluid for inv estigating the δ + dry out regime in refrigeration-scale ev ap orators. 6 Linear analysis In this section, the linear stabilit y analysis is presen ted up on the introduction of small amplitude p ertur- bations on the unkno wn fields. F ollo wing [16], also the interface is perturb ed by a p erio dic function of infinitesimally small amplitude. This will lead to a couple of forth-order differential equations that require four conditions at the interface and four b oundary conditions. 6.1 Linearised gov erning equations for perturbations Let us introduce a p erturbation on the steady solution m α b , defined as follows v α ( z , y , t ) = v α b ( z , y ) + b v α ( z , y , t ) , p α ( z , y , t ) = p α b ( z , y ) + b p α ( z , y , t ) (28) where b v α ( z , y , t ) = ( b u α , b w α ) . Substituting Eq. (28) in to Eq. (20), the resulting system is autonomous, therefore p erturbation fields take the following form: b v α ( z , y , t ) = e v α ( y ) e ikz − nt , b p α ( z , y , t ) = e p α ( y ) e ikz − nt (29) 10 A Preprint where e v α ( y ) = ( e u α , e w α ) , k is the wa ven umber and n = n r + in i with n r b eing the growth rate. By virtue of Eqs. (28)-(29), the contin uit y equation Eq. (20) 1 , b ecomes: e u α ( y ) = − ˙ e w α ( y ) ik (30) whic h is a useful relation b etw een the first and second comp onen ts of velocity p erturbation and will allow us to reduce the num b er of unknown fields. Moreov er, by virtue of (28), the momentum equation Eq. (20) 2 , dropping nonlinear terms, b ecomes b u α t + u α b b u α z + w α b b u α y + b u α ( u α b ) z + b w α ( u α b ) y = − ℓ α b p α z + 1 Re α ℓ α ∆ b u α , b w α t + u α b b w α z + w α b b w α y + b u α ( w α b ) z + b w α ( w α b ) y = − ℓ α b p α y + 1 Re α ℓ α ∆ b w α . (31) No w, by deriving (31) 1 with respect to y and (31) 2 with respect to z and, by summing the resulting equations, b y virtue of the incompressibilit y constraint, we get b u α ty + u α b b u α z y + w α b b u α yy + b u α ( u α b ) z y + b w α ( u α b ) yy − 1 Re α ℓ α ( b u α z zy + b u α yy y ) = b w α tz + u α b b w α z z + w α b b w α yz + b u α ( w α b ) z z + b w α ( w α b ) yz − 1 Re α ℓ α ( b w α xxx + b w α yy x ) . (32) By virtue of Eqs. (29)-(30), Eq. (32) m ultiplied by ik , b ecomes: n  D 2 − k 2  e w α + 1 Re α ℓ α  D 2 − k 2  2 e w α − ik (( w α b ) yz − ( u α b ) yy ) e w α − w α b  D 2 − k 2  D e w α + (( w α b ) z z − ( u α b ) z y ) D e w α − ik u α b  D 2 − k 2  e w α = 0 . (33) with D = d dy . Hence, the p erturb ed linear gov erning equations are, in compact form: ( X Y )  e w l e w v  = n ( W W )  e w l e w v  (34) where X = 1 Re l  D 2 − k 2  2 − ik  ( w l b ) yz − ( u l b ) yy  − w l b  D 3 − k 2 D  +  ( w l b ) z z − ( u l b ) z y  D − ik u l b  D 2 − k 2  Y = 1 Re v b ρ  D 2 − k 2  2 − ik (( w v b ) yz − ( u v b ) yy ) − w v b  D 3 − k 2 D  + (( w v b ) z z − ( u v b ) z y ) D − ik u v b  D 2 − k 2  W = −  D 2 − k 2  (35) 6.2 Linearised interfacial conditions for perturbations Let us introduce a p erturbation ζ = ζ ( z , t ) on the interface so that Eq. ( 1 ) is replaced by: S ( x , t ) = y − δ ( z ) − ζ ( z , t ) , (36) and the p osition of the p erturb ed interface is iden tified by S ( x , t ) = 0 . Hence, the interfacial condition (21) 3 b ecomes: ρ l  ∂ S ∂ t + v l · ∇ S  = F ( δ + ζ ) L . (37) It is p ossible to recov er the following expression for F ( δ + ζ ) : F ( δ + ζ ) L = β ζ (38) where β = H L  1 D/ 2 − δ ( z ) + 1 δ ( z )  . By virtue of Eqs. (28)-(29) and by assuming ζ ( z , t ) = ξ e ikz − nt , Eqs.(37)- (38) b ecome, dropping nonlinear terms: − nξ − ik ξ U l b + ( ik ) − 1 δ ′ ˙ e w l + e w l = β ξ . (39) 11 A Preprint Setting γ = − ( ik ) − 1 , we get: ξ = − γ δ ′ ˙ e w l + e w l β + n + ik U l b . (40) Moreo ver, by virtue of (21) 1 , following the same reasoning, w e hav e ξ = ¯ ρ ( − γ δ ′ ˙ e w v + e w v ) β + ¯ ρ ( n + ik U v b ) (41) and the following equation is immediately obtained: − γ δ ′ ˙ e w l + e w l β + n + ik U b = ¯ ρ ( − γ δ ′ ˙ e w v + e w v ) β + ¯ ρ ( n + ik U b ) . (42) Let us define λ = β + ik U b and ¯ λ = ¯ ρ − 1 ( β + ¯ ρik U b ) , so that (42) can b e written as follo ws − ¯ λγ δ ′ ˙ e w l + ¯ λ e w l + λγ δ ′ ˙ e w v − λ e w v = n h − γ δ ′ ˙ e w v + e w v + γ δ ′ ˙ e w l − e w l i . (43) In matrix compact form, it b ecomes:  A 1 B 1   e w l e w v  = n  C 1 D 1   e w l e w v  (44) where A 1 = ¯ λ − ¯ λγ δ ′ D B 1 = − λ + λγ δ ′ D C 1 = − 1 + γ δ ′ D D 1 = 1 − γ δ ′ D (45) Therefore, w e provide the final p erturb ed interfacial condition relative to momentum, where nonlinear terms ha ve b een dropp ed. In particular, the equation for the first comp onent of (21) 2 is  A 2 B 2   e w l e w v  = n  C 2 D 2   e w l e w v  (46) where A 2 = λ  [ C 1 + B ] + [ E 1 − γ δ ′ B ] D + F 1 D 2  B 2 = − λ  C 2 + E 2 D + F 2 D 2  C 2 = − C 1 − A + [ − E 1 + γ δ ′ A ] D − F 1 D 2 D 2 = C 2 + E 2 D + F 2 D 2 (47) and A = U b ( ¯ ρ − 1) , B = − ik ( U b ) 2 + ik ¯ ρ ( U b ) 2 + 1 Re l [2( U l b ) z ] ik − 1 Re v [2( U v b ) z ] ik − 1 W e k 2 δ ′ + 1 W e  1 R 1 + δ ′′  ik , C 1 = U b − 1 Re l ik , C 2 = ¯ ρU b − 1 Re v ik , E 1 = − 2 γ U b δ ′ + W b γ + 2 Re l ik γ δ ′ , E 2 = − 2 γ U b ¯ ρδ ′ + W b γ ¯ ρ + 2 Re v ik γ δ ′ , F 1 = − 1 Re l γ , F 2 = − 1 Re v γ . (48) While, the equation for the second comp onen t of (21) 2 is  A 3 B 3   e w l e w v  = n  C 3 D 3   e w l e w v  (49) 12 A Preprint where A 3 , B 3 , C 3 , D 3 are defined as in Eq. (47) but Eq. (48) is replaced by A = W b ( ¯ ρ − 1) , B = − ik W b U b + 1 Re l [( U l b ) y + ( W l b ) z ] ik + 1 W e k 2 + ik ¯ ρU b W b − 1 Re v [( U v b ) y + ( W v b ) z ] ik , C 1 = 2 W b − δ ′ U b + 1 Re l δ ′ ik , C 2 = 2 ¯ ρW b − ¯ ρU b δ ′ + 1 Re v ik δ ′ , E 1 = − γ W b δ ′ − 2 Re l , E 2 = − γ W b ¯ ρδ ′ − 2 Re v , F 1 = 1 Re l γ δ ′ , F 2 = 1 Re v γ δ ′ . (50) Let us remark that the last interfacial condition Eq. (21) 4 is simplified by virtue of Eq. (30) and b ecomes: ( D − D )  e w l e w v  = n ( 0 0 )  e w l e w v  (51) 7 Numerical metho d and results This section is in tended to comment on the numerical results obtained by solving the differen tial eigenv alue problem made of Eqs. (34), (44), (46), (49), (51) together with b oundary conditions. In a first instance, the base flo w is assigned to liquid and v ap our phases. Then, the numerical metho d is illustrated and v alidated against the Orr-Sommerfeld problem studied in [21]. Afterwards, the main steps of the numerical pro cedure are outlined and the comparison b etw een numerical and exp erimental data is discussed. 7.1 The base flow Let us define explicitly the base v elo cit y in the liquid and v ap our phases. The following hypotheses ha ve b een enforced: 1. The horizontal velocity in the v ap our core has a flat profile along the radius 2. The horizontal velocity in the liquid film is linear: 0 at the wall and equal to the v ap our velocity at the interface These profiles are consistent with exp erimental evidence. In particular, the Reynolds num b er in the v apour core is sufficien tly high that a uniform velocity profile is a reasonable approximation. Whereas, the liquid film is very thin, therefore, as a first approximation, the velocity can b e assumed to v ary linearly ov er the film thickness. Under these hypotheses, the base velocities v α b = ( u α b , w α b ) for α = l, v can be expressed in terms of the liquid film thickness δ and the v ap our quality x . The horizontal comp onen t of the v ap our phase velocity is obtained by definition: u v b = G ρ v x α (52) whic h, in dimensionless form, b ecomes: u v ∗ b = u v b ρ l G = ρ l ρ v x α (53) The horizontal velocity profile in the liquid film is linear, as mentioned in the previous hypothesis n. 2, u l ∗ b ( r ) = u v ∗ b D 2 − r δ , r ∈  D 2 − δ, D 2  (54) b eing δ = D δ ∗ . F rom the horizontal profile, by virtue of the contin uity equation and the condition stating that the velocity has to b e parallel to the interface: W b ( z , δ ( z )) U b ( z , δ ( z )) = d δ d z (55) the expressions for the v ertical comp onents can be determined. Once the v elocities profiles hav e b een iden ti- fied, their deriv atives appearing in Eqs. (34), (46) and (49) can b e easily computed, where the deriv ativ es of Γ α and K α can b e calculated explicitly from their expressions or numerically with finite difference method. 13 A Preprint -1 -0.8 -0.6 -0.4 -0.2 0 0 0. 2 0. 4 0. 6 0. 8 1 c i [-] c r [-] Dongar ra data In-house d at a Figure 3: Comparison b et ween the sp ectrum obtained in Figure 7 in [21] and the n umerical results obtained with the in-house co de. 7.2 V alidation A common tool to solve a differential eigenv alue problem originating from h ydrodynamic stability analysis is the Chebyshev- τ metho d [22]. Details ab out the implemen tation of the metho d and the adv an tage of using Cheb yshev p olynomials are o verlooked and the reader is referred to [23, 24]. Let us just summarise the main idea for implementation in case of tw o sup erp osed fluids o ccup ying a horizon tal region. This problem is addressed in [21], whose results are used to v alidate our in-house co de. F or tw o superp osed fluids, given that Chebyshev p olynomials are naturally defined in [ − 1 , 1] , the n umerical metho d is applied after a proper change of coordinates. Let us consider a tw o-dimensional domain where a fluid, marked with 2 , is o ccupying the region − 1 ≤ y ≤ 0 underlying a different fluid denoted by 1 , b ounded in 0 ≤ y ≤ m . In each region, the Orr-Sommerfeld equations hold and b oundary conditions and interfacial conditions are app ended to them. The differential eigenv alue problem to b e solved is given by Eqs. 6.4-6.7 in [21]. In order to employ the Cheb yshev- τ metho d in such a problem, the following change of co ordinates is p erformed z 1 = − 2 y − 1 z 2 = 2 m y − 1 (56) This c hange the co ordinates maps, resp ectively , y ∈ ( − 1 , 0) and y ∈ (0 , m ) in to ( − 1 , 1) , with the interface placed at z 1 = z 2 = − 1 for b oth fluids. The deriv ativ es hav e to b e rewritten in the following wa y d d y = − 2 d d z 1 y ∈ [ − 1 , 0] d d y = 2 m d d z 2 y ∈ [0 , m ] (57) Hence, the differential eigenv alue problem can be written as A ϕ = cB ϕ (58) b eing ϕ the eigenfunction, with A =        A 1 0 0 A 2 B C inter,c B C inter,c B C inter, 0 B C inter, 0 B C 1 0 0 B C 2        B =        B 1 0 0 B 2 B C inter,c B C inter,c 0 0 0 0 0 0        (59) where A i and B i are the op erator of the left and right side, resp ectively , of the Orr-Sommerfeld equation for the fluid i . B C inter,c are the rows for the b oundary conditions at the interface containing the eigenv alue c , 14 A Preprint -100 -80 -60 -40 -20 0 20 40 60 80 0 20 40 60 80 100 120 140 n r n i n r ,i [-] N [-] Figure 4: V alue of the real and imaginary part of the leading eigenv alue as a function of N , with T sat = − 15 ◦ C, q = 30 k W/m 2 , G = 1200 k g /m 2 , x = 0 . 7 and k = 0 . 1 . B C inter, 0 are the rows for the b oundary conditions at the interface not containing the eigenv alue and B C i are the b oundary conditions inv olving only one of the tw o fluids. The b oundary condition are used to make the matrix B no longer singular in order to av oid the o ccurrence of spurious eigenv alues [25, 26]. The in-house co de has b een tested and v alidated using the data presented in [21] as a b enc hmark and such a comparison is av ailable in Figure 3 that shows coherence b etw een the sp ectrum computed by the in-house co de and the b enchmark. 7.3 Results and discussion The n umerical method discussed in the previous subsection is used to solve the system made by the Eq. (34) and the b oundary conditions (44)-(46)-(49), once a base flow profile has b een chosen according to section 7.1. Unlik e the problem studied by [21], the origin of the cartesian frame is set on the pip e wall (see Figure 2 ). Therefore, in order to recov er the setting of [21], follo wing transformation has to b e applied ˜ y = y δ − 1 (60) In the new reference frame, the origin is placed on the liquid-v ap our interface and the pip e wall is iden tified b y ˜ y = − 1 . No w, the transformation prop osed in [21] can be p erformed for each fluid-o ccupied region. Hence, for the liquid-phase region ˜ y = − y l + 1 2 , y l ∈ [ − 1 , 1] d d ˜ y = − 2 δ d d y l (61) while for the v ap our-phase region ˜ y = m ( y v + 1) 2 , y v ∈ [ − 1 , 1] d d ˜ y = 2 mδ d d y v (62) where m = 1 2 δ − 1 . No w, let us enumerate the steps follow ed to solve the stability problem, once T sat , G, q and k are fixed: 1. The base flow and all its physical quan tities (co efficien ts of both phases v elo cities, void fraction, Boiling num ber, etc.) are computed in each p oin t along the pip e. 2. F or a prescrib ed v alue of the v apour quality x , and of the corresp onding δ ( z ) , the aforementioned c hange of reference is p erformed. 15 A Preprint -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 -1500 -1000 -500 0 500 n i [-] n r [-] Figure 5: Sp ectrum obtained with Chebyshev- τ method with N = 20 p olynomials, for T sat = − 15 ◦ C, q = 30 k W/m 2 , G = 1200 k g /m 2 , x = 0 . 7 and k = 0 . 1 . 3. The differen tial eigenv alue problem is written in a matrix form A x = nB x b y virtue of the Cheb yshev- τ metho d using N polynomials. 4. The algebraic eigen v alue problem is solv ed and the sign of the real part of the most unstable eigen- v alue n r is studied: if n r > 0 the liquid-v ap our interface is disrupted, triggering dryout. 5. The pro cedure contin ues from p oint 2 for a different v alue of v ap our qualit y . The generalized eigenv alue problem A x = nB x has the following form, using the notation shown in Section 6 :               X 0 0 Y A 1 B 1 A 2 B 2 A 3 B 3 D − D 0 1 0 D 1 0 D 0                e w l e w v  = n               W 0 0 W C 1 D 1 C 2 D 2 C 3 D 3 0 0 0 0 0 0 0 0 0 0                e w l e w v  (63) where matrices A and B are to be considered as algebraic matrices and not differential op erators. Being N the n umber of Chebyshev p olynomials emplo yed for discretization, dimension of A and B in Eq. (63) is (2 N + 8) × (2 N + 8) . Matrix B is made non-singular by following a standard pro cedure via the last five ro ws, reducing matrices dimension to (2 N + 3) × (2 N + 3) . The in-house co de has been tested by v arying the num b er of Chebyshev p olynomials N . Figure 4 shows one of the test cases. As usually happ ens for the Cheb yshev- τ metho d, the eigenv alue conv erges when N is appro ximately in the range (12 , 40) . Indeed, results from the numerical metho d oscillate when N is large. This numerical phenomenon is due to the truncation error done when appro ximating the eigenfunctions with a fixed n umber of Cheb yshev polynomials. In particular, the dimension of the matrix operators increases as N 2 and their entries b ecome notably large exceeding the precision of the computational unit [21]. In the presen t study , N = 20 has b een chosen. F or a prescribed set of parameters, running the co de to solv e the eigenv alue problem will pro vide the sp ectrum as in Figure 5 . F or this case, the leading eigenv alue has a p ositive real part, showing that the interface is unstable. Once the procedure has b een iterated o v er a range of x and the real part of the leading eigen v alue at eac h step is computed and stored, results of the procedure can b e display ed lik e in Figure 6 . The most unstable eigenv alue has a negative real part until a certain v alue of v ap our quality is reached. This v alue marks the transition b et ween a stable and an unstable interface and it is named the dryout vap our quality x dry . Therefore, when the v ap our quality x exceeds x dry , the liquid-v ap our in terface is unstable and dryout inception takes place. 16 A Preprint Figure 6: Plot of the real part of the most unstable eigenv alue for T sat = − 15 ◦ C, G = 2000 k g / ( m 2 s ) and q = 30 k W/m 2 . The green shaded area represents the condition of stable v ap our quality where the dryout do es not take place. 0 500 1000 1500 2000 2500 0 0.2 0.4 0.6 0.8 1 Experimental data Mathematical model k = 0.1 G [kg/(m 2 s)] x dry [-] Figure 7: Comparison b etw een exp erimen tal data and numerical results for T sat = − 15 ◦ C and q = 30 k W/m 2 . In Figure 7 , the b eha viour of the dryout v ap our quality x dry is rep orted as a function of the mass flux G . The authors in [13] hav e shown that in their exp eriments, increasing the inlet mass flux within the millichannel has a delaying effect of the o ccurrence of dry out. F rom our persp ective this should result in a stabilising effect of G on the liquid-v apour interface. The numerical pro cedure explained earlier in the pap er has b een rep eated for a set of prescrib ed v alues of the mass flux G . Our numerical results show go od agreement with exp erimen tal data, as shown in Figure 7 . Hence, the trend of δ + regime has b een captured by the present mo del, with the only exception of the low est v alue of the mass flux where exp erimen tal errors are larger and non-linear phenomena are triggered, making the linear analysis presented so far no longer suitable. Let us also remark that the w av en umber k has b een chosen in such a wa y that the resulting x dry w as as close as p ossible to the exp erimental v alue for G = 2000 k g / ( m 2 s ) . F urther, numerical results hav e b een compared with the whole exp erimental dataset obtained during the campaign at CERN facilities for [13]. The dryout instabilit y factor embo dies the comp ound effect of all the relev ant parameters on the onset of dry out and it can b e computed from the set of parameters a v ailable to run numerical inv estigations. As a result, once the dryout v ap our quality is obtained from the linear stabilit y analysis, it can b e plotted as a function of I δ + . Figure 8 summarises the results, sho wing again v ery go od agreement of numerical results with exp erimen tal data in the tested range of parameters. Our mo del has b een tested also for a different configuration in the δ + regime. This configuration is the ob ject of the 17 A Preprint 1.0e-10 1.0e-09 1.0e-08 1.0e-07 1.0e-06 1.0e-05 1.0e-04 1.0e-03 1.0e-02 0 0.2 0.4 0.6 0.8 1 Model D = 0.5 mm Experiments D = 0.5 mm Model D = 1mm Experiments D = 1 mm I δ [-] x dry [-] Figure 8: Comparison betw een results obtained via numerical inv estigations on the theoretical mo del and exp erimen tal data from tw o indep endent campaigns [13, 9 ]. Exp erimental data defining I δ + are plotted against the dryout v ap our quality x dry , determined b oth exp erimentally and numerically . researc h pap er [ 9 ] and inv olves a millichannel of diameter D = 0 . 5 mm. The agreement b et ween n umerical results and exp erimental data strengthen the v alidity of the approach taken in the present research pap er. In conclusion, the numerical predictions, along with the theoretical framework dev elop ed here, support the prop osed ph ysical mechanism whereby the dryout inception in millichannels for CO 2 in the δ + regime is a consequence of the disruption of the interface b etw een liquid and v ap our phases, and is therefore triggered b y interfacial instability . References [1] M. M. W eislogel and J. B. McQuillen, “Hydro dynamic dryout in tw o-phase flo ws: Observ ations of low b ond num ber systems,” in AIP Confer enc e Pr o c e e dings , vol. 420, pp. 413–421, American Institute of Ph ysics, 1998. [2] T. L. Bergman, A. S. Lavine, F. P . Incrop era, and D. P . DeWitt, Intr o duction to he at tr ansfer . John Wiley & Sons, 2011. [3] Collier and Thome, Conve ctive Boiling and Condensation . Oxford Universit y Press, 1994. [4] G. Hewitt, A nnular two-phase flow . Elsevier, 2013. [5] P . Gullo, K. T samos, A. Hafner, Y. Ge, and S. T assou, “State-of-the-art technologies for transcriti- cal R744 refrigeration systems–a theoretical assessment of energy adv antages for Europ ean fo od retail industry ,” Ener gy Pr o c e dia , v ol. 123, pp. 46–53, 2017. [6] B. V erlaart, “Co2 co oling is getting hot in high-energy physics,” CERN Courier , 2012. [7] R. Rev ellin, P . Haberschill, J. Bonjour, and J. Thome, “Conditions of liquid film dry out during saturated flo w b oiling in micro c hannels,” Chemic al engine ering scienc e , vol. 63, no. 24, pp. 5795–5801, 2008. [8] R. Y un and Y. Kim, “Critical quality prediction for saturated flow b oiling of CO2 in horizontal small diameter tubes,” International Journal of He at and Mass T r ansfer , v ol. 46, no. 14, pp. 2527–2535, 2003. [9] M. Ducoulom bier, S. Colasson, J. Bonjour, and P . Hab ersc hill, “Carb on dioxide flow b oiling in a single micro c hannel–P art I I: Heat transfer,” Exp erimental Thermal and Fluid Scienc e , vol. 35, no. 4, pp. 597– 611, 2011. [10] L. W o jtan, T. Ursenbac her, and J. Thome, “Inv estigation of flow b oiling in horizon tal tub es: P art IA new diabatic tw o-phase flow pattern map,” International Journal of He at and Mass T r ansfer , vol. 48, no. 14, pp. 2955–2969, 2005. [11] L. Cheng, G. Ribatski, J. Moreno Quibén, and J. Thome, “New prediction metho ds for CO2 ev ap oration inside tub es: Part I A t wo-phase flow pattern map and a flo w pattern based phenomenological mo del for tw o-phase flow frictional pressure drops,” International Journal of He at and Mass T r ansfer , vol. 51, no. 1, pp. 111–124, 2008. 18 A Preprint [12] D. Hellenschmidt and P . P etagna, “Effects of saturation temperature on the boiling prop erties of carb on dio xide in small diameter pip es at low v ap our qualit y: Heat transfer co efficient,” International Journal of He at and Mass T r ansfer , vol. 172, p. 121094, 2021. [13] G. Can tini, D. Hellenschmidt, C. P edano, P . Petagna, C. Sangan, and M. Carnev ale, “Inception of ev ap orative dryout for CO2 in milliscale pip e flo ws,” International Journal of He at and Mass T r ansfer , v ol. 251, p. 127299, 2025. [14] S. Chandrasekhar, Hydr o dynamic and hydr omagnetic stability . Do ver, 1981. [15] P . Drazin and W. Reid, Hydr o dynamic stability . Cam bridge universit y press, 2004. [16] K. Adham-Khodaparast, M. Kaw a ji, and B. Antar, “The Rayleigh–T aylor and Kelvin–Helmholtz sta- bilit y of a viscous liquid–v ap or interface with heat and mass transfer,” Physics of Fluids , vol. 7, no. 2, pp. 359–364, 1995. [17] D. Y. Hsieh, “Effects of heat and mass transfer on Rayleigh-Ta ylor instability ,” 1972. [18] D. Y. Hsieh, “Interfacial stability with mass and heat transfer,” The Physics of Fluids , vol. 21, no. 5, pp. 745–748, 1978. [19] T. Ruggeri and M. Sugiyama, R ational extende d thermo dynamics b eyond the monatomic gas . Springer, 2015. [20] National Institute of Standards and T echnology (NIST), “Nist chemistry w ebb o ok, srd 69,” 2025. [21] J. Dongarra, B. Straughan, and D. W alker, “Chebyshev tau-QZ algorithm methods for calculating sp ectra of hydrodynamic stability problems,” Appl. Numer. Math. , vol. 22, pp. 399–434, 1996. [22] C. Canuto, A. Quarteroni, M. Y. Hussaini, and T. A. Zang Jr, Sp e ctr al metho ds: evolution to c omplex ge ometries and applic ations to fluid dynamics . Springer. [23] D. Bourne, “Hydro dynamic stabilit y , the Cheb yshev tau method and spurious eigenv alues,” Continuum Me chanics and Thermo dynamics , vol. 15, no. 6, pp. 571–579, 2003. [24] G. Arnone, J. A. Gianfrani, and G. Massa, “Cheb yshev- τ metho d for certain generalized eigen v alue problems occurring in hydrodynamics: a concise surv ey ,” The Eur op e an Physic al Journal Plus , vol. 138, no. 3, p. 281, 2023. [25] G. McF adden, B. Murray , and R. Boisvert, “Elimination of spurious eigenv alues in the Chebyshev tau sp ectral metho d,” Journal of Computational Physics , vol. 91, no. 1, pp. 228–239, 1990. [26] G. Arnone, G. Can tini, F. Cap one, and M. Carnev ale, “The onset of p enetrativ e conv ection in an inclined p orous lay er,” International Journal of He at and Mass T r ansfer , vol. 216, p. 124532, 2023. 19