Numerical study of boiling of Liquid Nitrogen on a liquid-liquid contact plane
Numerical study of boiling of Liquid Nitrogen on a liquid-liquid contact plane
Rupak Kumar, Arup Kumar Das Department of Mechanical and industrial engineering, Indian Institute of Technology, Roorkee, 247667, India
Abstract
In this paper, a numerical study is conducted to investigate boiling of a cryogen on a solid surface as well as on a liquid surface. Both single mode and multi-mode boiling is reported for boiling on to a solid surface. In case of boiling on a liquid surface, liquid nitrogen (LN ) is selected as the cryogen (boiling fluid) and water is chosen as the base fluid (heating fluid). Different flow instabilities and their underlying consequences during boiling of a cryogen are also discussed. For the boiling on a solid surface, in the single mode, bubble growth, its departure, and area weighted average heat flux are reported, where they increase linearly with increase in the wall superheat. Asymmetry in the bubble growth and departure of 2 nd batch of the vapor bubbles have been observed due to local fluctuations and turbulence created just after the pinch off of the 1 st batch of vapor bubbles in case of multi-mode boiling on the solid surface. Boiling of LN on a liquid surface is reported for a base fluid (Water) temperature of 300 K. Vapor film thickness decreases with time and the minimum film thickness just before rupture is 7.62 ππ , dominance of thermocapillary over vapor thrust causes breaking of the vapor film at t = 0.0325 s. The difference in evaporation rate and vapor generation, before and after vapor film collapse is significant. Introduction
Boiling of a liquid (boiling fluid) in direct contact with another liquid (heating fluid) can be found in many industrial applications as well as in several natural incidents. One such scenario occurs during accidental spills of liquefied natural gases (LNG) while transportation via cargos in sea water. LNG, which is composed mostly of methane and other higher alkanes such as ethane, propane, butane and nitrogen, when comes in direct contact with water may transform its phase rapidly depending upon the local flow conditions and several other factors. This rapid transformation of phase also known as rapid phase transition (RPT) generally occurs whenever there is absence of any imperfections on the surface, it can be a mirror finished surface or a liquid-liquid contact plane. Boiling fluid gets superheated beyond its saturation temperature ( T S ), causing it to change its phase homogeneously at an instant. Rapid phase transfer and thereafter change of state, causes gas to expand quickly which resembles like a physical explosion. This violent explosion may initiate shock waves and can harm equipment and people in the vicinity, also, if the gas is flammable and comes in contact with a fire source, it can cause some serious damage. Therefore, there is a need to fully understand this complex phenomena of liquid-liquid boiling and the critical factors which causes this sudden change of phase and methods to control this behavior. Superheat limit of a fluid is the maximum attainable temperature at a given pressure. Liquids upon reaching this superheat limit changes its phase from liquid to vapor phase instantly which also expands simultaneously just after the phase change to ambient pressure which looks like a physical xplosion. Hence, in the past several attempts by different authors were made to calculate the superheat limit of the various fluids. Blander and Hengstenberg (1971) calculated the homogeneous nucleation temperature for different fluids and mixtures using steady state Zeldovich-Kagan theory and verified the results experimentally. They found that the theory predicts accurately the superheat limit for pure liquids but insisted on requirement of further investigation for binary mixtures. Rausch and Levine (1973) developed a thermodynamic and heat transfer model to predict the superheat limit of the cryogen to cause a shock wave, which they found out to be 0.84 times the critical temperature ( T C ) of the cryogen. Based on which they also provided an expression to predict the necessary temperature of the base fluid to superheat the cryogen up to its superheat limit ( T SL ) and eventually cause a shock wave. They verified the expressions with the experimental results of different fluids. They concluded that if the cryogen is a zeotropic mixture, the selective boil-off of the lighter hydrocarbon will make the mixture rich in the heavier hydrocarbons, which in turn will increase the mixture critical temperature, and hence the superheat limit and thus shock will be strongly dependent upon the composition of the mixture. Reid (1976) experimentally investigated superheating of a liquid droplet (n-pentane) using a bubble column filled with sulfuric acid and wrapped with a heating element, wrapped closely at the top compared to the bottom so that the top is warmer compared to the bottom, and a liquid drop is injected at the bottom which explodes the moment it reaches the warmer region where the local temperature is equal to its superheat limit. In this study, based on the kinetic analysis, he found that the liquids must be superheated in order to form vapor bubbles in the bulk fluid, since the small vapor bubbles are quite unstable and collapses readily. Nishigaki and Saji (1981) verified an empirical law with the experimental results for the superheat limit of oxygen and nitrogen. Sinha et al. (1987) experimentally determined the superheat limit of the liquid nitrogen using a transient heating method for different temperature and pressure range. After understanding the superheating, the primary reason behind the rapid phase transition (RPT), and its quantification, authors have focused on physical and thermodynamic understanding of this complex behavior of RPT. Enger et al. (1973) proposed a physical model describing the mechanism of boiling of liquid-liquid systems like liquefied gas-water systems, where they claimed that the lack of nucleation sites shifts the transition region towards higher temperature, superheating the boiling fluid. They concluded that the superheating of the boiling fluid depends upon the temperature difference between the fluids and the composition of LNG, also, vapor explosion would not occur in case the methane content of LNG is less than 40 mol %. They also claimed that the energy released upon vapor explosion of liquefied-gas water systems is limited by the superheat limit of the fluid. Drake et al. (1975) conducted and experimental study to compare the boil-off rates of different cryogenic fluids, liquid nitrogen (LN ), methane and ethane upon boiling on a water surface. They found that all the three cryogenic fluids behaved differently upon spillage onto water surface, LN boiled in a film boiling regime with minimum heat flux but maximum vapor superheat, ethane boiled in nucleate boiling regime with maximum heat flux and minimum superheat, methane on the hand found to be the intermediate case, where both nucleate and film boiling regime occur, heat flux and vapor superheat ranges between the data obtained for LN and ethane. apor film collapse between the base fluid and the cryogenic fluid leads to superheating of the cryogenic fluid and later shifting of boiling regime from film boiling to transition and finally nucleate boiling. Yalencia-Chavez (1973) studied the confined spillage of LNG on water surface and the effect of its composition on the evaporation rates. Experiments were performed to evaluate evaporation rates for pure methane, ethane, propane, their mixtures and LNG. He observed that methane initially boils in the film boiling regime, until the formation of an ice layer which promotes nucleate boiling, ethane boils in transition boiling regime and again formation of ice shifts the regime to nucleate boiling, but propane boils in nucleate boiling regime since beginning, whereas LNG mixtures evaporates in a preferential manner, where, more volatile component evaporates first, this causes rise in the saturation temperature of the residual liquid. LNG mixture starts to boil with film boiling regime but preferential evaporation causes methane to evaporate quickly, decreasing the vapor pressure in the vicinity of base of the bubble, which causes vapor film collapse. He concluded that, with the increase in amount of higher hydrocarbons in the LNG, probability of vapor film collapse also increases. Dincer at al. (1976) experimentally investigated the effect of initial temperature of water on evaporation rate of liquid nitrogen (LN ) and methane on water. They concluded that if the initial temperature of the water is low, the energy transfer from water surface to the cryogenic fluid is provided by latent heat of formation of ice at the surface with little effect on underlying water, whereas at higher temperature of water surface, energy is supplied through convection with homogeneous cooling of water, formation of ice surface also leads to shift in boiling mechanism from film to nucleate boiling, which was also observed by Yalencia-Chavez (1973). Vesovic (2006) studied the effect of ice formation on the evaporation of LNG over water body. He also concluded from this study that formation of ice leads to sharp change in temperature of water surface leading to vapor film collapse and change of mechanism from film to nucleate boiling. Vapor dispersion after the rapid phase transition is major cause of concern, hence, several studies have been conducted to predict the vapor dispersion, and various ways to control the dispersion. Benjamin et al. (2008) identified the key parameters which affect the vapor dispersion and used a CFD model to find their effect on the vapor dispersion upon spillage of LNG. Studies were conducted for two cases, in the first case LNG was released on a water surface and in the other case on a concrete surface. Some of the key parameters which had high impact on the vapor dispersion for the first case i.e. LNG released on water body, and in the presence of high as well as low wind velocity includes release rate, wind velocity and sensible heat flux from the ground, whereas some key parameters had high impact only at low wind velocity such as pool area, velocity profile of the gas phase and the turbulence at the source, while other parameters such as pool shape and geometry, ambient temperature, surface roughness and wind direction have relatively low impact as compared to other parameters discussed before. A similar study was also performed by Qi et al. (2010) where they also discussed the key parameters such as evaporation rate, pool area, atmospheric conditions, turbulence in the source term, ground surface temperature, height of the roughness and obstacles. They compared the numerical result with the experimental data and found to be in good agreement. Recently, Horvat (2018) developed a CFD based model to simulate spillage of LNG on water surface, its spreading, rapid phase transition and finally vapor dispersion. They discussed the effect of temperature of vapor cloud on its dispersion, initially just after the phase change, the temperature of vapor cloud at the saturation temperature i.e. -162 β, and hence t is negatively buoyant and vapor clouds floats just above the water surface, but after mixing with the surrounding air, when its temperature goes above -108 β, it become positively buoyant which lifts the vapor cloud from the water surface and helps in its dispersion. They observed that the effect of RPT event on vapor dispersion is limited but RPT has significant effect on the shock wave formed due to sudden expansion of vapor cloud. Some studied have also been conducted to contain the vapor cloud dispersion. Rana et al. (2010) conducted field experiments to understand the vapor cloud dispersion and ways to mitigate the dispersion upon LNG spillage on the land. They used two different techniques, water spray curtains and flat fans to control the vapor cloud dispersion, where they found water spray curtains provide insignificant dilution of vapor cloud. However, the flat fans Experimental studies have been conducted by numerous authors to model the source term in order to describe this complex phenomena of cryogenic boiling through numerical study. One such study was conducted by Vechot et al. (2013) where they investigated the effect of different heat transfer mechanism i.e. conduction, convection and radiation, and their individual contribution on the vaporization of cryogenic fluid, liquid nitrogen (LN ). They used a Dewar flask to to store the LN also minimize heat transfer through conduction and an electronic weighing machine to measure the vaporization rate. They conducted two sets of experiments to measure the vaporization rate of LN , in the first set, they divided the experiment into four groups, where in the first group they checked the effect of only conduction through walls and the lid at the top, then in the second group, along with the conduction, effect of convection in the room and natural radiation, then in the third group, effect of radiation from a bulb (200 W), convection in the room along with the conduction heat loss from all sides, and finally effect of forced convection where they used an electric fan, along with the natural radiation in the room and conduction through walls, in the second set of experiments, effect of forced convection are investigated, where they varied the wind speed from 1.9 m/s to 3.3 m/s, and also considered the conduction from all the sides. They observed that the vaporization rate increased significantly in case of forced convection and ration heat transfer, but contribution from convective heat transfer was significantly higher compared to other heat transfer mechanisms. Gopalaswami et al. (2015) experimentally investigated the vaporization heat flux of liquid nitrogen (LN ) continuously released onto water surface, where they used an electronic balance to measure the mass flux of vapor based onto the mass loss data for different spill rates. They predicted the heat flux using the energy balance for the nitrogen-water system. They also developed a boiling regime map for liquid nitrogen boiling on water using the experimental results, they compared the experimental results with the predicted results based on the expression developed by Berenson (1961) for film boiling heat transfer on a horizontal surface, which they found to be predominantly in the film boiling regime. They found the vaporization flux to be independent of amount of LN spilled, however the initial flux was dependent on the water temperature and spill rate. They also observed that the amount of vapor generation increases with spill area and the size of water body. Morse and Kytomaa investigated the effect of turbulence on evaporation rate of liquid nitrogen (LN ) and LNG. They calculated the turbulence intensity by measuring the temporal variation of the interfacial height at LN -water interface. The velocity magnitude variations were measured by comparing the potential and kinetic energy at the interface. They found that the rate of evaporation ncreases considerably with turbulence intensity for the spillage of liquid nitrogen (LN ) on water as compared to that of LNG. Stutz et al. (2013) experimentally investigated the role of evaporation waves in the propagation of boiling phenomenon in a superheated n-pentane on a copper block. They have used a mirror finished surface to reach the superheat limit. Thy concluded that in the case of a very high degree of liquid superheat, boiling on a mirror finished surface doesnβt follow the classical bubble growth geometry, in its place, a single bubble is seen first followed by phase change around that bubble which spreads around that bubble in all directions, which was referred as βStraw hat structureβ. Hence, they proposed a physical model based on the theory of evaporation waves which explains the phase change phenomena for this case. The proposed model divided the whole phenomena into three different physical phenomena, where at first the vaporization front was visualized as the evaporation wave, then the superheated liquid was partially vaporized in the vaporization front while remaining form small droplets and impinge on the surface which depending upon the heat flux go through partial phase change, and at last the development of straw hat structure drives the superheated liquid and changes the flow field nearby. Very few authors have worked on the numerical study of boiling of a cryogenic fluid, one such study was conducted by Liu et al. (2015) where they developed a model for cryogenic boiling on a solid surface, model includes both heterogeneous and homogeneous mode of boiling. The model successfully predicted different boiling regimes, nucleate, transition and film boiling for different wall superheats. They developed a boiling regime map for liquid nitrogen (LN ). They also investigated the effect of liquid level on the heat flux, which decreases with the increase in the liquid level. In a similar study, Ahammad et al. (2016 a&b) developed a CFD model for film boiling of a cryogenic fluid on a solid surface. They discussed the effect of wall superheat onto the average heat flux, bubble frequency and bubble diameter, they also proposed model based on the energy balance to predict the heat flux in the case of film boiling. Both the studies discussed were based on Rayleigh-Taylor instability (RTI), and are discussed in detail in the results and discussions sections. Considering the literature review, a very limited studies have been conducted for boiling of a cryogen on the solid surface, moreover all the studies have been performed for single-mode boiling assuming steady and uniform behavior throughout the domain. A multi-mode boiling for the cryogenic fluid on a solid surface has never been conducted before, moreover to the best of knowledge of the authors, interfacial study of boiling of a cryogenic fluid on a liquid surface has also not been performed before. Hence, in the present study, boiling of a cryogenic fluid (LN ) has been investigated on a solid (single-mode and multi-mode) as well as on a liquid surface.
2. Mathematical model
Following section describes the computational domain, mesh structure, governing equations, phase change model and surface tension model utilized for this study.
A computational domain of width of π π/2 and a length of π/2 (Figure 1) is chosen for the study of boiling of LN on solid surface, where π π/2 is the most dangerous wavelength (Carey, 1992), given by equation (1). π = 2π [ 3π(π π β π π£ )π] (1) Where, π is the interfacial tension between liquid nitrogen (LN ) an nitrogen vapor (N ), and π π , π π£ are density of LN and N respectively. A uniform structured mesh of 64 x 192 segments is selected for the present study based on grid independence test. Pressure outlet boundary condition is set up at the top, while symmetry boundary condition on either adjacent walls, at the bottom wall, no-slip boundary condition is set. Both liquid nitrogen (LN ) and nitrogen vapor (N ) are initially set at 77 K. The initial vapor film (Figure 1) is patched as a sinusoidal disturbance based on the Rayleigh-Taylor instability and given by equation (2), also. A linear temperature gradient is also applied in the initial vapor film (Ahammad et al., 2016) and given by equation (3). πΏ = π π
64 (4 + cos (2ππ₯π π )) (2) π π¦ = {π π β π₯π β π¦πΏ , πΌ = 1π π ππ‘ , πΌ = 0 (3) Figure 1.
Schematic of computational domain along with the boundary conditions
The equation for mass conservation (4), momemtum conservation (5), energy conservation (8), phase fraction transport equation (9), phase change model (13,14 & 15) and surface tension model (16) are as follows: LN (77 K)N T = T S T W S y mm e t r y S y mm e t r y Pressure outlet Ξ» d /2 Ξ» d / ( T y ) No-slip , Continuity equation ππππ‘ + π΅ β (π π) = π π (4) Continuity equation is by equation (4), where, π and π are density and velocity vector respectively, whereas π π is the source term to account for the mass transfer phase during boiling. Momentum equation πππ‘ (π π) + π΅ β (π ππ) = βπ΅π + π΅ β {π(π΅ π£ + π΅π£ π )} + ππ + π πΌπ½ (5) πΉ ππ = β π πΎπΏ πΌ πΎ π πΎ π πΏ π΅πΌ πΏ + πΌ πΏ π πΏ π πΎ π΅πΌ πΎ
12 (π πΎ + π πΏ ) π<π (6) π πΎ = ππΌ πΎ |ππΌ πΎ | , π πΏ = ππΌ πΏ |ππΌ πΏ | (7) Momentum equation is given by equation (5) where, π , π , π are pressure field, dynamic viscosity and gravitational constant respectively, whereas πΉ π π‘ is the surface tension force given by equation (6). It is based on the continuum surface force model, given by Brackbill et al. (1992). In equation (6), π ππ , π πΎ and π πΏ are interfacial surface tension between phase i and phase j , curvature of phase i and j respectively. The curvatures of phase i and phase j are given by equation (7). Energy equation π(ππΈ)ππ‘ + π β (π(ππΈ + π)) = π β (π πff π»π) + π β― (8) Energy equation is given by equation (8) where π πff , π β― are effective conductivity and volumetric heat source term to account for energy transfer during phase change respectively. Source term π β― is given by equation (15). Transport equation for phase fraction ππΌ πΎ ππ‘ + π΅ β (πΌ πΎ π) = 0 (9) β πΌ πΎπ©πΎ = 1 (10) nterfaces between different fluids are tracked by solving transport equations (9) individually for each phase which is called volume of fluid (VOF) method given by Hirt and Nicholas (1988). Here, πΌ πΎ is the phase fraction of πΎ π‘β phase. If πΌ πΎ is zero, in that case the cell is empty of πΎ π‘β phase and if πΌ πΎ is one, that cell is filled with πΎ π‘β phase alone and if πΌ πΎ varies between zero and one, there is an interface between πΎ π‘β phase and other phases. Sum of all the phases in a given cell is always equal to one (equation 10). π = β π πΎ πΌ πΎππΎ=1 (11) π = β π πΎ πΌ πΎππΎ=1 (12) Mixture properties used in continuity, momentum and energy equations is given by equation (11 & 12) which is phase fraction weighted linear sum of individual property of different phases. Phase change model
Phase change model utilized in the study considers both heterogeneous and homogeneous mode of boiling (Liu et al., 2015). Equation (13) and equation (14) has been employed to account for heterogeneous and homogeneous boiling mode respectively. π π2 = (π π πΌ π + π π£ πΌ π£ )(βπ β βπΌ π )πΏ (13) π π1 = ππΌ β π β (π β π π π π ππ‘ ) , π β₯ π π (14) π β― = βπ π . πΏ (15) Here, π , πΌ , πΏ , π , is thermal conductivity, phase fraction, latent heat of evaporation and evaporation frequency respectively where subscripts β²πβ² denotes liquid phase and β²π£β² denotes vapor phase. Surface tension model
Surface tension ( π ), has been utilized in the study so as to capture the thermocapillary effect, where π is linearly dependent upon the temperature given by equation (16), where, π is the surface tension at the saturation temperature and πΎ is the surface tension gradient (equation (17). π(π) = π β πΎ(π β π π ) (16) πΎ = β ππππ (17) A finite volume based commercially available computational fluid dynamics software ANSYS Fluent is used to carry out the numerical simulation for the problem. All the fluids are assumed to be incompressible, Newtonian and immiscible. Volume of fluid method (VOF) is employed to rack different interfaces. Pressure-Implicit with Splitting of Operators (PISO) scheme is used for pressure-velocity coupling. For spatial discretization, pressure staggering option (PRESTO) for pressure, least squares cell based for gradient and Geo-Reconstruct for volume fraction, second order upwind scheme for momentum is used. For transient discretization, a second order implicit scheme is used.
3. Results and discussions
An attempt is made to develop a model for boiling of liquefied Nitrogen (LN ) which has a very low boiling point of 77 K at atmospheric pressure on a solid as well as liquid surface. Boiling over a liquid medium is quite different from boiling on a solid surface. One of the important parameters in case of boiling on a solid surface, imperfections, is absent for the case of boiling on a liquid surface. In the absence of imperfections, boiling phenomena which generally occurs at the saturation temperature may be delayed up to an upper limit called superheat limit ( T SL ) of the liquid. Boiling on a liquid surface can happen anywhere between saturation temperature ( T S ) and superheat limit depending upon local flow conditions, thermal properties and interfacial interactions. Such superheating of the boiling fluid may lead to instantaneous phase change which may appear explosive some time in nature. A proper understanding of this superheating phenomena and consequently explosive boiling and the critical parameters which affect this behavior is of paramount importance for the safety concerns. Hence, in the following sections, boiling of a cryogenic fluid on a solid surface and subsequently on a liquid surface has been discussed. Numerical computation of boiling on a horizontal solid surface has been reported extensively before for various fluids such as water and others (Son and Dir, 1997; Tomar et al., 2005), but boiling of a cryogen has not been targeted by many. The very first attempt to investigate boiling of a cryogen on a horizontal surface was tried by Liu at al. (2015), where they have predicted the boiling nature of LN for various wall superheats. They have also developed a model for boiling of a cryogenic liquid where two different bubble formation routes have been explained. The first route considers the vapor generation at the horizontal surface which occurs in the case of heterogeneous boiling where the vapor present at the initial stage (t = 0 s) continues to grow due to evaporation at liquid-vapor interface. The other route takes into account the vapor generation in the bulk cryogenic liquid which generally occurs in case of homogeneous boiling where vapor is generated due to superheating of the bulk fluid whenever it comes in direct contact with the horizontal solid surface. The mass transfer models for the evaporation at the liquid-vapor interface (heterogeneous boiling) and evaporation in the bulk fluid (flash or homogeneous boiling) are expressed as: π π1 = π πΏ = β(π π πΌ π + π π£ πΌ π£ )βπ β βπΌ π πΏ (1) π π2 = π πΏ = βπΌ π π π πΆ π,π (π π β π π )πΏ π (2) where, π , πΌ , πΏ , πΆ π and π denotes heat conductivity, volume fraction, latent heat of evaporation, specific heat and time constant, respectively. Simulations have been performed by Liu et al. (2015) for 8 different wall temperatures covering nucleate, transition and film boiling which further resulted in a boiling regime map for LN . Significant presence of small bubbles (micro bubbles) in case of nucleate boiling have been reported which later diminishes and later on converted into film at higher temperatures. They have also investigated the effect of depth of LN and found that the heat flux generation is reduced with the increase in the level of LN . A similar study was also performed by Ahammad et al. (2016a) where they have employed CFD model to investigate saturated film boiling of LN as well as LNG, considering it as pure methane. Growth and departure patterns of vapor bubbles for various wall superheats for both i.e. LN and LNG were discussed. In their subsequent work (Ahammad et al. 2016b), film boiling heat transfer is predicted in case of boiling of a cryogen on a horizontal solid surface for both LN and LNG. In both the above studies, a Rayleigh-Taylor instability (RTI) based model (Zuber et al., 1959) for saturated film boiling on a horizontal solid surface have been employed. The most dangerous wavelength ( π π ), as selected by Liu et al. (2015) is β2 π π whereas the same has been taken as β3 π π by Ahammad et al. (2016). Both the studied were conducted for single mode boiling, and it was assumed that the behavior of the vapor bubbles will be repeated throughout the computational domain. No work has been reported to study multi-mode boiling heat transfer till date. In present work, boiling of a liquid nitrogen (LN ) is simulated on a horizontal solid surface for both single and multi-mode occurrences. Role of different flow instabilities such as Rayleigh-Taylor, Benard-Marangoni and Rayleigh-Benard is discussed, in the context of boiling. Calculation of wall heat flux and Nusselt number in case of boiling on a horizontal solid surface is discussed, in details. Bubble growth behavior and its departure pattern is also reported for different degree of superheats. In Multi-mode boiling, asymmetricity in the bubble growth and its departure is also highlighted during departure of 2 nd batch of vapor bubbles. Flow Instabilities in boiling of a cryogen
Flow instabilities play major role in boiling of a cryogenic fluid and hence, consideration of same in the development of a numerical model is essential. Two important hydrodynamic instabilities that have significant role in boiling are, Rayleigh-Taylor instability (RTI) and Kelvin-Helmholtz instability (KHI). Presence of a denser fluid over a lighter fluid causes instability along the interface which grows over time until the lighter fluid pinches off. This instability is observed during saturated film boiling and called Rayleigh-Taylor instability (RTI). As discussed before, based on the RTI, Zuber et al. (1959) was first to develop a model for saturated film boiling on a horizontal surface. He came up with two important wavelengths, most critical wavelength ( π π ), below which there will not be any bubble formation, most dangerous wavelength ( π π ), distance between two adjacent vapor bubbles for which the vapor generation will be maximum. In KHI, when a lighter fluid floats over a denser fluid and the fluids have different velocities, then instability is induced along the fluid-fluid interface which may grow in time, this can be observed in case of spillage and pool spreading of LNG over water surface, where LNG (lighter fluid) spreads over the water surface (denser fluid) upon spillage. part from hydrodynamic instability, temperature induced flow instability also plays important role in case of film boiling on a horizontal surface. Two important thermal induced flow instabilities are Benard-Marangoni instability (BMI) and Rayleigh-Benard instability (RBI). In BMI, when a very high temperature is applied to a very thin film of capillary length scale, a surface tension gradient is induced along the interface which tries to pull the interface from the region of low surface tension to high surface tension. Thermocapillary instability has a major role in vapor film collapse during saturated film boiling of a cryogenic fluid on a horizontal solid surface (Aursand et al., 2018). Marangoni number (Ma) is used to characterize the BMI and given by ππ = πΎββπππ , where πΎ is the surface tension gradient with respect to temperature, β is the film thickness, βπ is the temperature difference along the interface, π is kinematic viscosity and π is thermal conductivity. Unlike BMI, Rayleigh-Benard instability which occurs in the bulk fluid of considerable depth, here, a density gradient is induced along the depth of the bulk fluid induced by the temperature gradient. This instability causes a local Rayleigh-Benard convection in the region, as fluid particles with high temperature becomes less dense and tries to move up and replace the denser cold particles lying above. One must note that this instability doesnβt play any role in the film boiling of a cryogen on solid surface rather itβs important in the case of boiling on a liquid surface but discussed here for coherency. Rayleigh number characterizes the RBI, which is given by π π = πΌπβ βπππ . Here πΌ , π , β , βπ , π and π are thermal expansion coefficient, acceleration of gravity, film thickness, temperature difference in the region, kinematic viscosity and thermal conductivity, respectively. Wall Heat flux
In the context of random bubble formation in boiling over a surface, the area weighted average heat flux can be evaluated by:
1π΄ β« πππ΄ = 1π΄ β π π |π΄ π | ππ=1 (1) where, π΄ , π , π and π are total area of the wall, heat flux associated with the wall, and total number of facets and individual facet, respectively. Also. Nusselt number (Nu), which is an important dimensionless number characterizing the boiling heat transfer with respect to wave length ( π ) as length scale, can be calculated as: ππ’ = |π"|π π π (π π€πππ β π π ) (2) Figure 1. Shows the variation of area weighted wall heat flux derived from present simulation for a wall superheat of 103 K along with acceptable range of validation against the results of Ahammad et al. (2017). At t = 0.09 s and 0.12 s of first bubble cycle, film thickness is minimum and maximum, respectively, and subsequently these are the maximum and minimum values of heat flux and Nusselt number in Figure 1. Figure 1:
Wall heat flux variation for the present simulation as compared to Ahammad et al. (2017)
Bubble growth and departure
Transient evolution of vapor bubble for single mode boiling is shown in Figure 2 for a wall superheat of 223 K. Initially at t = 0.05 s, the vapor film is stable as the gravitational force dominates over the buoyancy force. With the temporal growth of the initial perturbation of nitrogen vapor due to evaporation of LN at the liquid-vapor interface as shown in Figure 1(a), nitrogen (a) (b) (c) (d) Figure 2 : Growth of vapour bubble and its departure for a wall superheat of 223K. q ( W / m ) Time (s)
Present Simulation
Time weighted average
Ahammad et al. apor moves towards the nodal point where the initial thickness is maximum. This movement of nitrogen vapor towards the nodal point helps in growth of initial disturbance as can be seen in Figure 1(b). Disturbance continues to grow as the heavier fluid LN lies above the lighter fluid N , a classic case of Rayleigh-Taylor instability (RTI). As the vapor bubble grow further in size (Figure 1(c)), necking appears at the base of the vapor bubble, which eventually leads to the pinch-off of the vapor bubble (Figure 1(d)). As discussed before, the film thickness is minimum just before the departure of the vapor bubble which also corresponds to the maximum Nusselt number and heat transfer coefficient. Similar behavior was also reported by Ahammad et al. (2016) and Liu et al. (2015) in their studies. Effect of wall superheat (a) (b) (c) Figure 3: (a)
Rate of bubble generation, (b) change in bubble diameter, and (c) heat flux generated due to phase transfer due to increase in wall temperature. Simulations are performed for different wall superheats and its effect on vapor bubble frequency, bubble diameter and heat flux are reported (Figure 3). Both bubble frequency and bubble diameter increases linearly with the increase in the wall temperature (Figure 3(a & b)), since increasing the
100 150 200 250 300 B ubb l e f re q e un c y , f ( H z ) Wall temperature, T w (K) PresentAhammad et al. (2016)
100 150 200 250 300 B ubb l e d i a m e t er , d ( m ) Wall temperature, T w (K) Present
Ahammad et al. (2016) H e a t f l u x , q ( W / m ) Wall temperature, T w (K) Present
Ahammad et al. (2016) all temperature increases the rate of evaporation at the liquid-vapor interface, which causes faster growth and departure of bubbles. Area weighted average heat flux, π , generated at the bottom wall also increases linearly with the increase in the wall temperature since the heat flux is directly dependent on the temperature difference between wall and the saturation temperature of boiling fluid. All the results, i.e. bubble frequency, bubble diameter and heat flux are compared against the results of Ahammad et al. (2016) and similar behavior were also reported by them. Multi-mode boiling and asymmetry in bubble departure (a) (b)
Figure 4:
Growth and behaviour of vapour bubbles for multi-mode boiling on a horizontal solid surface for two different wall temperature (a) 130 k, and (b) 230 K Most of the studies conducted (Liu at al., 2015; Ahammad et al., 2016) for boiling of a cryogenic fluid on a horizontal solid surface was in single mode. Present simulation, on the other hand, for the first time, simulates a multi-mode boiling of a cryogenic fluid on a horizontal solid surface. All the single mode studies assume that the behaviour of vapour bubbles will remain steady and uniform throughout, but that is not the case in experiments. Therefore, a multi-mode boiling case for the same conditions is simulated for a computational domain of size π/2 x π/2 mm , where π π is the most dangerous wavelength (Zuber et al., 1959). Boundary conditions are same as before, i.e., pressure outlet at the top, no-slip at the bottom wall, and symmetry at the adjacent walls. A constant temperature thermal boundary condition is given at the bottom wall. T w = 130 K T w = 230 K (a) (b) Figure 4: (a)
Asymmetry in bubble growth and the departure of 2 nd batch of vapour bubbles, and (b) chaotic behaviour before merging with the bulk fluid in multi-mode boiling, for a wall superheat of 223 K. Growth and departure for the 1 st batch of vapor bubbles and subsequent batches are shown in Figure 4 (a & b) for two different wall temperatures of 130 K and 230 K. As evident from the phase contours, the first batch of vapor bubble grow and depart in a similar fashion as compared to single mode without interfering with the behavior of their neighbors. The 2 nd batch of vapor bubbles have differences in their sizes and departure time, which is more evident for higher wall superheat. After the departure of 2 nd batch of vapor bubbles, it becomes more chaotic as can be seen in the Figure 4(a & b) just before the departure of 3 rd batch of vapor bubbles. (1) (2) (3) (4) (5)(6) (7) (8) (9) (10) he 2 nd batch of vapor bubbles compete with their neighbours for growth and departure. Figure 4(a) shows the behaviour of two neighbouring bubble vapour sites for various time steps which covers their growth and asymmetry in departure. The wall superheat for this case is 223 K. Growth and departure of 1 st batch of vapor bubbles occur in an idealized condition, where, each and every perturbation are of equal size, moreover the temperature distribution in the perturbations are also same. But, that is not the case for 2 nd batch of vapor bubbles, where initial perturbations are not of the same size and have different amount of vapor content, and the temperature distribution is also different for every vapor sites. (a) (b) Figure 5:
Presence of turbulence fluctuations causing asymmetry in bubble departure Apart from this, due to the pinch-off of the 1 st batch of vapor bubbles, a local turbulence fluctuation in the flow field is created which causes asymmetry in the bubble growth for the 2 nd batch of the vapor bubbles, which is also a reason of asymmetry in bubble growth and their departure (Ruth, 2019). Moreover, the temperature distribution in the film thickness after the departure of 1 st batch of vapor bubbles is not uniform anymore which may have also caused uneven vapor generation along the liquid-vapor interface. In the Figure 4(a), at t = 0.276 s, because of the local turbulence fluctuations created by the departure of 1 st batch of vapor bubbles, shape of bubbles A and B are quite different, bubble B is more spread out in the form of a plateau whereas bubble A has conical shape, which can also be seen at t = 0.288 s. At t = 0.324 s, bubble B also takes the conical shape, but bubble A clearly is more voluminous. With the progression in time, as the bubbles become larger in size, they also become asymmetric in the latitudinal direction (t = 0.348 s) due to the turbulence fluctuations. Figure 5 (a & b) shows the streamlines near the bubbles A and B at t = 0.348 s. Notably, the vorticities near the base of bubble A (Figure 5(a)) is stronger as compared to that of bubble B (Figure 5(b)). The same can be attributed to the local turbulent fluctuations. If we consider only bubble A, uneven fluctuation across it in latitudinal direction, as indicated by region R1 and R2 in Figure 5(a), may have caused the asymmetry for the bubble A. similar behavior can also be observed for bubble B, as specified by region R1β and R2β in Figure 5(b). At t = 0.348 s, necking also appears early for bubble A, necking grows further and becomes millimetric in size at t = 0.360 s, difference in size of neck between them is quite large and finally at t = 0.372 s, bubble A pinches off before the bubble B which takes 0.012 seconds more to pinch off afterwards. R1 R2 R1β R2β .2.
Boiling of a cryogen on liquid surface
Study of cryogenic fluid on a liquid surface is of paramount importance considering the variety of applications of concern such as cryogen spills on water. In such cases, it is crucial to avoid the collapse of stable vapor film, since collapse of vapor film leads to direct contact of boiling fluid and the heating fluid. In the absence of any bubble-forming nuclei, such contacts may lead to superheating of the boiling fluid. Superheating may cause rapid change of phase at an instant which appears explosive in nature, and if the vapor is flammable and there is also a presence of light source, it can be really detrimental for the safety concerns of general well-being. Hence, the boiling regime to be continuously in film boiling regime is desirous. So it is very important to identify the critical parameters which affect this collapse, also, in the case of any collapse in the vapor film, finding some innovative ways to either avoid this collapse or in the worst case cease the vapor diffusion. Boiling of a cryogenic fluid on a solid surface differs from the boiling on a liquid surface, and these differences can be categorized into two parts, (a) fluid dynamics, and (b) thermodynamics. Considering the fluid dynamics first, the normal thrust of the vapor phase will be different for boiling on a solid than that on a liquid surface since the liquid surface is deformable in nature. Secondly, considering the thermodynamics, the rate of diffusion in case of solid will be faster as compared to that of liquid. Apart from this, in case of liquid-liquid boiling, due to the presence of temperature gradient in the bulk fluid, a density gradient will be induced which will cause a local Rayleigh-Benard convection few millimeters below the hot liquid surface (Kumar et al., 2020). Boiling of a cryogen on a liquid surface is simulated here in single mode for a computational domain of width π π/2 and a length of π/2 . The liquid surface chosen for the study is water, and LN is the cryogen. Three phases in the study include water, LN and N vapor . Initially, water phase is filled up to a height of 15 mm from the bottom of the computational domain, over which a N vapor is patched given by equation (2). Rest of the space is occupied by the LN phase. Pressure outlet boundary condition is used on top, no-slip boundary condition at the bottom wall, symmetry at the adjacent walls, and all the walls are taken to be adiabatic. Initially, water phase and LN are taken to be of uniform temperature of 300 K and 77 K respectively, whereas the temperature distribution in the N vapour is given by: π π¦ = {π π β π₯π β π¦πΏ , π¦ β€ πΏπ π , π¦ > πΏ } (17) Figure 6 shows the transient evolution of the phase contours as obtained from simulation. At t = 0.021 s, initially patched vapor phase receives the heat diffused from the water surface which leads to evaporation of LN to at the liquid-vapor interface. The vapor film thickness continuously reduces, as can be seen from the Figure 6, and at t = 0.320 s, it reaches its minimum value before collapse, and finally at t = 0.325 s, the vapor film collapses. After the breaking of vapor film at t = 0.0325 s, the cryogen LN comes into direct contact with the liquid surface, water, and in the absence of bubble-forming nuclei, causes the superheating of LN , up to an upper limit called superheat limit ( π ππΏ ). Figure 6: Transient evolution of phase contours of boiling of LN on water surface (a) (b) Figure 7: (a) Change of vapor generation rate, before and after the vapour film collapse (b) low pressure zone created in the water phase during rise of vapor bubble On superheating, the liquefied nitrogen evaporates at a faster rate called homogeneous boiling. The amount of vapor generation before (inset) and after the vapor collapse as expected is very high which is shown in Figure 7(a). Due to homogeneous boiling, micro vapor bubbles are generated in the domain, which was also observed by Liu at al. (2015) in their study. These micro-bubbles either merge together to form a new macro bubble or join with the existing macro vapor bubble (Figure 6(g)). As the vapor bubble rises up, a low pressure zone is created near the water-N interface (Figure 7(b)), hence the water just below the vapor bubble also rises along the vapor bubble. Moreover, it also propels the LN to fill the void created by rising vapor bubble in the low pressure zone, which pulls the water phase along due to viscous drag. At t = 0.135 s, vapor bubble V o l u m e ( mm ) Flow time (s) t = 0.0325 sBefore vapor film collapse After vapor film collapse V o l u m e ( mm ) Flow time (s) [Pa] inally detaches with the water phase and moves towards the outlet, below which smaller bubbles are also rising slowly towards the outlet. Vapor film collapse
Vapor film collapse occurs whenever the vapor-liquid interface temperature becomes lesser than the Leidenfrost temperature of the boiling fluid. As discussed in the previous section, vapor film collapse leads to direct contact between the boiling fluid and heating fluid. Since the bubble-forming nuclei are absent in case of a liquid surface unlike solid surfaces, the contact leads to superheating of the boiling fluid, which causes rapid change of phase. Hence, predicting the Leidenfrost temperature for various flow conditions is very important. Several studies have been conducted in the past to predict the Leidenfrost point. Berenson (1961) developed a correlation for the minimum temperature difference required to avoid the film collapse in case of on a horizontal solid surface based on the Taylor-Helmholtz hydrodynamic instability, but it was verified with the experimental results of only two fluids, n-pentane and carbon tetrachloride within Β±10 %. Spiegler et al. (1963) derived an expression for Leidenfrost point using Van der Waals equation, which is given as π π , where π π is the critical temperature. Also it was assumed that the boiling fluid is already superheated to its upper limit i.e. superheat limit temperature ( π π π ) at the initial stage. (a) (b) (c) Figure 8: (a) Temperature field and streamlines just before and after the vapor film collapse (b) variation of minimum film thickness ( β πππ ) before rupture (c) change in evaporation rate with time h m i n ( mm ) Time (s)
Film collapse E va p o r a t i o n r a t e ( k g / s ) Time (s) E va p o r a t i o n r a t e ( k g / s ) Time (s)
Before vapor film collapse aumeister and Simon (1973) developed a method to predict the Leidenfrost point on a horizontal solid surface for cryogens, which included effects of critical temperature of the boiling fluid, thermal properties of the solid surface, surface energy of both, solid and liquid. The correlation given for the Leidenfrost point is (π π β π) Γ [0.16 + 2.4 Γ ( π π πΆ π,π π π π π€ πΆ π,π€ π π€ ) ] , where subscript β l β indicates boiling fluid and β w β indicates heating fluid, water, whereas, π π , π , πΆ π,π , π are critical temperature, density, specific heat and thermal conductivity, respectively. Recently, Aursand et al. (2018) developed a theoretical model to predict the Leidenfrost point for saturated film boiling on a horizontal surface, where they have used a non-equilibrium evaporation model based on kinetic theory which has also included thermocapillary effects along the evaporating liquid-vapor interface. They have applied a linear stability analysis on the developed model and found that, a film with a very small and a very large film thickness are always unstable, whereas in the case of intermediate film thickness, instability is dependent on the relative dominance of thermocapillary and vapor thrust, where thermocapillary is destabilizing and vapor thrust is stabilizing in nature. They concluded that the mechanism behind the film boiling collapse may be the thermocapillary instability at the liquid-vapor interface. In the present case of boiling of LN on water surface, the vapor film collapses at t = 0.0325 s. Figure 8 shows the temperature field and streamlines for three different instances. At t = 0.032 s, the fluctuation just below the N vapor in water phase, indicates the point of inception of film collapse. Reduction of vapor film thickness with time, and its collapse is shown in Figure 8 (b). The minimum thickness just before the rupture at t = 0.323 s is 7.62 ππ , which lies in the region of intermediate scales of film thickness as defined by Aursand et al. (2018). The intermediate film thickness scale lies in the range of 0.1 ππ < β πππ < 10 ππ . Since the vapor film collapses in the intermediate film thickness zone, we conclude that the thermocapillary dominates over the vapor thrust. Since the vapor thrust is related to the evaporation rate and the inertia of the vapor mass, and the vapor generation rate is very low in the region (left) before the vapor film collapse as shown in the inset (Figure 8(c)). Conclusions
A numerical study of boiling of a cryogen (LN ) on a solid surface, both single-mode and multi-mode, and also on a liquid surface has been computed. A numerical model based on Rayleigh-Taylor instability has been employed to perform the calculations. In case of boiling on a solid surface, bubble frequency, bubble diameter, heat flux increases linearly with increase in wall superheat. Presence of local fluctuations and turbulence since departure of the 1 st batch of vapor bubbles in case multi-mode boiling on a solid surface, causes asymmetry in bubble growth and departure of 2 nd batch of vapor bubbles. In the absence of sufficient heat flux, and dominance of thermocapillary over vapor thrust, thickness of initial vapor film decreases with time and collapses at t = 0.0325 s. Vapor film collapse leads to a sharp increase in vapor generation and evaporation rate. cknowledgements The authors would like to thank SERB, Department of Science and Technology, Government of India [grant number EMR/2017/000817] for sponsoring this research, and Mr. Lokesh Rohilla for his valuable insights.
References Ahammad, M., Olewski, T., VΓ©chot, L.N. and Mannan, S., 2016. A CFD based model to predict film boiling heat transfer of cryogenic liquids.
Journal of Loss Prevention in the Process Industries , , pp.247-254. 2. Ahammad, M., Liu, Y., Olewski, T., VeΜchot, L.N. and Mannan, M.S., 2016. Application of computational fluid dynamics in simulating film boiling of cryogens.
Industrial & Engineering Chemistry Research , (27), pp.7548-7557. 3. Aursand, E., Davis, S.H. and Ytrehus, T., 2018. Thermocapillary instability as a mechanism for film boiling collapse. 4.
Baumeister, K.J. and Simon, F.F., 1973. Leidenfrost temperatureβits correlation for liquid metals, cryogens, hydrocarbons, and water. 5.
Berenson, P.J., 1961. Film-boiling heat transfer from a horizontal surface. 6.
Blander, M., Hengstenberg, D. and Katz, J.L., 1971. Bubble nucleation in n-pentane, n-hexane and n-pentane+ hexadecane mixtures and water.
The Journal of Physical Chemistry , (23), pp.3613-3619. 7. Cormier, B.R., Qi, R., Yun, G., Zhang, Y. and Mannan, M.S., 2009. Application of computational fluid dynamics for LNG vapor dispersion modeling: A study of key parameters.
Journal of Loss Prevention in the Process Industries , (3), pp.332-352. 8. Carey, V.P., 2020.
Liquid-vapor phase-change phenomena: an introduction to the thermophysics of vaporization and condensation processes in heat transfer equipment . CRC Press. 9.
Drake, E.M., Jeje, A.A. and Reid, R.C., 1975. Transient boiling of liquefied cryogens on a water surface: I. Nitrogen, Methane and Ethane.
International Journal of Heat and Mass Transfer , (12), pp.1361-1368. 10. Enger, T., Hartman, D. E., & Seymour, E. V. (1973).
Explosive Boiling of Liquefied Hydrocarbon/Water Systems. Advances in Cryogenic Engineering, 32β41. doi:10.1007/978-1-4684-3111-7_4. 11.
Gopalaswami, N., Olewski, T., VΓ©chot, L.N. and Mannan, M.S., 2015. Small-scale experimental study of vaporization flux of liquid nitrogen released on water.
Journal of hazardous materials , , pp.8-16. 12. Horvat, A., 2018. CFD methodology for simulation of LNG spills and rapid phase transition (RPT).
Process Safety and Environmental Protection , , pp.358-369. 13. Kumar, R., Rohilla, L. and Das, A.K., 2021. Understanding interfacial behaviour during boiling of nitrogen from liquid-liquid contact plane.
International Journal of Heat and Mass Transfer , , p.120661. 14. Liu, Y., Olewski, T. and VΓ©chot, L.N., 2015. Modeling of a cryogenic liquid pool boiling by CFD simulation.
Journal of Loss Prevention in the Process industries , , pp.125-134. 5. Morse, T.L. and KytΓΆmaa, H.K., 2011. The effect of turbulence on the rate of evaporation of LNG on water.
Journal of loss prevention in the process industries , (6), pp.791-797. 16. Nishigaki, K. and Saji, Y., 1981. The superheat-limit of liquid oxygen and nitrogen at 1 atm.
Japanese Journal of Applied Physics , (5), p.849. 17. Qi, R., Ng, D., Cormier, B.R. and Mannan, M.S., 2010. Numerical simulations of LNG vapor dispersion in Brayton Fire Training Field tests with ANSYS CFX.
Journal of hazardous materials , (1-3), pp.51-61. 18. Reid, R., 1976, Superheated Liquids, American Scientist, 64, 146β156. 19.
Ruth, D.J., Mostert, W., Perrard, S. and Deike, L., 2019. Bubble pinch-off in turbulence.
Proceedings of the National Academy of Sciences , (51), pp.25412-25417. 20. Sinha, D.N., Brodie, L.C. and Semura, J.S., 1987. Liquid-to-vapor homogeneous nucleation in liquid nitrogen.
Physical Review B , (7), p.4082. 21. Son, G. and Dhir, V.K., 1997. Numerical simulation of saturated film boiling on a horizontal surface. 22.
Spiegler, P., Hopenfeld, J., Silberberg, M., Bumpus Jr, C.F. and Norman, A., 1963. Onset of stable film boiling and the foam limit.
International Journal of Heat and Mass Transfer , (11), pp.987-989. 23. Stutz, B. and SimΓ΅es-Moreira, J.R., 2013. Onset of boiling and propagating mechanisms in a highly superheated liquid-the role of evaporation waves.
International Journal of Heat and Mass Transfer , (1-2), pp.683-693. 24. Tomar, G., Biswas, G., Sharma, A. and Agrawal, A., 2005. Numerical simulation of bubble growth in film boiling using a coupled level-set and volume-of-fluid method.
Physics of Fluids , (11), p.112103. 25. VΓ©chot, L., Olewski, T., Osorio, C., Basha, O., Liu, Y. and Mannan, S., 2013. Laboratory scale analysis of the influence of different heat transfer mechanisms on liquid nitrogen vaporization rate.
Journal of Loss Prevention in the Process Industries , (3), pp.398-409. 26. Vesovic, V., 2007. The influence of ice formation on vaporization of LNG on water surfaces.
Journal of hazardous materials , (3), pp.518-526. 27. Zuber, N., 1958. On the stability of boiling heat transfer.
Trans. Am. Soc. Mech. Engrs. ,80