CFD Analysis of Phase-Change Heat Transfer in Heat Pipes Using a Combined Boiling Model (CBM)

1. Introduction

Cooling systems used in modern applications are becoming increasingly complex, reflecting the growing sophistication of the components they serve. Adequate and efficient heat dissipation is essential for ensuring the reliable operation of devices such as electric motors, hydrogen fuel cells, control units, and motherboards, batteries [1,2], particularly in industries like automotive and aerospace. As the demand for compact and powerful electronic systems rises, traditional single-phase cooling approaches are no longer sufficient. Instead, modern thermal management relies on processes that combine heat transfer with phase-change phenomena.

Among these, heat pipes have been studied extensively due to their ability to transfer large heat loads with minimal temperature gradients. First investigated in the 19th century [3] and formalized by Grover in 1963 [4], heat pipes exploit evaporation and condensation to achieve exceptionally high effective thermal conductivity. Predicting their thermal–fluid behaviour, however, is challenging due to the interplay of multiphase flow, compressibility, and conjugate heat transfer.

Early analytical models treated heat pipes as one-dimensional conduction systems with simplified assumptions for fluid transport [5], but modern studies increasingly employ computational fluid dynamics (CFD) to resolve vapour flow, liquid film dynamics, and interfacial phase change [6].

In its simplest form, a heat pipe is a closed thermosiphon as shown in Figure 1.

Despite their apparent simplicity as closed thermosiphons, accurate simulations must account for nucleate boiling, condensation, and evaporation processes. Operation in the nucleate boiling regime is desirable for high heat transfer, yet the narrow temperature window (e.g., ΔT ≈ 1–30 °C for water at atmospheric pressure) imposes strict modelling requirements [7,8].

Various CFD studies on subcooled wall boiling have been reported in the literature [9,10,11]. These reviews indicate that researchers primarily employ interfacial surface tracking methods, such as Volume-of-Fluid (VOF), or the Euler–Euler framework to model boiling phenomena. VOF enables detailed bubble- and droplet-scale studies but is computationally expensive [12,13,14,15] , whereas Euler–Euler, based on averaged field equations, allows for coarser grids and larger time steps, making it computationally more efficient and thus more suitable for practical engineering applications.

Within the latter, the widely used RPI wall-boiling model [16] often requires empirical recalibration and may under predict flow boiling [17,18,19,20]. Improvements are needed, particularly for heat pipes where circulation arises from combined pressure and gravity effects.

Evaporation and condensation models are also essential for capturing vapour transport, with Lee’s approach [13,21,22,23,24] showing strong agreement with experiments [8]. For realistic systems, conjugate heat transfer (CHT) between fluid and solid domains must be included, though only a few studies have addressed multiphase CHT comprehensively [11,25].

In previous work, we proposed the Combined Boiling Model (CBM) and validated it in a simplified CHT case [26]. Here, we extend its application to heat pipes, focusing on phase-change heat transfer under realistic operating conditions.

2. Model Description

2.1. Background

The CBM model was developed and preliminary validation on a simplified case employed in our previous study [26] . That case represented an open system with a velocity inlet and a pressure outlet boundary condition, consisting of two domains: the fluid region and a solid heating plate. A detailed grid-independence and time-step sensitivity analysis was performed, from which the optimal numerical parameters were established. These parameters were subsequently applied to the more complex heat pipe model investigated in the present work.

In contrast to the Robinson-type duct of the previous study, the present model represents a closed vertical heat pipe filled with working fluid. It comprises three main components: the fluid domain of the pipe, the solid copper casing, and the coolant region at the condenser side Figure 2.

Geometrically, this configuration is not significantly more complex than the Robinson test duct; however, from a numerical perspective it is considerably more challenging. In particular, the closed domain requires consideration of vapour compressibility. Furthermore, bubble dynamics are no longer dictated by inlet velocity, but instead result from buoyancy forces, pressure gradients, and bubble drag within the flow, all of which must be explicitly resolved.

In addition, heat pipe simulations must account for condensation processes. For this reason, a simplified film-condensation model has been incorporated into the CBM. The key extensions relative to the original analysis are summarized as follows:

• closed system,
• weakly compressible vapour flow,
• inclusion of condensation,
• three interacting domains: the solid domain (single-phase, energy equation), the heat pipe fluid domain (multiphase), and the coolant domain (single-phase flow).

2.2. Governing Equations

The Euler–Euler multi-fluid framework treats each phase as an interpenetrating continuum, with interphase exchange terms added to the conservation laws [27]. Ensemble averaging removes microscopic interfaces and yields macroscopic equations analogous to single-phase flow, extended by phase volume fractions and interfacial transfer terms [28,29,30]. Based on the formulation of Drew and Passman [31], the governing equations are:

Continuity equation (phase k):

$${\displaystyle \frac{\partial {\alpha }_k{\rho }_k}{\partial t}}+\nabla \cdot {\alpha }_k{\rho }_k{v}_k={\displaystyle \sum _{l=1,l\ne k}^N{\Gamma }_{kl}}, k=1,\dots N$$

(1)

with the compatibility condition

$${\displaystyle \sum _{k=1}^N{\alpha }_k=1}.$$

(2)

Momentum equation (phasek):

$$\begin{array}{l}{\displaystyle \frac{\partial {\alpha }_k{\rho }_k{v}_k}{\partial t}}+\nabla \cdot {\alpha }_k{\rho }_k{v}_k\otimes v_k=\\ -{\alpha }_k\nabla p+\nabla \cdot {\alpha }_k\left({\tau }_k+T_k^t\right)+{\alpha }_k{\rho }_{k}f+{\displaystyle \sum _{l=1,l\ne k}^N{\mathrm{M}}_{kl}}+{\displaystyle \sum _{l=1,l\ne k}^N{v}_{ki}{\Gamma }_{kl}}, k=1,\dots N\end{array}$$

(3)

where f includes gravity and inertial forces. Pressure is assumed equal for all phases:

$$p=p_k, k=1,\dots N$$

(4)

Energy equation (mixture form):

$$\begin{array}{c}{\displaystyle \frac{\partial {\alpha }_k{\rho }_k{h}_k}{\partial t}}+\nabla \cdot {\alpha }_k{\rho }_k{\mathrm{v}}_k{h}_k=\nabla \cdot {\alpha }_k\left({\mathrm{q}}_k+{\mathrm{q}}_k^t\right)+{\alpha }_k{\rho }_k{q^{‴}}_k+{\alpha }_k{\rho }_k\mathrm{f}\cdot {\mathrm{v}}_k+\\ \nabla \cdot {\alpha }_k\left({\tau }_k+{\tau }_k^t\right)\cdot {\mathrm{v}}_k+{\alpha }_k{\displaystyle \frac{\partial p}{\partial t}}+{\displaystyle \sum _{l=1,l\ne k}^N}{\mathrm{H}}_{kl}+{\displaystyle \sum _{l=1,l\ne k}^N{h}_{ki}}{\Gamma }_{kl}\end{array}$$

(5)

with h enthalpy, q_k conductive flux, and H_kl interfacial energy exchange. Liquid and vapour are assumed to be in local thermal equilibrium [29]. A closed system is obtained once interfacial exchange terms and boundary conditions are specified [27,30].

2.3. Interfacial Momentum Exchange

The interfacial momentum source term accounts for drag, turbulent dispersion, lift, and wall lubrication:

$$M_{kl}=F_D+F_{TD}+F_L+F_{WL}.$$

(6)

The drag contribution is given by

$$F_D={\displaystyle \frac{3}{4}}{C}_D{\rho }_c{\alpha }_d{A^{‴}}_i\hspace{0.33em}\left|v_r\right|v_r$$

(7)

where v_r = v_d - v_c.

In AVL FIRE™, the latest Gas-Liquid System 3 [27,32] is used. It blends bubbly and droplet drag correlations based on dispersed-phase volume fraction, using Schiller–Naumann [33] or Tomiyama [34] models depending on bubble size and Eötvös number. Turbulent dispersion accounts for vapour diffusion in mixing flows.

2.4. Wall Boiling Model

The near-wall heat transfer is modelled using the RPI framework of Kurul and Podowski [16], extended for forced convective boiling following Gilman [19], with further simplifications by Shi and Zhang [20]. The complete formulation, including detailed derivations, was presented in our previous work [26], therefore only the governing partitioning and key correlations will be presented in this section .

The total wall heat flux is expressed as:

$${\dot{q}}_w={\dot{\mathrm{q}}}_{\mathrm{conv}}+ {\dot{\mathrm{q}}}_{\mathrm{sta}}+{\dot{\mathrm{q}}}_{\mathrm{sli}}+{\dot{\mathrm{q}}}_{\mathrm{eva}}$$

(8)

The evaporative component is:

$${\dot{q}}_{eva}=\dot{m}{h}_{fg}, \dot{m}={\displaystyle \frac{\pi }{6}}{D}_l^3{N}_w{\rho }_{v}f$$

(9)

with bubble departure frequency f [35], nucleation site density N_w [36], and departure diameter D_l [37] obtained from standard correlations.

The convective contribution is:

$${\dot{q}}_{conv}=C_m{h}_{conv}\left(T_w- T_l\right)\left(1-A_b\right), A_b=A_{sta}+ A_{sli}$$

(10)

Stationary and sliding bubble quenching fluxes are:

$${\dot{q}}_{sta}=h_{quen}\left(T_w- T_l\right)A_{sta}, {\dot{q}}_{sli}=h_{quen}\left(T_w- T_l\right)A_{sli}$$

(11)

where h_quen is the transient conduction coefficient [38].

The sliding bubble model accounts for bubble nucleation, growth, sliding, merging, and eventual lift-off along the heated wall (Figure 3).

Bubble motion is governed by the balance of buoyancy, pressure-difference, gravitational, and drag forces. Nucleated bubbles first grow and slide along the wall, during which they exchange heat with the liquid and contribute to the sliding heat flux. As they move, bubbles may coalesce with neighbouring bubbles at a distance defined by the nucleation site spacing g_s, leading to accelerated growth. The sliding distance l_sl therefore depends on the number of merged bubbles n_nr and the remaining distance before lift-off [19,20].

Lift-off occurs once the bubble reaches the departure diameter D_l, determined by buoyancy–drag balance and surface tension constraints. The corresponding sliding velocity is obtained iteratively since the drag coefficient C_d depends on the bubble Reynolds number. In this way, the model explicitly captures the four key stages of bubble dynamics: (1) nucleation, (2) sliding and growth, (3) coalescence, and (4) lift-off.

2.5. Film Condensation Model

A simplified film condensation model based on Nusselt theory for laminar condensation on vertical walls [39] was applied in the present investigation. The model assumes: T_w < T_sat, a liquid film is always present, liquid velocity is negligible, and all wall heat flux is converted into condensation.

The film thickness is approximated as:

$$h_l={\displaystyle \frac{{\alpha }_l{V}_P}{A_b}}$$

(12)

where α_l is the liquid volume fraction, V_P is the cell volume, and Ab is the cell wall area exposed to condensation.

The temperature gradient across the film is:

$${\displaystyle \frac{\partial T}{\partial x}}\cong {\displaystyle \frac{\left(T_{sat}-T_w\right)}{h_l}}$$

(13)

The wall heat flux is:

$$q_w=A_b{k}_l{\displaystyle \frac{\partial T}{\partial x}}$$

(14)

where kl is the liquid thermal conductivity.

The corresponding cell condensation rate is:

$${\dot{m}}_{cond}={\displaystyle \frac{A_b{k}_l\left(T_{sat}-T_w\right)}{h_l{h}_{fg}}}$$

(15)

where h_fg is the latent heat of vaporization.

2.6. Flow Evaporation and Condensation

Phase change in the bulk flow is modelled using the simplified approach of Lee [40], which assumes saturation-temperature equilibrium and employs empirical coefficients for calibration [40,41,42]. The corresponding mass source terms are:

$${\mathrm{Г}}_l=-{\mathrm{Г}}_v=f_c{\alpha }_v{\rho }_v{\displaystyle \frac{\left(T-T_{sat}\right)}{T_{sat}}} for condensation (T < T_{sat})$$

(16)

$${\mathrm{Г}}_v=-{\mathrm{Г}}_l=f_e{\alpha }_v{\rho }_v{\displaystyle \frac{\left(T-T_{sat}\right)}{T_{sat}}}for evaporation(T>T_{sat})$$

(17)

The corresponding energy source terms are:

$$H_c={\mathrm{Г}}_M\Delta h\mathrm{f}\mathrm{o}\mathrm{r} \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}(T < T_{sat})$$

(18)

$$H_e={-\mathrm{Г}}_M\Delta h\mathrm{f}\mathrm{o}\mathrm{r} \mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\left(T>T_{sat}\right)$$

(19)

where Γ_M is the phase-change mass flux and Δh the latent heat. Empirical coefficients f_c and f_e were set to 0.1, following recommendations from Lee [40] and subsequent studies [41,42].

2.7. Fluid–Solid Conjugate Heat Transfer

Conjugate heat transfer (CHT) between solid and fluid domains is modelled using an interfacial diffusion coefficient. The coefficient is derived from the heat flux ${\dot{q}}_w$, obtained from the modified RPI model, and the interface temperature difference ${∆T}_{int}=\left(T_w- T_l\right)$:

$$D={\displaystyle \frac{{\dot{q}}_{w}A}{{∆T}_{int}{c}_{p,l}}}$$

(20)

where A is the heated surface area and c_p,l the liquid specific heat capacity.

The computed diffusion coefficient is assigned to the fluid domain to evaluate the local heat transfer coefficient. This is then imposed as a boundary condition on the solid domain, where the wall temperature T_w is recalculated. The updated T_w is in turn enforced at the fluid side, ensuring coupled exchange. The process is repeated iteratively within each time step until thermal equilibrium is achieved.

3. Computational Model

The updated Combined Boiling Model (CBM) simulations were compared against experimental data reported by Fadhl, Wrobel, and Jouhara [12] and detailed in Fadhl’s dissertation [8]. The dataset provides thermal conditions in a closed thermosyphon at various power inputs; in the present study, two representative operating conditions were selected. The experimental setup is shown schematically in Figure 4.

Eight thermocouples were mounted along the outer wall of the heat pipe, with measurement positions distributed from the condenser to the evaporator section.

3.1. Computational Mesh and Time-Step

Based on the test rig geometry, a multi-domain block-structured mesh was created, comprising 161.280 cells (Figure 5).

The fluid domain contained 56.000 cells (blue), the solid copper domain 59,200 cells (orange), and the coolant domain 46.080 cells (grey).

A full mesh/time-step sensitivity study was not repeated here, as this was performed in earlier work [26]. That study showed that meshes finer than 1 mm adequately resolved vapour distribution, with interface temperature deviations within ±5 percentage of experiments. Wall resolution yielded low y+ values (8–20), which, although below the values often suggested for Euler–Euler boiling simulations (y+ > 40 [43]), ensured consistent temperature prediction. A 1.0 mm mesh with Δt = 0.005 s was identified as the optimal configuration.

In the present work, these settings were adopted with additional near-wall refinement (0.4 mm boundary layer and five solid layers of 0.18 mm each) to capture condensation films without artificial breakup in the Euler–Euler framework.

3.2. Numerical Settings

All simulations were performed in transient mode using the Euler–Euler framework with the models described above. Three computational domains were considered: fluid, solid, and coolant. In the solid domain, only the energy equation was solved, while the fluid domain included continuity, momentum, turbulence, energy, and volume fraction equations. The Combined Boiling Model (CBM) was applied to capture near-wall phase change. Key simulation parameters for the multiphase fluid domain are summarized in Table 1.

The interfacial exchange models (Table 2) employed the CBM for mass and energy transfer, “Gas–Liquid System 3” for momentum exchange, Tomiyama’s correlation for bubbles (1 mm) and Schiller–Naumann for droplets (0.5 mm).

A sensitivity study on bubble and droplet departure parameters, was not conducted in present study those the default values provided by AVL FIRE™ were retained.

Water properties were taken from the software database; vapour was modelled as a weakly compressible ideal gas. Phase change was restricted to the saturation temperature, with vapour pressure following the standard saturation curve (Figure 6).

Solver stability and accuracy were maintained using the SIMPLE algorithm for pressure–velocity coupling. Convergence criteria were set to 1 × 10⁻⁴ for continuity, energy, and volume fraction, and 5 × 10⁻⁴ for momentum and turbulence. The MINMOD differencing scheme with 40% upwind blending [44] was used on momentum, volume fraction and energy. Using these settings, the number of iterations fluctuated between 5 and 20 per time step, with higher values occurring during boiling or condensation onset.

3.3. Initial and Boundary Conditions

Boundary conditions were specified for all three computational domains. In the condenser heat exchanger domain, a water inlet mass flow rate of 0.003 kg/s at 293.75 K was imposed, while the outlet was set to ambient pressure (Figure 5). The external walls were treated as adiabatic. The solid domain was defined as copper, and a prescribed heat rate boundary condition was applied at the evaporator section (Figure 5). Two heat load cases were investigated, 173 W and 376 W. The heat pipe fluid domain was initialized with a mixture of liquid water and water vapour, with the evaporator section fully filled with liquid water, corresponding to a filling ratio of 1, and the upper region filled with water vapour.

4. Results and Discussion

4.1. Initial Flow Simulation

The simulations were conducted in two stages. In the first stage, the coolant flow field was solved until a quasi-steady state was reached. This solution was subsequently used to initialize the full conjugate heat transfer (CHT) model, thereby ensuring a physically consistent start-up condition for the coupled phase-change simulation.

The computed coolant velocity field is illustrated in Figure 7 (a), where the colour contours represent the velocity magnitude and the streamlines indicate the flow direction. The coolant flow simulation was performed for a total duration of 60 s; however, a quasi-steady state was reached after approximately 30 s, as confirmed by the outlet mass flow rate presented in Figure 7 (b).

This preliminary step provided a stable and realistic hydrodynamic field for the subsequent CHT simulations, minimizing numerical transients during model initialization.

4.2. CHT Full-Model Simulation

In current preliminary validation study five CHT simulations were performed and analysed, of which three were compared directly with available experimental and numerical data [12]. The coupled vapour–liquid flow dynamics developed primarily during the initial transient period, while the overall thermal response required a longer time to stabilize.

The transient development of the liquid and vapour phases is best observed for the 100% filling ratio, in which the entire evaporator section was initially filled with liquid water, and the remaining adiabatic and condenser sections contained pure vapour. The evolution of the liquid volume fraction for two representative heat loads (173 W and 376 W) is shown in Figure 8.

At the beginning of the simulation, a rapid evaporation rate was observed in both cases, stabilizing after approximately 10 s. The larger vapour generation at higher heat input (376 W) became evident after this period. The formation of large vapour pockets was observed in the adiabatic section before collapsing into the condenser region. Within the condenser, the simulations consistently predicted the formation of a thin liquid film on the upper wall surface, as illustrated in Figure 9.

The fluid flow in this region clearly show vapour transport toward the upper condensation surface, followed by flow reversal and downward liquid-film drainage. This behaviour is consistent with experimental observations and confirms that the present model is capable of capturing the essential evaporation–condensation mechanisms occurring within the heat pipe.

4.3. Preliminary Model Validation

Experimental validation data were available for a 50% fill ratio [8], corresponding to half of the evaporator section being initially filled with liquid water , while the remaining volume contained vapour. Three heat-input cases were investigated: 173 W, 225 W, and 376 W.

Due to the low amount of liquid phase, the closed thermosyphon exhibited slightly different behaviour compared to the 100% fill-ratio case discussed earlier. The flow evolution during the initial transient period for the two representative heat inputs (173 W and 376 W) is illustrated in Figure 10.

At start-up (up to approximately 8 s), both cases showed similar evaporation patterns within the evaporator section. However, after a longer time (approximately 30 s), the 376 W case displayed significantly less residual liquid in the evaporator section and a noticeably thicker condensate film in the condenser section. The stronger evaporation and higher vapour production observed at 376 W are consistent with the expected temperature increase under higher thermal loading. The influence of heat input on the overall temperature distribution is shown in Figure 11.

The temperature gradients within the heat pipe are clearly visible, with the highest temperatures occurring in the evaporator section at 376 W. A slight flow asymmetry was observed, resulting from the lower position of the coolant inlet. The colder coolant entered from the bottom, absorbed heat as it travelled through the condenser, and gradually increased in temperature before exiting from the opposite upper side.

Quantitative comparisons further support the validity of the proposed CBM (Combined Boiling Model) approach. The results show improved agreement with experimental measurements and with previous numerical studies using the VOF–Lee model [8]. Temperatures were monitored at eight locations, as shown in Figure 4, and the computed values were compared with both experimental data and Lee model predictions.

4.3.1. Heat Input: 173 W

At 173 W, deviations in the condenser remained below 3.5%, despite the use of a simplified condensation model. The largest improvement compared to the Lee results occurred in the evaporator, where deviations were reduced to 1.27% and 1.79%, compared with 9.42% and 12.14% in the previous model. In the adiabatic section, the deviation decreased from 10.68% (Lee) to 3.13% (CBM).

4.3.2. Heat Input: 225 W

At 225 W, the Lee model achieved slightly better agreement in the condenser (errors typically below 1%), whereas the CBM model exhibited deviations between 2–4% in this region. However, the CBM model significantly outperformed the Lee model in the evaporator and adiabatic sections, where deviations were reduced to 2.30%, 1.66%, and 3.38%, compared with 2.72%, 10.49%, and 10.32% for the Lee model.

4.3.3. Heat Input: 376 W

At 376 W, the CBM model again demonstrated better agreement in the evaporator and adiabatic regions compared to the Lee model, although slightly larger deviations were observed in the condenser section, where the Lee simulation produced marginally closer results.

4.4. Overall Model Performance

The overall performance of the model, together with its comparison against experimental data and the Lee model results, is presented in Figure 12. Although detailed deviations from the experimental measurements were provided in the previous tables, the overall agreement can be more clearly visualized by examining the trend line curves. The markers in the figure correspond to the individual measurement points, while the solid lines represent the fitted trend lines. From these results, it can be concluded that the proposed Combined Boiling Model (CBM) provides a markedly improved prediction of the thermal behaviour of the analysed heat pipe across all investigated heat loads.

While the condensation sub-model integrated within the CBM framework remains simplified, the achieved level of agreement with experimental data indicates that the CBM approach successfully captures the dominant physical mechanisms governing evaporation–condensation interactions. Further refinement of the condensation treatment is expected to reduce residual deviations, particularly under varying thermal loads and transient operating conditions.

In summary, the proposed model demonstrates robust predictive capability for coupled phase-change dynamics in heat pipes and provides a solid foundation for the future development of high-fidelity, physics-based phase-change models in CFD.

5. Conclusions

A comprehensive numerical investigation of a closed two-phase thermosyphon (basic heat pipe) was performed using the Combined Boiling Model (CBM) [26] within a conjugate heat transfer (CHT) framework. The study focused on the preliminary validation of the model against experimental data and previously published numerical results for different heat loads ranging from 173 W to 376 W. The main conclusions can be summarized as follows:

• Numerical stability and initialization:

The two-stage simulation strategy—first solving for steady coolant flow and then initializing the CHT model—provided a stable and physically consistent start-up condition, reducing numerical transients and improving convergence of the coupled phase-change solution.

• Prediction of flow and phase distribution:

The CBM model successfully reproduced the key flow features of the thermosyphon, including vapour generation in the evaporator, vapour-pocket formation and subsequent collapse in the adiabatic section, and thin-film condensation along the upper condenser wall. The predicted flow reversal and liquid-film drainage showed strong agreement with theoretical expectations and experimental observations, confirming that the model accurately reproduces the essential evaporation–condensation mechanisms.

• Validation against experimental data:

Across all investigated heat loads, the CBM model shows close agreement with measurements, with the strongest accuracy in the evaporator and adiabatic regions. Compared with a traditional VOF-Lee formulation, temperature deviations are markedly lower, in some cases by up to an order of magnitude, which highlights the improved predictive performance of the proposed approach. A key practical advantage is efficiency: CBM attains these results with time steps about ten times larger and with coarser meshes than VOF Lee, which leads to substantially shorter runtimes. The two approaches however target different use cases. CBM is well suited for fast, system scale simulations of engineering devices where overall thermal behaviour is the priority, whereas VOF-Lee remains preferable for interface resolved, microscale phenomena such as geyser boiling where an explicitly captured vapour liquid interface is required.

• Influence of fill ratio and heat load:

The simulations accurately captured the sensitivity of thermosyphon operation to the initial liquid inventory. Higher heat inputs led to stronger evaporation, larger vapour volumes, and thicker condensate films, consistent with experimental trends.

• Model limitations and future work:

While the current condensation sub-model remains simplified, the achieved level of accuracy indicates that the CBM framework effectively captures the dominant heat and mass transfer mechanisms. The present work represents only a preliminary validation study and a demonstration of the feasibility of applying the CBM model to heat pipe simulations. Further, more extensive validation under a wider range of operating conditions and with various working fluids is required to extend the model’s applicability and reliability.

These findings highlight the potential of the CBM framework as a reliable and computationally efficient tool for modelling two-phase heat transfer in advanced thermal management systems. This represents the best agreement with experimental data currently reported in the literature for such simulations problems.