Interface Capturing Flow Boiling Simulations in a Compact Heat Exchanger

Forced convection combined with boiling permits effective heat removal from surfaces. High heat transfer coefficients are vital for flow boiling phenomenon especially in nuclear reactors but have relevant applications in steam generators, microchannels, compact heat exchangers, etc. In the latter examples, different fin shapes (rectangular, triangular, wavy, offset strip fin (OSF), etc.) can be added to increase heat removal from the surface. In addition, refrigerant coolant (i.e., R113) as a working fluid has also been used in conjunction with fins to increase heat transfer rate, especially in compact heat exchanges for low temperature/pressure applications [1,2]. However, accurately predicting the two-phase heat transfer coefficient of novel configurations can be challenging.

From an experimental perspective, several studies have been conducted to predict two-phase heat transfer coefficient in mini channels with OSF. In Ref. [3], R113 behavior is examined in a vertical channel with OSF over a wide range of Reynolds numbers and applied heat fluxes. The authors proposed a correlation to find the two-phase forced convective part of the local two-phase heat transfer coefficient. The experimental results showed that the heat transfer coefficient values may lie within ±25% of their predictions. More recently, in Ref. [4], an experimental study was performed with R134a in a brazed heat exchanger with vapor quality up to 0.7. A correlation for Reynonds number factor $F$ was proposed (used for two-phase forced convective heat transfer coefficient evaluation) as well as two-phase frictional multiplier $ϕf$ (applied in two-phase frictional pressure drop correlations). In Ref. [5], the boiling heat transfer characteristerics of HFE-7100 were investigated in a vertical channel with OSF and without them. The results showed that for low heat loads (less than 80 kW/m²), convective boiling was the dominant process in two-phase heat transfer for which OSF may have had significant impact. For higher heat loads (more than 120 kW/m²), the influence of inner channel geometry is negligible as nucleate boiling governs the two-phase heat transfer. In Ref. [6], the authors conducted an experimental study with serrated fins in a horizontal rectangular channel using R134a as the working fluid. They varied heat flux and captured flow patterns from dispersed bubbly to annular flow followed by intermittent dry-out. The results showed that for low vapor quality (less than 0.6) and under high heat fluxes, the mass flux influence on heat transfer coefficient is almost negligible. Conversely, for high vapor quality flow, the heat transfer coefficient increases with mass flux growth.

Due to the high computational costs, numerical investigations are mostly conducted to determine the heat transfer coefficient for single-phase flows. In Ref. [7], the SST $k−ω$ turbulence model is used to predict Colburn factor $j$ for a 3D geometry with OSF. The obtained results lie within ±20% range comparing to a well-known correlation proposed by Manglik and Bergles [8]. Later, in Ref. [9], the researchers proposed a new correlation for $j$ that includes the Prandtl number. This was required as the existing correlations were limited to a working fluid with a low $Pr$ (e.g., for air $Pr$= 0.7). In Ref. [10], it was discovered that the laminar flow heat transfer coefficient is mainly influenced by the growth of the boundary layer. In turbulent flow, vortex shedding also becomes meaningful. 3D simulations using Fluent were performed in Ref. [11] in which the standard $k−ϵ$ model was used for “far from wall” calculations, and the SST $k−ω$ turbulence model was applied for “near wall” calculations. The authors discovered a noticeable difference in results when compared to published work in Ref. [8]. Specifically, differences (∼15%) were discovered for the laminar regime (⁠$Re$ < 120) and a significant difference (∼35%) in the transitional (500 < $Re$ < 1250) and turbulent (⁠$Re$ > 1250) regimes. In both Refs. [12] and [13], 3D simulations in Fluent were conducted and compared with the existing correlation [8]. More geometry parameters were analyzed in Ref. [12], a wider range of $Re$ was investigated in Ref. [13].

To model boiling flows with interface capturing, researchers employ various interface-tracking methods (ITM). For example, a volume of fluid (VOF) based flow solver was proposed in Ref. [14] that reduces computational cost by keeping interfacial temperature at saturation and explicitly coupling it with source terms. This solver was validated, and its capability was demonstrated on 3D flow boiling simulation in a rectangular microchannel. In work done by Ref. [15], a new ITM is used that couples both VOF and level set methods (VOSET). Researchers implemented a microlayer model to take into consideration the evaporation in the microlayer underneath a bubble and conjugate heat transfer between the wall and fluid region. They investigated subcooled boiling of water in a vertical rectangular channel in 2D. To ensure multiple bubble nucleation, 15 nucleation cavities are introduced in the domain. Varying heat fluxes from 100 to 500 kW/m² allowed them to observe flow regime transition from bubbly to elongated bubbly flow. In Ref. [16], the influence of surface characteristics on bubble dynamics is investigated in a vertical rectangular channel. Bubble departure from a solid surface is modeled in 2D. The authors concluded that a hydrophilic surface increased the heat transfer from the wall due to the evaporation of the microlayer. On amphiphilic surfaces, the waiting period (between the departure of one bubble and nucleation of a subsequent bubble) is essentially eliminated which results in an effective contact line heat transfer. Numerical simulations conducted alongside flow boiling experiments are presented in Ref. [17]. In that study, FC-72 is used as the working fluid in a vertical rectangular channel in 2D. The presented computational results showed good agreement with experimental data (including flow visualizations, cross section averaged void fraction, velocity profiles, wall temperature, and heat transfer coefficient). This prediction diverges from known solutions when the mass velocity increases, in part due to 2D modeling limitations. In Ref. [18], the authors investigated flow boiling of FC-72 in a vertical rectangular microchannel in 3D. The results revealed good agreement of numerical predictions for wall temperature variations and overall flow patterns with experimental results. This study also obtained fluid temperature distributions that cannot be measured in small channel experiments.

From the literature review it can be seen that due to high complexity of the flow boiling phenomenon, the majority of high-resolution simulations are conducted in 2D. Only a few cases are conducted in 3D, which impedes comprehensive understanding. In addition, of the studies published, most consider simple rectangular channel geometries which are not necessarily realistic. Thus, we present high fidelity 3D simulations to provide a deeper understanding of the boiling in a vertical mini channels with OSF. We use finite-element code PHASTA to solve Navier–Stokes equations (Sec. 2). The code performance is validated against experiments for a single nucleation site in a vertical rectangular channel [19,20] (Sec. 3). Afterwards, studies for a bubble departing from the wall in a mini channel with OSF are also carried out (Sec. 4). Influence of surface characteristics such as contact angle and liquid superheat on bubble dynamics are analyzed. The results show that there is a small difference in bubble departure diameters (∼4%) for contact angles 30 deg and 45 deg. Local two-phase heat transfer coefficient is also investigated. To achieve higher void fractions, two nucleation cavities are introduced in the same channel with OSF (Sec. 5). For these studies, liquid superheat influence on bubble dynamics is also evaluated. Heat transfer enhancement due to bubble presence is quantified on planes along the bubbles' path.

where $ui$ is the fluid velocity, $ρ$ is the density, $p$ is the static pressure, $τij=μ(∂ui∂xj+∂uj∂xi)$ is the viscous stress tensor, $fi$ is the body force, $cp$ is the specific heat at the constant pressure, $T$ is the temperature, $k$ is the thermal conductivity of a fluid, and $q$ is the dissipation function. For surface tension, the continuum surface force model is used introduced by Ref. [21]. Discretization schemes that are implemented in PHASTA are described in more details in Refs. [22,23].

The Navier–Stokes equations are solved using the “one-fluid” approach for which fluid properties (density, viscosity, specific heat, and thermal conductivity) are found using $Hε$ (e.g., $ρ(φ)=ρlHε(φ)+ρg(1−Hε(φ)) ρl$ and $ρg$ are the liquid and vapor phase densities respectively).

In the above expression, $F1$ is the parameter that changes $FCA$ with respect to the desired contact angle value [26], $TCA$ and $HCA$ are the thickness and height of $FCA$ application region, respectively, $and dwall$ is the distance from the wall. Current contact angle is calculated using interface normal vector (which is found using distance function as $∇φ$⁠) and wall normal vector in wall-adjacent elements. More details on contact angle algorithm implementation in PHASTA and its verification could be found in Ref. [25].

Here, $q¯$ is the volume-averaged heat flux, and $hfg$ is the latent heat of evaporation. This algorithm was verified against analytical bubble growth model [29] and experimental data for pool and flow boiling [30].

In order to validate the performance of the boiling and contact angle algorithms for flow boiling conditions, a separate study is conducted. Experiments in a vertical rectangular channel with HFE-301 were chosen as the references [19,20]. These experiments were conducted for the single nucleation site and such parameters as bubble departure diameters and departure frequencies are available for comparison.

The experiments [19,20] were carried out in a wide range of parameters. However, to match the desired flow conditions, the following flow characteristics are considered: mass flux 140 kg/m²·s (⁠$Re$ = 3125), wall heat flux 9.7 kW/m², inlet subcooling 13.5 K, and atmospheric pressure.

The simulation domain constructed for this study is shown in Fig. 1. The channel dimensions are 324 × 10.1 × 10 mm. It has a total of 264 mm before the test section (including 50 mm unheated and 214 mm heated parts of the channel). The test section is 10 mm long and has a nucleation site of a trimmed conic shape placed in the middle of the section. An additional 50 mm are added after the test section for outlet treatment (to avoid backflow issue and dissipate bubbles since mesh is comparatively coarse there). The heated area in the simulations matches the experimental one.

The real cavity size is of the order of nm which is computationally expensive to model. Instead, three bigger cavity sizes are considered: 0.288, 0.18, and 0.1125 mm of the cavity diameter (the depth is the same for all three cases 0.128 mm) (see Fig. 2). The mesh used near cavities will be reused for OSF simulations (discussed in Secs. 4 and 5).

Exact surface characteristics are unknown [19,20]. However, using information from measured contact angles for refrigerant with similar properties [31] as well as magnified views of high-resolution camera [20], a contact angle 15 deg is set for these validation studies.

where $fi$ bubble departure frequency for $ith$bubble, $tdi$ is the time of departure of $ith$bubble, and $tni$ is the time of nucleation of $ith$bubble.

It was observed that as the cavity diameter decreases, the departure diameters of bubbles also become smaller. This happens as the cavity is always filled with vapor and thus, the bubble interface is attached to the cavity's edges (constant boundary conditions are set for level set function on cavity walls). These boundary conditions are needed to ensure multiple bubble nucleation. Otherwise, there would be no vapor phase left in the cavity. Due to this, a smaller cavity diameter leads to faster bubble departure.

The results for the smallest cavity agree well with experimental data [20]. From bubble departure frequency plot (Fig. 3, top right) it can be seen that the data is consistent (i.e., smaller bubbles depart faster from the cavity). It is also noted that the modeled values are somewhat lower than observed experimentally (e.g., $MSE$ for the smallest cavity is 146.2 Hz).

where V is bubble volume $(πD36)$ in cubic meters and $tb$ is the nucleation period in seconds. The results are shown in the bottom of Fig. 3. From this plot it could be seen that the amount of steam is consistently lower for the smallest cavity (as bubble departure frequency results discussed before). This could be attributed to microlayer evaporation that is not considered in this study. In future work, additional source terms will be considered to account for the amount of steam generated in this tiny liquid film close to the wall.

To model the flow through OSF, a part of a vertical rectangular channel used in experiments [3] is considered (Fig. 4, left). In the domain, only the flow region is modeled. Therefore, there is a “step” of 0.2 mm corresponding to the fin thickness (Fig. 4, right) that is present in the $y$-direction. The flow is resolved for the first 12 mm of the domain, other parts are created to increase the code robustness (coarse meshes are used to prevent backflow and dissipate bubbles, the outlet treatment box includes last 6 mm (not shown in Fig. 4, right)).

A no-slip boundary condition is imposed at the top and bottom of the domain where the flow is periodic in the $z$-direction. Inflow/outflow boundary conditions are set for velocity and temperature profiles in $x$-direction. Gravity is applied in the opposite direction of the flow. There is an uncertainty in the velocity and temperature profiles in the channel. Therefore, it is assumed that the flow is fully developed in the modeled part of the channel. Also, based on experimental data from Ref. [3], $Re$ = 200 is chosen for these simulations. Parabolic velocity and temperature profiles are considered to be representative for this flow. Constant temperature $Tw$ = 323.3 K (2 deg superheat) is set at the right and left walls of the domain (which represent inner channel walls assuming highly conductive fin material). According to FLUENT simulations [32], this value corresponds to heat flux of 1000 W/m², which was used in the experiment in Ref. [3]. Since there is an uncertainty in exact temperature values at the inner walls, superheat studies were performed. Fluid properties for R113 at 102.4 kPa are presented in Table 1. Initial time-step is set 5 × 10⁻⁶ s and was automatically adjusted during the simulation to keep Courant–Friedrichs–Lewy $(CFL)$ number = 1.0.

In such a complex geometry, considerable friction losses are present when the flow passes through the OSF. To properly resolve these losses, boundary layer meshing structure (BLs) is created on wall surfaces (Fig. 5). To accurately represent a bubble (around 25–30 elements across the bubble diameter are needed [33]) as well as to minimize the computational cost, several mesh refinement regions are introduced (Table 2). The bubble with an initial radius of 0.3 mm is placed at the location (0.00825; 0.000135; and 0.004775) m. This size was chosen to save computational resources (as smaller bubble would require higher mesh resolutions). The bubble is located at the bottom step of the fin. Due to stagnant liquid at the step corners, more superheated liquid is accumulated here. Therefore, bubble departure diameter is expected to be bigger rather than when the bubble is initialized at the top of the fin. In addition, the target contact angle of 30 deg is set. This value was also a source of uncertainty (since exact surface characteristics are unknown), therefore, a contact angle study was performed. For the contact angle algorithm to perform as designed, BLs are needed on the wall (Fig. 5, bottom).

It was mentioned that parabolic velocity and temperature profiles are imposed at the channel inlet to achieve laminar (⁠$Re$ = 200) flow. However, due to the staggered geometry of the channel (see Fig. 5, top), initially these profiles are somewhat different near the right and left walls. To overcome this discrepancy, single-phase simulations are conducted first to make sure that the temperature and velocity gradients are consistent near the walls. The time step size was adjusted automatically keeping maximum $CFL$ number 20.0. The resulted velocity and temperature profiles are shown in Fig. 6. Single-phase simulations are carried out for all three bulk superheats considered (0, 2, and 4 deg). Then, the solutions at the last time step are set as initial conditions (ICs) for flow boiling simulations.

Contact angle is an important parameter for boiling simulations. To investigate the influence on bubble growth and departure, two contact angles are considered: 30 deg and 45 deg. Bubbles at the initial time step are presented in Fig. 7, top. The resulting bubble characteristics are shown in Table 3. Using the bubble tracking algorithm [28], bubble growth statistics are obtained (see Supplemental Figure S2 available in Supplemental Materials on the ASME Digital Collection).

From Fig. 7, one may see that the bubble with a higher contact angle reaches smaller departure diameter and departs faster. Overall, the difference in bubble departure diameters for two contact angles considered is ∼4%. The results are obtained for the single bubble departure event only. To verify this observation, multiple bubble nucleation events are needed (to have more reliable bubble statistics). For that purpose, Sec. 5 discusses flow boiling simulations conducted in the same geometry with OSF, but with two nucleation sites.

Parabolic temperature profile as well as wall temperature are the sources of uncertainty in the current study. To investigate the influence of wall temperature values on bubble growth, three superheat studies are conducted (Table 4). The simulations were conducted until a bubble reached the coarser mesh region (mesh #6, Fig. 5). The average distance from the nucleation site that the bubble covered is around 0.00375 m. Bulk superheat is chosen to keep the same temperature difference between the wall and the channel center. Mesh design and other model conditions are the same as in the contact angle study. For the current simulation, the target contact angle is 30 deg, advancing contact angle is 25 deg and receding contact angle is 35 deg. The resulted bubble growth statistics are shown in Supplemental Figure S3 available in Supplemental Materials, and bubble positions at the departure moments are presented in Supplemental Figure S4 available in Supplemental Materials on the ASME Digital Collection.

Higher superheat values at the surface lead to a more intensive evaporation process. Because of that, bubbles reach bigger sizes for cases #4–5 compared to case #3. However, since vapor prevents liquid from coming closer to the surface, it exacerbates the heat transfer for case #4–5 relative to case #3. On the other hand, the small-size bubble (as in case #3) departs much sooner from the wall and, thereby, heat transfer coefficient increases. Another observation is that a larger bubble (cases #4–5) slides longer along the wall (see bubble departure locations in Supplemental Figure S4 available in Supplemental Materials on the ASME Digital Collection: in case #4 the bubble slides to the fin corner and then departs; in case #5 the bubble departs only in the middle of the next fin since it was sliding along the surface before that). Such bubble behavior tends to disturb thermal boundary layer, hence, increasing heat transfer coefficient [34,35]. More details on heat transfer coefficient evaluation will be presented in Sec. 5.3.3.

Since only one bubble is present in each simulation, a smaller part of the channel is considered (Fig. 8, left) for local void fraction evaluation. The resulting void fraction statistics for each case are shown in Fig. 8, right. It is noted that since for case #5 the bubble size is larger, its corresponding void fraction has higher values compared to cases #4 and 3. The bubble growth is somewhat similar for each case: as the bubble attaches to the surface, its volume grows nonlinearly. However, when the bubble departs and starts to move toward the channel center, the growth rate becomes more linear. Since for case #3, there is no superheat at the channel center, bubble growth rate in this region is considerably smaller than, for example, in case #4. For case #5, there is 4 deg superheat at the channel center which, consequently, leads to a much higher bubble growth rate compared to cases #3 and 4. The final local void fraction achieved for cases #3–5 are 0.57%, 1.42%, and 4.78%, respectively.

Bulk superheat studies show that the temperature distribution in the channel plays a crucial role in heat transfer evaluation. The difference in 2 deg (cases #3 and 4) results in ∼50% increase of the bubble departure diameter and the change in 4 deg (cases #3 and 5) — in nearly 170% growth.

In the equation above, values are found from Fourier's Law using area-averaged normal to the wall temperature gradient values (⁠$1S∫S∇TdS$ in Eq. (11)), thermal conductivity of R113 liquid phase $k$ (see Table 1 for reference), and known temperature differences at the walls and the centerline (⁠$ΔT$ = 2 deg). More details on that could be found in Ref. [36]. These temperature gradients were obtained at 0.04 mm from the wall. This value is chosen to make sure that temperature changes linearly over this distance. The calculations are performed for each plane on the bubble path (see Supplemental Figure S6 available in Supplemental Materials) at four time intervals: 0, 9, 18, and 27 ms. Points at $t=$ 0 s correspond to single-phase heat transfer coefficients.

Bubbles with low density vapor move faster than liquid. Correspondingly, a bubble of a bigger size travels further downstream (see Supplemental Figure S1 available in Supplemental Materials). The results show that local heat transfer coefficient for case #5 is generally higher than for cases #3–4. This could imply that the nucleate boiling part of the two-phase heat transfer coefficient that depends on the liquid superheat plays an important role (for case #5 the superheat rate at the wall is the highest). At the same time, convective boiling heat transfer part is also significant for case #5: from snapshots of bubbles' location, one could notice a larger disturbance of the thermal boundary layer caused by the biggest bubble sliding along the wall (compared to cases #3–4).

By analyzing the change in the heat transfer coefficient on plane #1, it was concluded that as the bubble grows, it acquires more energy due to phase change. Then as the bubble starts to move downstream following the flow velocity, it leaves the plane, and the heat transfer coefficient correspondingly decreases for all cases. When the bubble comes closer to plane #2, the heat transfer coefficient is increased significantly. It should be noted that since for cases #3–4 bubbles traveled smaller distance over this plane, the growth of the heat transfer coefficient is not as significant as for case #5. For plane #3, the heat transfer coefficient starts to increase as the bubble would closer approach this surface. For planes #4 and #5 (side walls), the heat transfer coefficient trend can be explained as follows: when the bubble is relatively far away from the side surface (as, for example, in cases #3 and 4), due to phase change, it consumes the amount of liquid in the vicinity of the wall, thereby, degrading convective heat transfer. However, when the bubble comes closer to the wall (or when its size increases), the boiling heat transfer prevails. Observed local heat transfer coefficient values are obtained for the section of the channel where the bubble is nucleated. Overall, the heat transfer enhancement effect is more pronounced for plane #1 where the bubble is nucleated. Once it starts to grow (at $t=$ 9 ms), the heat transfer coefficient is increased by 44.2% (case #3), 86.1% (case #4), and 315.7% (case #5). To validate the results, the authors intend to compare it with the experiments [3] as the part of the future work. However, for this a higher void fraction should be achieved which imposes additional numerical challenges.

To increase void fraction, two nucleation cavities are introduced in the same channel on the right and left walls (see Supplemental Figure S7 available in Supplemental Materials on the ASME Digital Collection). Cavities have conic shape (diameter 0.288 mm, depth 0.128 mm) and are trimmed at the bottom. The two cavities are located at (8.25; 0; 4.675) and (9.75; 2.8; 4.675) mm. Exact locations of nucleation sites in experiments [3] are unknown as they depend on surface imperfections and may change for different heat flux conditions. Information on mesh used for these studies is available in Supplemental Materials on the ASME Digital Collection.

As for the previous runs, since the exact temperature values at the inner channel walls likely depend on the streamwise location in the full heat exchanger design (in experiments the size of the channel is 750 × 3 × 100 mm), several simulations were conducted with three different bulk superheat values: 0, 2, and 4 deg. These simulations are conducted with $Re$ = 200. For these simulations the target contact angle is set to 15 deg, advancing contact angle is 10 deg, and receding contact angle is 20 deg. The same working fluid is considered R113 (Table 1). Simulations are carried out with the initial time step 2.1 ms, and it was adjusted during the runs to keep flow $CFL$ number = 0.5.

It has been observed that when a bubble grows at the wall, there is a thin liquid film exists between the bubble and the wall called microlayer [37,38]. The contribution of the microlayer to bubble evaporation is substantial (heat fluxes reach 1 MW/m² [39]) especially at atmospheric pressure. However, as was shown in Ref. [40], heat flux from liquid evaporation in the microlayer region is negligible for low heat fluxes. Since in the current simulations, applied wall temperature corresponds to heat flux 1000 W/m², the contribution from the microlayer is not considered.

Bulk superheat plays a major role in the bubble growth, and it is important to accurately represent it in the simulations. There are several examples in which conjugate heat transfer between the solid and fluid domains was accounted for in boiling simulations, including [41–44]. However, in the current studies, wall thickness was small (0.2 mm), and the material is highly conductive (copper) [3], therefore, conjugate heat transfer is neglected.

To quantify the difference in the temperature distribution in the channel, sensitivity studies with three different bulk superheats are conducted. Departure diameters of the first-generation bubbles (e.g., the first bubble produced at each site) are summarized in Table 5. The case with 0 deg superheat is considered as a baseline for this study. From Table 5 it can be seen that a 2 deg variation in the bulk superheat causes a significant difference between departure diameters. This effect is more pronounced for 4 deg superheat difference. In Supplemental Figure S5 available in Supplemental Materials on the ASME Digital Collection, bubble shapes for the first-generation bubbles at the right wall cavity are shown at simulation time 6 ms. It could be seen that for 2 and 4 deg bulk superheat, bubbles readily depart, whereas for the 0 deg bulk superheat case, the bubble is more closely attached to the wall. Compared to the bubble departure from the wall simulations (the results are shown in Table 4), current bubble departure diameters are smaller since a smaller contact angle (i.e., 15 deg) is imposed.

Some of the bubble statistics are shown in Fig. 10, bubble trajectories are available in Supplemental Material. From these plots it can be seen that departure diameter is generally higher for the higher bulk superheat. For 2 and 0 deg bulk superheat cases, the obtained departure diameter values are more consistent whereas for 4 deg a higher discrepancy was observed (which was attributed to coalescence between just departed bubbles and bubbles left in cavities). Bubble departure frequency is more consistent when there is no coalescence between already departed bubbles and newly nucleated ones.

Snapshots of simulations for 4 deg bulk superheat are shown in Fig. 11; for the 0 and 2 deg bulk superheats in Fig. 12. One could notice that for higher temperature gradients values near the walls, bubbles grow larger (for example, compare all three cases at 8 ms). For 4 deg bulk superheat, bubbles reach around 1 mm in diameter (more detailed data is presented below) and they extensively coalesce with each other (see snapshot for 32 ms, Fig. 11). At the end of the simulation, the vapor phase occupies more than a half of the channel cross section which could later result in a flow transition (from bubbly to slug flow).

Results shown in Fig. 12 indicate that for lower bulk superheat values (e.g., 0 deg) bubbles have smaller departure diameters. Due to the size, bubbles coalesce less often and stop growing once they are near the centerline (there is no superheat here). For 2 deg bulk superheat, more bubbles are observed compared to the previous case due to higher temperature values in the channel.

Heat transfer coefficient is calculated locally on planes along bubbles' path (see Fig. 13) for three considered cases. The same procedure was used to calculate local heat transfer coefficients as was discussed previously in Sec. 4.7. The results are shown in Fig. 14. For all three plots presented, points at $t$ = 0 s correspond to single-phase heat transfer (these are solutions at the last time-step of single-phase simulations which are recycled as initial conditions for the two-phase flow runs).

If one considers heat transfer evolution for bulk superheat of 0 deg, it may be seen that for plane #1, the value does not change significantly. This implies that conditions for bubble growth and departure stay consistent over time. The plots from Fig. 10 for nucleation site #1 also support this point. The same conclusion can be made for heat transfer coefficient behavior on plane #1 for bulk superheat 2 deg. For 4 deg bulk superheat, there is a decline in heat transfer after ∼10 ms. This happens due to coalescence between the departed bubble and the newly nucleated one. The newly merged bubble stays longer near the wall consuming more energy from the surrounding liquid. For the subsequent bubble (point at ∼16 ms), no coalescence event was observed. At ∼32 ms, an increase in heat transfer for plane #1 attributes to the fact that the cavity is full of liquid (no vapor phase is left at this time). The situation is different for plane #6: increases in heat transfer coefficient (for bulk superheats 0 and 2 deg after ∼60 and 30 ms correspondingly) suggest that small-size bubbles are nucleated (could be seen in Fig. 10 as well). Thus, more liquid occupies the cavity causing the rise of heat transfer. For 4 deg bulk superheat, there is a coalescence event again (at ∼10 ms) happening between the departed and newly nucleated bubbles which contributes to a decrease in heat transfer coefficient.

For plane #2, a decrease in heat transfer coefficient is observed for all three cases (starting from ∼8 ms). This is the result of the numerical assumption from the boiling algorithm — vapor phase is considered to be at saturation temperature. During single-phase simulations, bubbles initially nucleated in the domain work as heat sinks. As the code runs, this saturated temperature affects thermal boundary layer (areas of lower temperature near walls could be seen at 8 ms in Figs. 11 and 12). Although this effect distorts heat transfer evaluation at the beginning, it disappears as boiling simulations are run longer.

For plane #3 one could clearly see the effect of bulk superheat on heat transfer: before bubbles reach this plane, the average heat transfer coefficient values are 114.7 (0 deg), 189.5 (2 deg), and 253.5 (4 deg) W/(m²·K).

For planes #4 and 5, the situation is similar: values are nearly constant until bubbles come closer to these planes causing disturbances in thermal and hydrodynamic boundary layers (e.g., see the growth on plane #4 after ∼60 ms for bulk superheat 0 deg, Fig. 14).

Interface capturing flow boiling simulations were conducted for a vertical mini channel with offset strip fins (OSF) using code PHASTA. To establish the credibility of the results, model validation is conducted for a single nucleation site in a vertical rectangular channel. The observed bubble departure diameter agrees well with the experiments for the smallest cavity size considered. The steam generation and bubble departure frequency are underpredicted which could be attributed to microlayer evaporation (not considered in the current study).

Second, a complex geometry with OSF under single-bubble flow boiling was analyzed. Bubble dynamics characteristics are investigated with varying contact angle and bulk superheat. It was observed that temperature distribution in the channel plays a major role in determining heat transfer process for this laminar flow (⁠$Re$ = 200). The increase in bulk superheat by 2 deg leads to a nearly 50% growth of the bubble departure diameter and the increase by 4 deg allows to achieve almost a 170% growth.

Third, to observe multiple consequent bubble nucleation, two nucleation cavities were introduced in the same channel with OSF. Single-phase solutions are used as initial conditions for boiling simulations. Bulk superheat sensitivity tests are conducted to analyze bubble dynamics. Analysis of departure diameter values shows that the 4 deg temperature difference in the bulk fluid may result in a 40% difference in bubble departure diameter value. Local heat transfer coefficient is evaluated on each plane along the bubbles' path. It is shown that for some cases (e.g., for bulk superheat 0 deg, plane #3) bubbles approaching the wall may increase heat transfer coefficient by up to 2.5 times. With a smaller nucleation diameter, more liquid is cooling the cavity which, in turn, results in nearly twice as high heat transfer coefficient (e.g., for bulk superheat 2 deg, plane #6).

Future work will be focused on improving model predictions of bubble departure frequency and amount of steam released from a cavity as well as testing PHASTA capabilities on modeling flows with a higher void fraction.

ACUSIM linear algebra solution library was employed which was provided by Altair Engineering Inc. along with mesh and geometry modeling libraries provided by Simmetrix Inc. for simulations. We also acknowledge the computing resources provided on Henry2, a high-performance computing cluster operated by North Carolina State University. Authors greatly appreciate the constructive and useful comments from the reviewer as well as valuable feedback from Edward Chen.

• Current work is performed under Two-Phase Flow DNS Phase 2 Project supported by Mitsubishi Heavy Industries, Ltd.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

$cp$ =

specific heat at the constant pressure, J/(kg·K)

$dwall$ =

distance from the wall, m

$fi$ =

body force (including gravity and surface tension), Pa/m

$FCA$ =

contact angle force, H

$F1$ =

contact angle algorithm parameter

$hfg$ =

latent heat of evaporation, J/kg

$HCA$ =

height of contact angle force application region, m

$Hε$ =

smoothed Heaviside function

$k$ =

thermal conductivity of a fluid, W/(m·K)

$p$ =

static pressure, Pa

$q$ =

dissipation function, J/(m³·s)

$q¯$ =

volume-averaged heat flux, W/m⁴

$R0$ =

bubble radius, m

$R1$ =

radius of temperature gradient collection outer shell, m

$t$ =

time, s

$tb$ =

nucleation period, s

$T$ =

temperature, K

$TCA$ =

thickness of contact angle force application region, m

$ui$ =

fluid velocity, m/s

$V$ =

bubble volume, m³

$ρ$ =

density, kg/m³

$τij$ =

viscous stress tensor, Pa

$ε$ =

interface thickness, m

$ρl$ =

liquid density, kg/m³

$ρg$ =

vapor density, kg/m³