A heat transfer model for liquid film boiling on micro-structured surfaces

ABSTRACT High heat transfer coefficient (HTC) and critical heat flux (CHF) are achieved in liquid film boiling by coupling vibrant vapor bubbles with a capillary liquid film, which has thus received increased interest for thermal management of high-power electronics. Although some experimental progress has been made, a high-fidelity heat transfer model for liquid film boiling is lacking. This work develops a thermal-hydrodynamic model by considering both evaporation atop the wick and nucleate boiling inside the wick to simultaneously predict the HTC and CHF. Nucleate boiling is modeled with microlayer evaporation theory, where a unified scaling factor is defined to characterize the change of microlayer area with heat flux. The scaling factor η is found to be independent of wicking structure and can be determined from a few measurements. This makes our model universal to predict the liquid film boiling heat transfer for various micro-structured surfaces including micropillar, micropowder, and micromesh. This work not only sheds light on understanding fundamental mechanisms of phase-change heat transfer, but also provides a tool for designing micro-structured surfaces in thermal management.


INTRODUCTION
Thermal management is becoming increasingly important for advanced electronic and energy systems, such as radars, microprocessors, power inverters, and space systems, where a large amount of heat needs to be dissipated from a limited space and area, and most importantly with a small temperature difference [1 -5 ].Among various thermal management techniques, liquid-vapor phase-change-based cooling strategies, such as capi l lary evaporation [6 -9 ] and immersion cooling with pool boiling [10 -13 ], have attracted great attention due to their excellent heat transfer performance.However, conventional cooling methods utilizing liquid-vapor phase-change processes encounter challenges in simultaneously improving the heat transfer coefficient (HTC) and critical heat flux (CHF) on the same wicking structure [4 ,9 ,14 ,15 ].Recently, capi l larydriven liquid film boiling where vapor bubbles are generated within a wicked liquid film has shown some promise in simultaneously enhancing HTC and CHF [16 ,17 ].Despite the experimental progress on enhancing the heat transfer performance of liquid film boiling on novel micro-structured surfaces [18 -22 ], there sti l l lacks a high-fidelity model that can capture the detailed physical phenomena and predict the CHF and HTC, and guide the design of wicking surfaces.During liquid film boiling, the heat first passes from the heated surface to the solid-fluid matrix of the wick, and then dissipates either by nucleate boiling through vapor bubbles inside the wick or by evaporation through the thin-film region atop the wick [16 ].As the heat flux increases, the nucleate boiling becomes dominant, leading to corresponding changes in the heat transfer ratio between nucleate boiling and evaporation atop the wick, as well as the dynamic characteristics of the vapor bubble and liquid film meniscus.Characteristics of such coupling phenomena are further complicated by the variety of wicking structures, which vary both the two-phase capi l lary delivery and the thermal transport in the solid-fluid matrix [16 ,23 ,24 ].It is thus challenging but highly desirable to develop a high-fidelity heat transfer model that can account for the thermal-hydraulics and interfacial processes on micro-structured wicking surfaces.
Recently, empirical correlations for CHF of liquid film boiling on the surfaces of micropi l lar arrays [25 ] or copper inverse opals [26 ] have been proposed.It is not known whether these specific fitting formulas based on their own experiments can be used for other structures.Determining the HTC of liquid film boiling under varying heat flux can be much more challenging, as it dynamically varies with the expanded liquid-vapor interfacial area of vapor bubbles coupled with the evaporative meniscus [16 ,18 ,20 ].Recently, a model was developed for liquid film boi ling [27 ] by assuming a constant liquid-vapor interfacial area, i.e. not vary ing w ith heat flux and modeling nucleate boiling inside the wick through pore-scale evaporation.The evaporation atop the wick is neglected.Under these assumptions, the model has good agreement with experimental data after an empirical parameter was adjusted based on the specific wicking structure [27 ].
In this work, we develop a model to predict both HTC and CHF of liquid film boiling for various micro-structured wicking surfaces.Both evaporation atop the wick and nucleate boiling inside the wick are taken into consideration.Evaporation atop the wick is determined by analyzing the liquid meniscus curvature variation and the thermal resistance network of the thin-film region close to the tri-phase line [6 ,28 ].Nucleate boiling inside the wick is analyzed at pore-scale through the microlayer evaporation borrowed from pool boiling [29 -32 ], with a scaling factor related to heat flux being introduced.Liquid film boiling experiments on micromesh surfaces were conducted along w ith w ick wettability, wicking capability, and structural characterizations, to determine the scaling factor η. According to our experiments, the scaling factor η is found to be independent of structural parameters, which makes our model universal to predict the liquid film boi ling heat transfer for various micro-structured surfaces.Our model predictions are in good agreement with the reported experimental data on various types of microstructured surfaces, including micromesh, micropowder, and micropi l lar.

MODELING AND ANALYSIS
This work aims to develop a unified model for the prediction of HTC and CHF of liquid film boiling on various micro-structured surfaces with different liquid supply methods.Typical microstructures include micromesh, micropowder, and micropi l lar (Fig. 1 A-C), while the liquid can be delivered with one-side [16 ], two-side [33 ], and all-around [21 ] supply methods (Fig. 1 D-F).Liquid film boiling with one-side liquid supply on a generalized microstructured surface is i l lustrated in Fig. 1 G.A set of key structural, thermophysical, and wicking parameters including porosity ε w , permeability K w , effective pore radius r eff , effective thermal conductivity k eff , and characteristic radius of wick skeleton r c relating to three typical microstructures, including micromesh, micropowder, and micropi l lar, are listed in Note S1.Similarly, different liquid supply methods alter liquid deliver y boundar y conditions as seen in Note S2.
Figure 1 G i l lustrates transport in the x to z direction, i.e. heat is transferred mainly in the thickness direction ( z -direction), while the liquid flows along the wicking direction ( x -direction).On the one hand, heat is absorbed from the heated surface to the liquid film, and transfers inside the liquid film through heat conduction or activates nucleate boiling.On the other hand, the flow of thin liquid film inside the wicking structures depends on both the wicking resistance usually characterized by liquid permeability and capi l lary pressure, in addition to boiling and vaporization.Here, liquid permeability depends on vapor bubble distribution during the boiling process.
The following assumptions are made to model the physical processes: (1) nucleate boiling occurs uniformly under a constant heat flux across the surface since most wicks are planar-isotropic, perpendicular to the direction of thickness.
(2) Heat conduction through the wicked liquid film is one-dimensional along the z -direction, since wick thickness along the z -direction is generally much smaller than its length along the x -direction.
(3) The temperature of the vapor within the micropore is uniform along the z -direction and equal to the saturated vapor temperature.This is reasonable since the liquid film thickness is typically less than hundreds of micrometers, leading to negligible pressure variation that drives the temperature change [34 ]. (4) Liquid can be well wicked to maintain a liquid meniscus on top of the wick; the equilibrium liquid meniscus is used to calculate the evaporation atop the wick.(5) For nucleate boiling inside the wick, the microlayer evaporation model can be adopted as has been widely used in pool boiling [29 ,32 ,35 ]. (6) The evaporation atop the wick and nucleate boiling inside the wick can be characterized as heat transfer coefficients h e and h bv , respectively.
In the following, we first find the temperature distribution across the wicking structure as a function of the heat transfer coefficient based on the heat balance equation.The heat transfer coefficient h e of Figure 1.Capillary-driven liquid film boiling on the micro-structured surface.Typical wicking structures include (A) micromesh [16 ], (B) micropowder [21 ], and (C) micropillar [25 ] studied in this work.Common liquid supply methods for liquid film boiling include (D) one-side [16 ], (E) two-side [33 ], and (F) all around [21 ]. (G) Schematic of liquid film boiling heat transfer on a micro-structured surface, where nucleate boiling occurs inside the wicking structure and the evaporation dissipates heat atop the wick.(H) Control volume for heat transfer modeling in the thickness direction ( z -direction) with a unit width.(I) Control volume for liquid transport along the wicking direction ( x -direction).Panel (A) is reproduced with permission from ref. [16 ], Copyright 2018 Elsevier.Panel (B) is reproduced with permission from ref. [21 ], Copyright 2019 Elsevier.Panel (C) is reproduced with permission from ref. [25 ], Copyright 2014 Elsevier.
evaporation atop the wick is derived by analyzing the thermal resistance network relating to the curvature of the liquid meniscus.The heat transfer coefficient h bv of nucleate boiling inside the wick is determined by a pore-scale analysis using the microlayer evaporation theory.The transport of the liquid film inside the wicking structure is modelled using the Brinkman equation.The CHF of boiling can be due to many mechanisms; we consider here the CHF of liquid film boiling to be determined by liquid wicking failure assuming the liquid on the heated surface is completely supplied by surface wicking [26 ].By solving the coupled momentum and energy equations, we obtain theoretical expressions for the HTC h t and the CHF q CHF of liquid film boi ling on various micro-structured surfaces.

Temperature distribution across the wicking liquid film
Considering a control volume with a unit width inside the wick (Fig. 1 H), the heat balance equation can be expressed as q c ( z ) = q c ( z + d z ) + d q b , where q c = q c (d x 1) is the heat conduction rate in the z -direction, q c is the conduction heat flux, and q b is the nucleate boiling heat transfer rate.Given that the conduction heat flux at z + d z can be expressed as q c (z + d z ) = q c (z ) + dq c (z ) / d z , the heat balance equation can then be expressed as: We note that an equivalent homogeneous wick is adopted in Fig. 1 H for generalization, and this control volume is arbitrarily independent of wicking structure.The conduction heat flux in the zdirection is given by: where T w is the temperature of the wicking structure.Similar to Ref. [27 ], we assume that heat is dissipated by nucleate boiling with a uniformly volumetric heat transfer coefficient h bv by averaging the size of bubbles along the wick.The heat transfer rate by nucleate boiling can then be obtained as: where T v is the vapor temperature.Substituting Eqs. ( 2) and ( 3 ) into Eq.( 1), the governing equation of T w across the wick can be expressed as: Equation ( 4) is essentially a fin equation, which is a second-order ordinary differential equation (ODE).Two boundary conditions at the bottom ( z = 0) and top ( z = δ l ) of the liquid film must be imposed for solving Eq. ( 4).The heat flux at the heated surface ( z = 0) is equal to the total heat flux, and this boundary condition can be expressed as [27 ]: where q t is the total heat flux.Considering evaporation happens at the top surface of the wick ( z = δ l ) [16 ], the boundary condition at the top surface should be written as: where q e and h e are the heat flux and heat transfer coefficient of the evaporation atop the wick, respectively.δ l is the thickness of the wicked liquid film.
We note that the evaporation atop the wick was neglected in Ref. [27 ] by writing the boundary condition k e f f (dT w / d z ) | z =δ l = 0 .Integrating Eq. ( 4) with boundary conditions of Eqs. ( 5) and ( 6), the temperature distribution of the wicking structure can be obtained as: where m = h bv /k e f f .In the forthcoming sections, δ l and h e are determined through analysis of the transport and evaporation of the wicking liquid film, while h bv is determined by pore-scale analysis of microlayer evaporation.T w ( z ) can be determined according to Eq. ( 7) with known δ l , h e , and h bv , and thus the HTC of liquid film boiling can then be obtained by writing

Transport and evaporation of the wicking liquid film
The liquid film transport inside the wicking structure can be described by the Brinkman equation as [36 ]: where u l is the liquid flow velocity, ε l is the effective porosity, μ l is the liquid dynamic viscosity, K w is wick permeability, K rl is the relative liquid permeability [37 ,38 ], P l is the local liquid pressure, ρ l is the liquid density, and g is the gravitational acceleration.The boundary conditions for Eq. ( 8) can be obtained with the non-slip velocity at the substrate and shear-free stress atop the wick as: Equation ( 8) is a second-order linear ODE, which can be solved with the boundary conditions of Eq. ( 9) to obtain liquid wicking velocity u l ( x , z ).Then, the average liquid wicking velocity for arbitrary cross-section at x can be expressed as ūl (x ) = δ l 0 u l (x, z )d z/δ l [6 ], and thus the liquid pressure gradient d P l / d x driving the liquid film can be obtained as: where λ = √ ε l /K rl K w .The relative liquid permeability K rl in Eq. ( 10) is computed using the empirical relation as K rl = s 3 [38 ], where s is the liquid saturation.This liquid saturation s is related to nucleate boiling and can be determined by analyzing vapor flow based on the Ergun equation [27 ] (see Note S3).A control volume, shown in Fig. 1 I, is adopted to depict liquid flows along the x-direction due to capi l lary wicking.The energy balance equation due to phase change can be written as q t (d x • 1) = h fg d ṁ v , where h fg is the latent heat of vaporization, and ṁ v is the mass flow rate of vapor.Combining the mass conservation as d ṁ l + d ṁ v = 0, where ṁ l ( x ) = ρ l ūl (x )δ l (x ) • 1 is the mass flow rate of the liquid, we obtain: Since all the liquid evaporated from the surface should enter at x = 0, the boundary condition can be set as: where L is the wicking length as shown in Fig. 1 G. Integrating Eq. ( 11) with the boundary condition of Eq. ( 12), we obtain the average liquid flow velocity as: Substituting Eq. ( 13) into Eq.( 10), the pressure gradient d P l / d x can then be obtained as: Here, δ l is related to the curvature of liquid meniscus H , which depends on liquid pressure P l and can be determined by the Young-Laplace equation as H = ( P v − P l )/2 σ [39 ], with P v being the vapor pressure and σ being the surface tension.Equation ( 14) is a non-linear first-order ODE for P l , which can be solved by the four-order Runge-Kutta method with boundary condition P l | x =0 = P sat .The liquid pressure P l and liquid film thickness δ l can be obtained accordingly (see Note S4).
With the curvature of the liquid meniscus determined, the heat transfer coefficient h e of evaporation atop the wick can then be calculated as: Here, A unit is the cross-sectional area of the one-unit cell for the wick, R tf is thermal conduction resistance through the thin-film region and R i, e is the thermal resistance for interfacial evaporation [28 ].The heat transfer coefficient h e of evaporation atop the wick is developed based on the saturated vapor environment, which is often seen in liquid film boiling experiments [16 ,21 ,23 -25 ].The detailed calculations of A unit , R tf and R i, e for different wicking structures can be found in Note S5.

Heat transfer coefficient of nucleate boiling inside the wick
For capi l lary-driven liquid film boi ling, there are two t ypical t ypes of microstructure: porous structures with interconnected micropores and micro-pi l lared surface.For the former, we use a representative micropore unit within the wick to analyze the nucleate boiling inside the wick.This micropore unit can be conceptualized as an annular pore with an effective pore radius of r eff , as shown in Fig. 2 A. This assumption is reasonable when the effective structural characteristics, including effective pore radius and porosity, are used to ensure the same capi l lary force and permeability for liquid wicking and bubble escaping.For the micro-pi l lared surface, we use a micropi l lar unit to analyze the nucleate boiling inside the wick.Since the analysis of both types of structures is similar, we will present here only the model derivation for the porous structures with interconnected micropores.The derivation for the micropillar can be found in Note S6.The micropore unit approximation for porous structures was also used in earlier studies [27 ] to analyze the nucleate boiling, where the nucleate boiling was assumed to be thinfilm evaporation with a uniform liquid film thickness δ tl independent of heat flux (Fig. 2 B).In Ref. [27 ], a constant liquid-vapor evaporation area is assumed, resulting in a constant heat transfer coefficient h bv as: Indeed, under higher heat flux, the effective evaporation area due to nucleate boiling should change with the heat flux due to the different fractions of bubbles activated.Similar to microlayer evaporation analysis for pool boiling widely used in the literature [29 ,32 ,35 ], we adopt microlayer evaporation to model nucleate boiling inside the wicking structures at various heat fluxes.As shown in Fig. 2 C, the heat within the micropore mainly transfers through a microlayer region characterized by a thickness of δ ml .For saturated water, the microlayer has a thickness δ ml that typically ranges between 1 and 10 μm [29 ], which is similar to the thin-film region for evaporation atop the w ick w ith a thickness < 0.15 r c [7 ], where r c is the characteristic radius of the solid skeleton of the wick (see Note S1).An averaged thickness δ ml of microlayer region is used here to estimate h bv by setting δ ml = 0.15 r c /2 in our calculations.
The nucleate boiling heat transfer can be considered as heat transfers through a microlayer conduction resistance R l from the wall to the liquid-vapor interface, then evaporates with an interfacial resistance R i, b , given by: where φ ml is the center angle corresponding to the microlayer region in the annular micropore, A ml is the liquid-vapor interfacial area of the microlayer region (or microlayer area), k l is the thermal conductivity of the liquid, and h i is the heat transfer coefficient for the liquid-vapor interface [40 ].Thus, h i can be obtained from the Schrage equation [41 ], with an accommodation coefficient of 0.04 chosen for water [28 ,40 ].The total heat dissipated by microlayer evaporation q b within the micropore can be expressed as Assuming the microlayer evaporation dissipates heat with an equivalent volume-averaged heat transfer coefficient h bv , q b can also be expressed as , where V up = V pore / ε w is the total volume of unit, V pore = π r eff 2 1 = r eff A total /2 is the volume of the micropore unit, and A total is the surface area.The volumetric heat transfer coefficient of nucleate boiling h bv can then be obtained as: Combining the analysis for the micro-pi l lared surface (see Note S6), a general form of h bv can be then obtained as: where χ is a structural factor: χ = −1 is used for the porous medium such as micropowder and micromesh, while χ = ε w /(1 − ε w ) is used for micro-pi l lared surface.The term A ml / A total represents the area fraction of the microlayer region to the total heated surface.As shown in pool boiling, the microlayer area increases with the heat flux [29 ,31 ,32 ,42 ], and an approximately linear function between A ml / A total and heat flux was recently proposed [32 ], i.e.A ml / A total ∝ q t .Since the microlayer evaporation in pool boiling and liquid film boiling are similar, this linear function between the microlayer area fraction ( A ml / A total ) and the total heat flux is also adopted in this work, i.e.A ml / A total = ηq t by introducing η as a scaling factor.Equation ( 20) can then be reformulated as: Equation ( 21 ) shows that larger heat flux q t , higher porosity ε w , and smaller effective pore radius r eff can lead to higher heat transfer coefficient of nucleate boiling h bv .In this work, the empirical parameter η is obtained by fitting the model predictions to the experimental results of liquid film boi ling, with the details provided in Section 3.

HTC and CHF of liquid film boiling
With the obtained h e and h bv , the temperature distribution of the wicking liquid film along the zdirection can be calculated from Eq. ( 7).Specifically, the wall temperature T w | z =0 can then be expressed as: Since h e varies w ith w icking distance x , an average wall temperature is adopted as: The HTC of liquid film boiling can then be calculated as We note that the previous model [27 ] simplified the liquid film thickness δ l as the wick thickness δ w (i.e.δ l ≈ δ w ) and ignored evaporation atop the wick, i.e. h e = 0.By applying these simplifications to Eq. ( 24), we have h t = mk eff tanh( m δ w ), which is the same result as derived in Ref. [27 ].When no bubbles nucleate (i.e.h bv = 0), Eq. ( 24) reduces to h t = 1/(1/ h e + δ w / k eff ), which is consistent with the HTC correlation developed for capi l lary evaporation [6 ].By such limit analysis, we can see that our model is more general and can be reduced to those in Ref. [27 ] and for capi l lary evaporation.
The CHF of liquid film boiling occurs when the liquid pressure drop is equal to the maximum capi l lary pressure P c, max in the wick [6 ], i.e. − L 0 (dP l / d x ) d x = P c,max .By combining this equation with Eq. ( 14), the relation of CHF for liquid film boiling can be obtained as: From Eq. ( 25), it is evident that increasing capi l lary force P c , max , increasing wick permeability K w , and augmenting relative liquid permeability K rl can enhance the CHF.Since K rl is an exponential function of the liquid saturation s [27 ,37 ], increasing the liquid saturation by promoting vapor escape can also enhance the CHF.Neglecting substrate friction (i.e.neglecting no-slip boundary condition and assuming uniform velocity along the z -direction at each x with ∂u l /∂z = 0 ), assuming uniform vapor distribution across the liquid film (i.e. the relative liquid permeability K r l to be constant), Eq. ( 25) reduces to q CHF = 2P c,max ρ l h f g δ w K rl K w /μ l L 2 , which is similar to the CHF correlation for liquid film boiling proposed by Zhang et al .[26 ].When simplifying liquid film thickness to become the wick thickness ( δ l ≈ δ w ), neglecting gravitational effect (i.e.g = 0), and assuming the K rl equals 1 (indicating that no bubbles exist inside the wick), Eq. ( 25) can be reduced to q CHF = 2P c,max ρ l h f g δ w K w [1 − tanh (λδ w ) /λδ w ] /μ l L 2 , which is consistent with the CHF correlation to predict the dry-out heat flux of capi l lary evaporation [43 ].In Note S7, we also compare the heat transfer performance between liquid film boiling and capi l lary evaporation using our model.

EXPERIMENTAL DETERMINATION OF THE SCALING FACTOR η
The scaling factor η in Eq. ( 21) characterizes the relationship between microlayer area fraction and heat flux for liquid film boiling, which is difficult to determine theoretically due to a lack of understanding of the bubble dynamics of liquid film boi ling.Instead, we theoretically analyze the η value of pool boiling in Note S8.Considering the difference in the transport processes between the pool boiling and liquid film boiling, we here obtain the η value by fitting our model predictions with our experimental data on sintered multi-layer micromesh surfaces.In this work, the physical properties of the micromesh, such as porosit y, permeabilit y, and wettability, which are not fully presented in the existing literature, are characterized (see Note S9).Liquid film boiling experiments are conducted using a custom-made experimental system [16 ,28 ], which provides a saturated vapor environment for liquid film boi ling (Fig. 3 A).The tested sample is mounted vertically to the heating block with a 10 × 10 mm 2 heating area and immersed in degassed deionized water with one-side liquid supply (see Note S10), consistent w ith our prev ious works [16 ,20 ,28 ].The heat flux q t and surface temperature T w | z =0 are obtained from the onedimensional temperature distribution in the heating block.The measured HTC is calculated as The experimental data on micromesh samples s1 and s2, which have similar wire diameter ( ∼50 μm) and thickness ( ∼265 μm) but different spacing width (77 and 160 μm, respectively), is used to obtain the η value.The model prediction of HTC is obtained from Eq. ( 24) based on the structure properties as shown in Table S5.To maintain consis-tenc y w ith the experiment, we include an extra heat conduction resistance of the substrate in the modelpredicted HTC calculation, as 1/(1/ h t + δ sub / k sub ), where h t is calculated by Eq. ( 24), δ sub is the thickness of the substrate in the experiment, and k sub is the thermal conductivity of the substrate.This calculation method is also used in the following sections when comparing model-predicted HTC with experimental data.The factor η in Eq. ( 21) is determined to be η = 2.15 × 10 −3 cm 2 W −1 by fitting the model prediction of HTC with our experimental data on such uniform micromesh samples, with a mean absolute percentage error (MAPE) of 5.19% (Fig. 3 B).Such a unified value of η could then be used to estimate the microlayer area fraction of the uniform wicking surface according to A ml / A total = ηq t .Figure 3 C shows the model-predicted CHF for different uniform micromesh samples, with thickness ranging from 260 to 377 μm, porosity ranging from 0.57 to 0.76, and spacing width ranging from 77 to 205 μm.With the same value of η = 2.15 × 10 −3 cm 2 W −1 , our CHF model agrees well with experimental data, with ±15% accuracy (Fig. 3 C).
Figure 3 D-F show the effects of thickness δ w , porosity ε w , and spacing width s w on CHF and maximum HTC (HTC at CHF).The CHF is enhanced with increased micromesh thickness (Fig. 3 D), due to enhanced wick permeability and decreased resistance of the liquid flow [6 ].However, increasing the wick thickness could also increase the vapor escaping resistance and sightly decrease the CHF.Although increasing the porosity of the micromesh also increases the liquid wicking velocity, there exists an upper limit for the wick porosity for the micromesh with the same wire diameter and the same spacing width (Fig. 3 E).Increasing the spacing width enhances wick permeability, however, the maximum capi l lary force is reduced due to the enlarged effective pore radius [44 ], thereby yielding an optimal spacing width s w for CHF (Fig. 3 F).The maximum HTC is mainly controlled by the maximum microlayer area A ml, max for nucleate boiling, calculated as A ml, max = 2.15 × 10 −3 q t A total , where is the total area of the micropore with N being the number of micropores and V total being the total wick volume.As shown in Fig. 3 D and E, the maximum HTC is increased with larger thickness and increased porosity since the maximum microlayer area for evaporation is expanded with an increased CHF q CHF and the expanded total micropore area A total .The effect of spacing width s w on the maximum HTC is multifaceted.Though increasing s w augments porosity to expand the total micropore area A total , it simultaneously increases the effective pore radius r eff and may reduce the total micropore area A total .Besides, q CHF initially increases with spacing width and subsequently decreases, which also accordingly changes A ml, max .Consequently, the competition results in an optimal s w to maximize effective h bv for nucleate boiling and HTC (Fig. 3 F).From Fig. 3 D-F, we can clearly see that our model with the same value of η = 2.15 × 10 −3 cm 2 W −1 can predict both maximum HTC and CHF with high fidelity for micromesh with different structural parameters.Although changing the wick geometry changes both the total surface area A total and the effective microlayer region area A ml , the area fraction A ml / A total is effectively estimated by the scaling factor η in our model.

Model validation
Figure 4 A compares the experimental data on micropowder surface ( d w = 250 μm, δ w = 900 μm and ε w = 0.64) from Weibel et al .[24 ], our ex-perimental data on micromesh of sample s2 ( d w = 50 μm, s w = 160 μm, δ w = 264 μm and ε w = 0.71), and our model predictions for HTC using Eq. ( 24) on the same sample.The modelpredicted HTC from the Sudhakar et al .model [27 ] for our micromesh sample s2 is also presented in Fig. 4 A. It is shown that the model by Sudhakar et al .[27 ] predicts a constant HTC, due to its underlying assumption of a constant evaporative area inside the micropore for nucleate boiling and adiabatic condition atop the wick (i.e.h e = 0 due to the neglect of evaporation atop the wick).In contrast, our model can predict the increasing trend of HTC with increasing heat flux, which agrees well with results from previous liquid film boiling experiments [16 ,18 ,20 -24 ].This is because our model considers both the evaporation atop the wick and the variation in the nucleate boiling inside the wick under varying heat flux.
To validate our model and understand the coupling roles between evaporation atop the wick and nucleate boiling inside, we develop two simplified models (RM) based Eq. ( 24 ).RM1 is a simplification with the approximations of Sudhakar et al .[27 ] with h e = 0, δ l = δ w and ηq t = 0.9, and RM2 is a reduced model by setting h e = 0 (neglecting evaporation atop the wick).The value ηq t = 0.9 in RM1 is derived from the assumption according to Ref. [27 ], where evaporation occurs uniformly at the entire liquid-vapor interface (i.e.φ ml = 2 π ), and the microlayer thickness is assumed as 0.1 times the effective pore radius (i.e.δ ml = 0.1 r eff ), giving ηq t = A ml / A total = [( r eff − δ ml ) φ ml ]/(2 π r eff ) = 0.9.As shown in Fig. 4 A, RM1 can predict nearly identical results with Ref. [27 ], and a constant HTC is obtained.At low heat flux, the predictions of HTC from RM2 deviate from our original model.This is because at low heat flux, the heat transfer fraction of evaporation atop the wick is relatively high and RM2 neglects the contribution of the evaporation atop the wick.In our model, the relationship between evaporation atop the wick and nucleate boiling inside is self-consistent.The heat dissipation fraction of evaporation atop ( q e /q t ) from our model decreases as the heat flux fraction ( q t /q CHF ) increases (Fig. 4 B).The evaporation heat dissipation fraction exceeds 20% when heat flux fraction ( q t /q CHF ) is below 30%, which indicates that evaporation atop the wick should be considered in liquid film boiling heat transfer modeling, and its neglection may cause a large deviation, especially for low heat fluxes.
As shown in Fig. 4 C, our model also exhibits good agreement with experimental data of sample s2 when predicting CHF.To further validate our model prediction, we develop RM3 based on our CHF relation of Eq. ( 25) by neglecting the friction of substrate  [24 ].Model prediction based on Sudhakar et al .[27 ] for our sample s2 measurement is also presented.A simplified model, RM1 (with h e = 0 and ηq t = 0.9), is developed based on our model by setting consistent conditions with the assumptions of Sudhakar et al .[27 ].Another simplified model, RM2, is developed based on our model with the assumption of h e = 0.The heat flux fraction q t /q CHF is adopted since the CHF is varied for different microstructures in the literature.(B) The fraction of the evaporation heat flux atop the wick to the total heat flux q e /q t as a function of heat flux fraction q t /q CHF .(C) Comparison of the model-predicted CHF with the experimental data and predictions from other models.The red solid line, purple triangle, and blue diamond represent our experimental result of sample s2, model prediction from Sudhakar et al .[27 ], and model prediction from Zhang et al .[26 ], respectively.(D) Comparison of measured HTC on micro-pillared surface with the predictions of this work and from pool boiling models [32 ,45 ]. (E) Comparison of measured CHF on micro-pillared surface with the predictions from this work and pool boiling models [10 ,46 ,47 ].The red circles in (D) and the red solid line in (E) are the experimental data on micropillar ( d w = 60 μm, δ w = 320 μm and ε w = 0.6) from Cai and Bhunia [25 ].
(i.e.liquid wicking velocity is uniform along the zdirection at each x with ∂u l /∂z = 0 ) and neglecting the gravitational effect (i.e.g = 0), to be consistent with the case in Ref. [27 ].Notably, the CHF prediction from RM3 is very close to the model by Sudhakar et al. [27 ].Besides, a simplified model (RM4) is developed with ∂u l /∂z = 0 , g = 0, and K rl = 0.15, according to Zhang et al .[26 ] and, as expected, the prediction is consistent.
Here, we also compare our model on liquid film boiling with the models developed for pool boiling.Since most pool boiling models are developed based on the micropi l lar, we select micropi l lars as the micro-structured surface to compare the predictions.The micropi l lar used in the modeling has the following geometric parameters: a pi l lar diameter of 60 μm, a porosity of 0.6, a thickness of 320 μm, and a wicking length of 2 mm, which are same as those used in the experiment conducted by Cai and Bhunia [25 ].As shown in Fig. 4 D, our model prediction for HTC shows better agreement with the experiment compared to the HTC models developed for pool boiling in literature [32 ,45 ].This is because natural convection contributes significantly to pool boiling  [25 ,48 ], micropowder [21 ,24 ], and micromesh [16 ,23 ]. (A) h t , model vs. h t , exp .(B) q CHF ,model vs. q CHF ,exp .
heat transfer, but is not as important in liquid film boiling due to the thinness of the capi l lary film.Although there are many CHF models for pool boiling that consider surface properties and liquid wicking [10 ,46 ,47 ], they are not suitable for predicting CHF in liquid film boiling (Fig. 4 E), due to the different mechanisms aforementioned.For pool boiling, the relatively large bubbles easily coalesce to form a vapor blanket, leading to the thermal-hydraulic CHF.However, in liquid film boiling CHF occurs because of surface dry-out, which is caused by capi l lary wicking failure.

Model prediction
In this section, we use our model to predict both HTC and CHF with unified microlayer area factor η = 2.15 × 10 −3 cm 2 W −1 for liquid film boi ling on various types of uniform micro-structured surfaces.The wicking structures include silicon micropi l lar arrays [25 ,48 ], packed copper micropowders [21 ,24 ], and sintered copper micromeshes [16 ,23 ] (Table 1 ).The liquid supply methods for these experiments include one-side, two-side, and all-around, with a heated area ranging from 4 to 100 mm 2 , a wick thickness ranging from 200 to 1200 μm, and a wick porosity ranging from 0.51 to 0.76. Figure 5 A shows that the HTC predictions from our model are in good agreement with experimental results for various uniform wicking structures, achieving a MAPE of 13.7%.For the HTC ranging from 3 0-3 00 kW m −2 K −1 , the model predicts the experimental data well, with an accuracy of ±30% due to its ability to calculate the heat dissipated by both nucleate boiling inside the wick and evaporation atop the wick.The experimental CHF with a range of 40-1200 W cm −2 could also be well predicted by our model with an error of 20% (Fig. 5 B).The MAPE for CHFs of all samples is calculated to be 11.9%.We note that HTC and CHF are affected by many aspects, including surface geometries, thermal conductivity, liquid supply methods, surface wettability, and total heat flux input.Our model has taken these aspects into account and experimentally determines a factor η to model the variation in microlayer area fraction during liquid film boiling.Therefore, our model can predict both the HTC and CHF for various wicking structures within a spread of ±30%.

CONCLUSION
In this work, a heat transfer model for predicting both the CHF and HTC is developed for liquid film boiling on various micro-structured surfaces.The model accounts for both evaporation atop the wick and nucleate boiling inside the wick.Evaporation atop the wick is calculated by analyzing the thermal resistance network at the thin-film region.Nucleate boiling is modeled using microlayer evaporation theory with an empirical factor η to describe the relationship between microlayer area fraction and heat flux as A ml /A total = ηq t .The scaling factor η is found to be independent of structural parameters and the value is determined to be η = 2.15 × 10 −3 cm 2 W −1 by our liquid film boi ling experiments.Our model can be reduced to prev ious models w ith simplifications.With the same value of η, the model predictions of both HTC and CHF are in good agreement with the reported experimental data for various uniform micro-structured surfaces, including micropillar, micropowder, and micromesh, with a spread of ±30%.This work provides a tool for designing micro-structured surfaces for advancing thermal management applications.

Figure 2 .
Figure 2.Pore-scale analysis of nucleate boiling inside the wicking structure.(A) Nucleate boiling inside the wick micropores, which can be equivalent to annular interconnected micropores[27 ]. (B) Nucleate boiling in a micropore unit is modeled as evaporation with uniform thickness of the thin liquid film in the previous model[27 ]. (C) Nucleate boiling in a micropore unit is modeled as microlayer evaporation at the microlayer region with an average thickness of δ ml in this work.

Figure 3 .
Figure 3. Determination of effective microlayer evaporation factor η by experimental results.(A) Schematic of the custom-made experimental setup for liquid film boiling measurement.(B) Determination of η = 2.15 × 10 −3 cm 2 W −1 with our experimental data on copper micromesh samples s1 and s2.(C) Comparison of experimental CHF and model-predicted CHF with η = 2.15 × 10 −3 cm 2 W −1 .The CHF and the maximum HTC as the function of (D) thickness δ w , (E) porosity ε w , and (F) spacing width s w .The red circles in (D-F) represent the experimental data and the black solid lines are the modeling results with η = 2.15 × 10 −3 cm 2 W −1 .

Figure 4 .
Figure 4. Model validation with experimental data and simplifications.(A) Comparison of model-predicted HTC with the experimental data and prediction from other models[27 ].The experimental data on micromesh sample s2 ( d w = 50 μm, s w = 160 μm, δ w = 264 μm and ε w = 0.71) are collected in this work, and the experimental data on micropowder ( d w = 250 μm, δ w = 900 μm and ε w = 0.64) are taken from Weibel et al .[24 ].Model prediction based on Sudhakar et al .[27 ] for our sample s2 measurement is also presented.A simplified model, RM1 (with h e = 0 and ηq t = 0.9), is developed based on our model by setting consistent conditions with the assumptions of Sudhakar et al .[27 ].Another simplified model, RM2, is developed based on our model with the assumption of h e = 0.The heat flux fraction q t /q CHF is adopted since the CHF is varied for different microstructures in the literature.(B) The fraction of the evaporation heat flux atop the wick to the total heat flux q e /q t as a function of heat flux fraction q t /q CHF .(C) Comparison of the model-predicted CHF with the experimental data and predictions from other models.The red solid line, purple triangle, and blue diamond represent our experimental result of sample s2, model prediction from Sudhakar et al .[27 ], and model prediction from Zhang et al .[26 ], respectively.(D) Comparison of measured HTC on micro-pillared surface with the predictions of this work and from pool boiling models[32 ,45 ]. (E) Comparison of measured CHF on micro-pillared surface with the predictions from this work and pool boiling models[10 ,46 ,47 ].The red circles in (D) and the red solid line in (E) are the experimental data on micropillar ( d w = 60 μm, δ w = 320 μm and ε w = 0.6) from Cai and Bhunia[25 ].

Table 1 .
Summary of the experimental conditions and sample characteristics in the literature.
Figure 5.Comparison between model predictions and experimental data on various types of wicking structures: micropillar