- Email: [email protected]
A CFD based model to predict film boiling heat transfer of cryogenic liquids
Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Contents lists available at ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
A CFD based model to predict film boiling heat transfer of cryogenic liquids chot b, Sam Mannan a, * Monir Ahammad a, Tomasz Olewski b, Luc N. Ve a b
Mary Kay O'Connor Process Safety Center, College Station, TX, USA Mary Kay O'Connor Process Safety Center Extension, Texas A&M University at Qatar, Doha, Qatar
a r t i c l e i n f o
a b s t r a c t
Article history: Received 7 July 2016 Received in revised form 8 September 2016 Accepted 24 September 2016 Available online 28 September 2016
Following an accidental spill of cryogenic liquid (e.g., LNG) on a solid substrate (e.g., concrete), the vapor generation corresponds to different boiling regimes i.e., film boiling, transition boiling, and nucleate boiling. As film boiling phenomena dictate the vapor generation in the early stage of the spill, it is considered as the most important boiling regime in the context of cryogenic (e.g., LNG) source-term estimations. This paper presents CFD simulations of cryogenic film boiling for liquid nitrogen (LN2) and LNG as pure methane. Different aspects of CFD modeling such as vapor-liquid interface morphology, the behavior of heat flux at the heated surface, the effect of wall superheats on bubbles generation frequency and bubbles departure diameter are presented. Based on the results of CFD simulations, a first principle model is applied to correlate the wall heat flux in the film boiling regime. This model can be used to enable a faster estimation of wall heat flux when CFD simulations and use of empirical correlations are not feasible. © 2016 Elsevier Ltd. All rights reserved.
Keywords: Computational fluid dynamics (CFD) Source-term modeling Film boiling Liquefied natural gas (LNG) Bubble generations
1. Introduction Due to the technological break-through of shale gas exploration, LNG industries in the USA and other countries are booming (FERC, 2016). Once an importing country, USA is turning to be LNG exporting nation. As a result, LNG manufacturing activities related to liquefaction, storage, transportation and re-gasification have increased. Experts (Havens and Spicer, 2007) and government agencies (United States Government Accountability Office, 2007) emphasized the need for better risk assessment procedures associated with these operations. Standards such as (NFPA 59A, 2013), therefore, recommend accurate estimation of LNG spill consequence via validated models. Consequences of a LNG release include asphyxiating vapor cloud, fire (flash and pool), vapor cloud explosion (in highly congested area), and rapid phase transition (RPT). It is to be noted that a RPT occurs due to very fast evaporation of LNG on water, causing a pressure wave, is unique to the spill on water and therefore will not occur during a spill on land (Bubbico and Salzano, 2009). A loss of containment of storage tanks, pipes or hoses carrying
* Corresponding author. E-mail address: [email protected] (S. Mannan). http://dx.doi.org/10.1016/j.jlp.2016.09.017 0950-4230/© 2016 Elsevier Ltd. All rights reserved.
cryogenic liquid e.g., LNG will result in the formation of the liquid pool on spilled substrate. Heat transfer from the surrounding to the pool results in vapor formation which disperses and may cause flammable and/or asphyxiating cloud in the dispersed area. Ignition condition and congestion level in the dispersed area dictates whether the dispersed cloud may escalate to consequence type such as vapor cloud explosion, flash fire or pool fire. The size of the cloud, or in other words, the amount of vapor formed from the boiling pool, determines the severity of the consequence phenomena. Thus, estimations of vapor formation from a boiling pool, also known as source-term, an input to dispersion model, is the key to accurate estimation of consequence severity of a loss of containment event (Webber et al., 2010). The sources of heat to an on-land cryogenic pool are conductive heat from the ground, convective heat from the atmosphere and solar radiation or radiative heat from the fire. Among the different heat sources, conductive heat contributes the most to vapor forchot et al., 2013). In mation during the early stage of the spill (Ve general, three different regimes of the conductive heat are recognized during a spill of a cryogenic liquid: film boiling at the beginning, nucleate boiling towards the end of the spill and transition boiling between these two. A pictorial depiction of boiling phenomena, owing to the conductive heat transfer from the substrate, is presented in Fig. 1. A large temperature gradient during
248
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Fig. 1. Cryogenic spill phenomena and associated boiling regimes.
the early stage of the spill (189 K for LNG assuming ground at 300 K), result in a persistent vapor film between the substrate and the boiling liquid, and is known as film boiling. As shown in the schematic diagram under curve AB, generated bubbles completely cover the substrate surface forming a continuous film between the liquid and solid, which restrict the heat flux. With time, the surface temperature falls, reaching a lowest heat transfer rate (Leidenfrost point) as shown by point B in Fig. 1. Further decrease in temperature gradient cannot sustain a continuous film, the film breaks off, and boiling liquid comes in contact with the heated substrate, which is known as transition boiling regime. The behavior of heat flux as a function of wall superheat and schematic diagram of transition boiling are presented by the curve BC and corresponding area below, as shown in Fig. 1. With the contraction of intermittent vapor film, more and more liquid participates in convective heat transfer from the heated surface, causing an increase in heat flux as a function of decreasing wall superheat. The duration of boiling in this regime is very short compared to film boiling and nucleate boiling (CD). As a result, the cumulative amount of vapor generation in transition boiling regime is quite small compared to other boiling regimes. With the further decrease of wall superheats, in the nucleate boiling regime, bubbles become isolated from each other and detach from the nucleation sites of the substrate. Geometrical properties of the substrate, such as surface roughness play a major role in estimating vapor generation due to nucleate boiling. Point D in Fig. 1 represents the onset of nucleate boiling (ONB), below which the heat transferred to the pool contributes in free convection. In the context of consequence analysis, consideration of film boiling regime is very important because of the fact that it corresponds to the highest temperature gradient at the beginning of spill. In 2004, U.S. Federal Energy Regulatory Commission (FERC) contracted ABS consulting to develop a case, to identify appropriate methods for performing consequence analysis of LNG release in water to estimate vapor dispersion and radiation hazards distance (FERC, 2004). The proposed methods have later become a “de facto” standard for consequence analysis of LNG release on water (Johnson and Cornwell, 2007). The detail computational methods have considered only the film boiling regime in estimating the vapor generation source-term. For the LNG release on land, 1-D heat transfer model is the most popular among the industries.
However, this model does not consider the effect of heat transfer resistance due to the presence of different boiling regimes in the liquid phase. Thus, the estimation of vapor source-terms results inaccurate estimation of LNG spill consequence severity. Therefore, this study focuses in film boiling regime to facilitate better consequence analysis.
2. Literature review The simplest of the vapor source-term model is based on 1-D heat conduction from the solid substrate. Reid utilized transient 1-D heat conduction equation to propose four different models as each model corresponds to different boundary conditions in estimating vapor generation rates (Reid, 1980). The most ideal of these models assumes a perfect contact between the liquid and solidchot et al., substrate. However (Reid, 1980), himself and others (Ve 2013) criticized this model for unrealistic assumptions and proposed to consider surface-thermal resistance due to the presence of film in the liquid solid interface. This is accounted via estimating heat transfer co-efficient from empirical film boiling correlations. The first notable correlation to estimate film boiling heat transfer co-efficient on a flat boiling surface is credited to (Berenson, 1961). His formulation was based on the original work of Rayleigh-Taylor (R-T) hydrodynamic instability model of (Zuber, 1959) and expressed as Equation (1).
Nu ¼ 0:425
rv ðrl rv Þghfg kv mv DT
1=4
s
3=8
g ðrl rv Þ
(1)
Experimental works of (Holster and Westwater, 1962) confirmed the coherency of R-T instability during the film boiling on a horizontal surface. Water and Freon 11 (CCl3F) systems were studied to validate the predictions of Berenson correlation. A generalization of film boiling correlations on flat surface were proposed by (Klimenko, 1981). This proposed correlation was validated against different liquids including cryogens (Klimenko and Shelepen, 1982). The previously mentioned FERC “de facto” consequence analysis method has considered this correlation to estimate vapor formation source-term. For an upward facing geometrical system, according to this correlation, Nusselt number is predicted as follows:
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Nu ¼ 3:02 102 Ar1=3 Pr 3 f 1 ðbÞ;
for
Nu ¼ 1:37 103 Ar1=2 Pr 3 f 2 ðbÞ;
for
1
1
Ar < 108 Ar > 108
(2) (3)
Where,
f1
¼ 1 for b > 0:71 ¼ 0:89 b1=3 for b < 0:71
f2
¼ 1 for b > 0:5 ¼ 0:71 b1=2 for b < 0:5
C DT Cp;v mv g L3 r ðr rv Þ b ¼ p; v ; Ar ¼ c v 2l ; Pr ¼ ; hfg kv mv rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h Lc sl and characteristic length; Lc ¼ Nu ¼ kv gðrl rv Þ In spite of the fact that both Berenson and Klimenko correlations are very useful in estimating vapor generation source-term, both have fundamental limitations. For example, it is unrealistic to consider a constant film thickness over the solid surface whereas in practical case it is continuously changing owing the bubble generations and the movement of boiling the liquid. Other limitations include assumptions of constant bubble diameter and height, periodic bubble generation from the heated surface. As a result, a temporal aspect of heat flux is lost when correlations were used. Furthermore, local physical properties of the liquid and vapor were replaced by the mean properties. Computational fluid dynamics (CFD) is a great tool for overcoming these limitations and therefore utilized in this research to capture the dynamic nature of film boiling and associated wall heat flux. At the end of the XX century, CFD simulation of film boiling on a horizontal surface was pioneered by (Son and Dhir, 1997). A moving mesh method is applied to study bubble and film dynamics for water (Dhir, 2001). Water at near critical conditions was also studied by the authors to provide steady-state bubble release pattern (Son and Dhir, 1998). A sub-cooled film boiling of water on a horizontal disk was studied by (Banerjee and Dhir, 2001). To overcome the computational difficulties in simulating both low and high density ratio fluids (Juric and Tryggvason, 1998), implemented interfacial source terms in the continuity equations on a Eulerian grid. This numerical recipe was further improved by (Esmaeeli and Tryggvason, 2004a, 2004b). Volume of fluid (VOF) method (Youngs, 1982) was used by (Agarwal et al., 2004; Welch and Rachidi, 2002; Welch and Wilson, 2000) to conjugate heat transfer system for studying saturated horizontal film boiling of water on steel. The film boiling of water was studied at saturation condition and near critical point by (Tomar et al., 2008, 2005) using coupled level-set and volume of fluid (CLSVOF) method. Film boiling of water on a sphere was studied by (Yuan et al., 2008) using a non-orthogonal body fitted coordinates. The first attempts to study boiling of cryogenic liquid were attempted by (Liu et al., 2015) and (Ahammad et al., 2016) using a commercial CFD package, ANSYSFLUENT. Liu proposed different boiling routes to model three boiling phases for a pool boiling system of liquid nitrogen. In spite of a decade progress in CFD application to model film boiling, no attempt has been made to study LNG system. This study is primarily focused on CFD applications to study film boiling cryogenic liquids i.e., LNG (as pure methane) and liquid nitrogen (LN2).
3. The film boiling model The Rayleigh- Taylor (R-T) instability approach (Zuber, 1959) is
249
adopted in modeling film boiling using a two-dimensional computational domain. In this approach, bubbles generated under the liquid phase tries to invert phases due to gravity. Bubble spacing, bubble size and bubble frequency can be used to estimate the vapor formations. These parameters are governed by the hydrodynamic factors of the boiling systems, therefore CFD method is utilized. A perturbation analysis of R-T instability (Carey, 1992) suggested that for a horizontal system, the highest amount of bubble generation would occur if the distance between two nearest bubble is equal to the “most dangerous wavelength” as given by the following Equation.
ld ¼ 2p
3s ðrl rv Þg
1=2 (4)
Fig. 2 graphically present the underlying concepts of R-T approach when adopted for CFD simulation of film boiling. On the left hand side, it shows that bubble generation from the sustained vapor film which rises against the gravity force. It is assumed that the vapor film thickness is large in compared to the substrate surface roughness. The location of the surface where a bubble evolves and detaches from the film is considered as a node point. Between two nearest bubble, the point where the film thickness is lowest is considered as an anti-node point. The dotted rectangle contains four nodes at its vertices and one in the center. Among these five nodes, there are four anti-nodes as shown on the right hand side of Fig. 2. Assuming uniformity of surface, i.e., bubble generation is consistent through-out the heated surface, a domain enclosed by a node point and anti-node point would suffice for this study. Thus, a domain of ld/2 3/2 ld is used where a constant wall temperature is applied as a boundary conditions in the heated surface. The initial vapor film is a sinusoidal perturbation to induce R-T instability and expressed as Equation (5). A linear temperature gradient (Equation (6)) is applied in the initial vapor film and saturation temperature is used elsewhere in the initial temperature field.
d¼
ld 64
Ty ¼
4 þ cos
2px
ld
Twall DT·y=d Tsat
for a ¼ 1 for a ¼ 0
(5)
(6)
3.1. Simulation conditions The computation of the film boiling model is performed by a commercial CFD package ANSYS FLUENT, version 14.0. A 2D computational domain of grid size of 5.79 mm 17.37 mm for LN2 and 9.78 mm 29.33 mm for LNG were used. Grid and time step size independence test were performed for both cases and it is found that a structural grid of 64 192 segments with a time step of Dt ¼ 104 s was sufficient to capture the bubble dynamics with reasonable accuracy. Temperature dependent physical properties of LN2 (Barron, 1999) and LNG as pure methane (Rowley et al., 2003) are used. Symmetric boundary conditions were applied on the both sides of the computational domain as shown in Fig. 2. The volume of Fluid (VOF) method is used to track the vapor-liquid interface. In this method, mass, momentum and energy equations are solved once with and the additional equation for computing the volume fraction (a) of each computational cell. Vapor is considered as the primary phase. Geometric reconstruction algorithm is used to better estimate the bubble interface. To account the surface tension force which stabilizes the interfacial instability, a continuum surface force (CSF) model is used (Brackbill et al., 1992). First order
250
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Fig. 2. Two-dimensional CFD set-up for studying film boiling of cryogenic liquids.
implicit discretization is used for transient time formulations and VOF equation. Energy and momentum equations were discretized using QUICK algorithm. PRESTO procedure is used to interpolate pressure values. PISO algorithm is used to solve the corresponding governing equations. 4. Results and discussions Liquid nitrogen (LN2) system with wall superheats of 32 K, 73 K, 103 K, 153 K and 196 K and LNG boiling system with wall superheats of 43 K, 83 K, 123 K and 163 K were simulated by following
aforementioned approach. The results of these simulations are discussed as follows: 4.1. Individual bubble growth Fig. 3(a) e (d) present evolution of bubble from the initial vapor film for a LN2 system of 103 K wall superheat. In these Figures, x ¼ 5.79 mm corresponds to node point. The initialized vapor film (by Equations (5) & (6)) stabilizes at the beginning and then starts growing. Vapor generates at the interface as the liquid evaporates. As vapor generates, it increases the film thickness and moves
Fig. 3. Evolution of bubbles for a wall superheat of 103 K for LN2 system.
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
4.2. Behavior of wall heat flux The heat flux from the solid wall during the film boiling strongly dependent on the thickness of the film as also observed in other studies (Agarwal et al., 2004; Tomar et al., 2005). The behavior of heat flux can be explained by the behavior of the sustained film thickness during the bubble evolution. As vapor phase presents greater thermal resistance than the liquid phase; and thus these locations where the film thickness is smaller have higher heat flux than that of a thicker film. From Fig. 4, the heat flux decreases as the film thickness increases steadily in the beginning of bubble generation. At 0.05 s the vapor reaches at the most stable uniformly distributed condition. Vapor start pushing the bubble at this time and slowly builds up the bubble. As vapor from the film moves towards the bubble, the thickness of the film decreases on the surface because of the vapor deficiency. As a result the area weighted average heat flux increases with time. At around 0.18s (see Fig. 4 (a)) the bubble grows to its maximum size and is about to leave the film. At this point, the average film thickness on the surface is minimum, resulting maximum heat flux. Immediately after bubble detachment, the retreated vapor from bubble stem takes some time to stabilize. Because of the presence of wake behind the bubble, the interface again necks slightly and generates a smaller peak as seen in the heat flux curve seen around time from 0.21s to 0.31s. This phenomenon repeats causing large peaks representing the bubble release and troughs as points of highest film thickness, (d) in Fig. 4. The CFD simulated heat flux of this work is
compared with corresponding estimations using Berenson and Klimenko correlations in Fig. 4. It is found that the time weighted average (TWA) simulated wall heat flux is bounded by the correlations for the case of wall superheat of 103 K. Fig. 5 depicts the trend of wall heat flux for LNG at a wall superheat of 83 K. Similar to LN2, the peaks represent the liberation of bubbles which follows a period of bubble growth. The average wall heat flux is obtained from the simulations is 23 kW/m2 which is very high compared to the estimations of Klimenko and Berenson correlations.
4.3. Effect of wall superheats on bubble generation frequency and bubble diameter Fast Fourier Transform (FFT) is performed on 16,384 data points of the simulated heat flux to determine the frequency of the bubble generation. Fig. 6 depicts the results of the FFT analysis for the simulated case of LNG with a wall superheat of 43 K. The dominant frequency of the bubble generation which corresponds to the x-axis value (Fig. 6) of the highest peak has been found to be 4.272 Hz. Similar analysis has been performed for LN2 and LNG for all simulated wall superheat conditions. Fig. 7 shows bubble generation frequency for various simulated wall superheats. It is clear that the bubble formation rate is strongly
50000 Simulation
45000
Berenson
Klimenko
TWA-Simulation
40000 35000
q (W/m2)
towards the node as shown in Fig. 3(b). Further generations of vapor in the film tend to move towards to the node, withdrawing a large quantity of vapor from the film, and results in a visible bubble at Fig. 3(c). Then, the average film thickness decreases as the bubble grows and starts departing. When the buoyancy force inside the bubble cannot be compensated by surface tension and gravity, bubble detach from the film Fig. 3(d). Due to surface tension force and capillary action, some portion of vapor from the bubble stem returns to the initial film thus increases the film thickness. The average minimum film thickness is usually occurring just before the liftoff of the bubble. Individual bubble growth mechanism for LN2 and LNG are similar and therefore is not described.
251
30000 25000 20000 15000 10000 5000 0
0
0.1
0.2
0.3
0.4
0.5
0.6
Time, s Fig. 5. Surface heat flux for the film boiling of LNG at a wall superheat DT ¼ 83 K.
Fig. 4. Trend of wall heat flux for LN2 film boiling at a wall superheat of 103 K.
252
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Bubble diameter, db (m)
0.01 LN2
0.009
LNG
0.008 0.007
db = 0.0025ΔT0.2465 R² = 0.8578
0.006 0.005
db = 0.001ΔT0.2766 R² = 0.7258
0.004 0.003 0.002 0.001 0 20
60
100
140
Wall Superheat, ΔT (K)
180
Fig. 8. Bubble diameters of the simulated film boiling cases at various wall superheats.
4.4. Proposed model to predict film boiling heat transfer
Bubble Frequency, f (Hz)
Fig. 6. FFT of heat flux variation for LNG at a wall superheat of DT ¼ 43 K.
12 LN2 10
LNG f = 0.4179ΔT0.6036 R² = 0.9514
8 6
f = 0.1233ΔT0.7863 R² = 0.9933
4 2 0 20
60
100
140
180
Wall Superheat, ΔT (K) Fig. 7. Bubble generation frequency of simulated film boiling cases at various wall superheats.
dependent on the wall superheats. Within the simulated window, a model is fitted to predict the bubble generation frequency as a function of wall superheats, and it is found that, bubble frequency of LN2 follows a power of 0.8 whereas LNG follows a power of 0.6. However, the co-efficient is higher for LNG in compared to LN2 in an order of 3.4. Fig. 8 shows the dependency of bubble diameters, estimated at the time of departure from the vapor film. For a particular wall superheat, the bubble size is larger for LNG in compared to LN2. A power model is fitted (shown on Fig. 8) for both LN2 and LNG to predict the bubble diameter as a function of wall superheats. For LNG the exponent of the model has found to be 0.25 and for LN2, its value is 0.28. However, the co-efficient for LNG is 2.5 times higher than that of LN2.
Fig. 9(a) compares the time weighted average (TWA) simulated heat fluxes with the heat fluxes estimated by using Berenson and Klimenko correlations, for a LN2 film boiling system. Similarly, Fig. 9(b) compares that of for LNG film boiling system. For the case of LN2, at a lower wall superheat, the simulated heat flux is closer to the estimation of Klimenko correlations and at higher wall superheat; the simulated heat flux is closer to estimation of Berenson correlations. The trend of the simulated results, depict that a logarithmic slope is much higher than the slope of both correlations. For the case of LNG, the TWA simulated heat fluxes are much higher than the estimations from Berenson and Klimenko correlations. It is important to mention that, the simulated results are not subjected to the assumptions such as constant film thickness, constant bubble diameters, absence of dynamic movement of vapor and liquid, usage of mean fluid properties etc. Thus, the simulated heat fluxes are more realistic than the aforementioned correlations. Film boiling simulation for a single wall superheat takes approximately 800 min of computer time (Processor: 3.33 GHz; RAM: 32 GB) to simulate two consecutive bubble detachments. Therefore, performing simulations for a large range of wall superheats might not be always possible when the time does not allow. Therefore, a first principle model based on simple heat balance is used to predict the wall heat flux for non-simulated wall superheats. The proposed model is based on the heat balance as the amount of heat transferred from the wall is approximately equal or proportional to the heat taken by the bubbles that were formed in the film boiling regime. This can be mathematically expressed as Equation (7).
qfVb $f $m$rv $hfg
(7)
For a two-dimensional CFD simulation, the above equation can be modified as follows:
q ¼ C·
pd2b 4
·f ·m·rv ·hfg
(8)
where, it is assumed that C is a dimensionless constant equals to 1.0. For some cases C might not be 1.0. For an example, if the local physical property data such as r, m and hfg strongly varies with the temperature field, C ¼ 1.0 may not be the proper representation. Furthermore, the bubble diameter db will also vary with the time of simulation and wall superheats, etc. The departure diameter of the first bubble for a particular wall superheat might not be equal to
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
253
q (W/m2)
50000
5000
q (W/m2)
5000
TWA Simulation Klimenko Berenson Proposed model
500 20
Wall superheat, ΔT (K)
Proposed model TWA Simulation Klimenko Berenson
500
200
20
(a) LN2
Wall superheat, ΔT (K) (b) LNG
200
Fig. 9. CFD based correlation to predict heat flux of (a) LN2 and (b) LNG film boiling.
that of the second bubble. The shape of the bubble can be different from spherical shape. Wake formation in the liquid and liquid depth may also affect the size and shape of the bubbles. Therefore, adjustment of the constant C might be needed to accurately predict the heat flux of the non-simulated cases. The parameters in Equation (8) can be estimated using the simulation results. Fourier transformation is employed to determine the bubble generation frequency f from the transient wall heat flux for different wall superheat as depicted in Figs. 6 and 7. The departure diameter of the simulated bubble is considered to estimate the volume of the bubble (for 3D simulation) or area of the bubble (for 2D simulation). Therefore, by replacing the term, f and db using the fitted model as shown in Figs. 7 and 8 and replacing the parameter m, i.e., number of nodes per unit area, with 1/L2c , the proposed model for LN2 and LNG are as follows:
q ¼ 8:02 109
q ¼ 1:7 107
rl rv rv Lv ·DT 1:34 ðfor LN2 Þ s
(9)
Acknowledgments This paper was made possible by a National Priority Research Project (NPRP) award [NPRP 6-425-2-172] from the Qatar National Research Fund (a member of The Qatar Foundation). Nomenclature
s r g
l
rl rv rv Lv ·DT 1:1 ðfor LNGÞ s
the heat flux in the film boiling regime. The proposed model can be used along with in 1-D standard heat transfer model to determine the liquid phase resistance, and thus, enabling more accurate estimation of vapor formation source-term.
(10)
Fig. 9 depicts the fitness of the correlation for the proposed correlation for the two-dimensional simulations performed in this study. It is shown that first principle based correlation predicts the simulated heat fluxes to an excellent agreement. 5. Conclusions Film boiling of cryogenic liquids, i.e., LN2 and LNG, is simulated by adopting Rayleigh-Taylor (R-T) approach in a commercial CFD environment (ANSYS Fluent). It is found that the dynamics of film thickness affects the wall heat flux and therefore the frequency of bubbles evolution and bubbles diameter. The bubble generation frequency is estimated from the transient wall heat flux by decomposing via FFT and the most dominant frequency is considered as the average bubble frequency for that particular wall superheat. In the simulated window, LN2 and LNG have shown different trends in bubble generation frequency and bubble departure diameter resulting from increased wall superheats. The simulated film boiling heat fluxes for LN2 were higher in comparison to the estimates by Klimenko correlations but it underestimates when compared to Berenson correlation. On the other hand, the LNG simulations over-estimates the wall heat fluxes when compared against both Berenson and Klimenko correlations. Finally, based on the trends of bubble generation frequency and bubble diameters, a first principle model is applied for correlating
hfg k
m p Cp T DT t Dt Lc
a h
d x, y q Vb f n db m
Surface tension, N/m Density, kg/m3 Acceleration of gravity, 9.81 m/s2 Wavelength, m Heat of vaporization, J/kg Thermal conductivity, W/m.K Dynamic viscosity, Pa.s Pressure, Pa Heat capacity, J/kg.K Temperature, K Wall superheat, K Time, s Time step, s Characteristic length, m Volume fraction; vapor volume fraction Average local heat transfer co-efficient, W/m2.K Thickness of the initial vapor film, m Axis of the simulated plane, m Wall heat flux, W/m2 Bubble volume, m3 Frequency of bubble generation, Hz Number of nodes per square area Diameter of the bubble just before detachment from the film, m Number of nodes per square meter
Subscripts v Vapor phase l Liquid phase wall, sat Indication of wall condition and saturation condition b Bubble y Y axis direction
254
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Dimensionless numbers f1, f2, b Dimensionless number Nu Nusselt number Ar Archimedes number Pr Prandtl number References Agarwal, D.K., Welch, S.W.J., Biswas, G., Durst, F., 2004. Planar simulation of bubble growth in film boiling in near-critical water using a variant of the VOF method. J. Heat. Transf. 126, 329. http://dx.doi.org/10.1115/1.1737779. Ahammad, M., Liu, Y., Olewski, T., Vechot, L., Mannan, M.S., 2016. Application of computational fluid dynamics (CFD) in simulating film boiling of cryogens. Ind. Eng. Chem. Res. 55, 7548e7557. http://dx.doi.org/10.1021/acs.iecr.6b01013. Banerjee, D., Dhir, V.K., 2001. Study of subcooled film boiling on a horizontal disc: Part IdAnalysis. J. Heat. Transf. 123, 271e284. http://dx.doi.org/10.1115/ 1.1345889. Barron, R.F., 1999. Cryogenic Heat Transfer. Taylor & Francis, Philadelphia, PA, USA; London, UK. Berenson, P.J., 1961. Film-boiling heat transfer from a horizontal surface. J. Heat. Transf. 351e356. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335e354. Bubbico, R., Salzano, E., 2009. Acoustic analysis of blast waves produced by rapid phase transition of LNG released on water. Saf. Sci. 47, 515e521. http:// dx.doi.org/10.1016/j.ssci.2008.07.033. Carey, V.P., 1992. Liquid-vapor Phase Change Phenomena: an Introduction to the Thermophysics of Vaporization and Condensation Process in Heat Transfer Equipment. Hemisphere Publishing Corp, Washington , D.C. Dhir, V.K., 2001. Numerical simulations of pool-boiling heat transfer. AIChE J. 47, 813e834. http://dx.doi.org/10.1002/aic.690470407. Esmaeeli, A., Tryggvason, G., 2004a. Computations of film boiling. Part I: numerical method. Int. J. Heat. Mass Transf. 47, 5451e5461. http://dx.doi.org/10.1016/ j.ijheatmasstransfer.2004.07.027. Esmaeeli, A., Tryggvason, G., 2004b. Computations of film boiling. Part II: multimode film boiling. Int. J. Heat. Mass Transf. 47, 5463e5476. http://dx.doi.org/ 10.1016/j.ijheatmasstransfer.2004.07.028. FERC, 2016. Approved North American LNG Import and Export Terminals [WWW Document]. http://www.ferc.gov/industries/gas/indus-act/lng/lng-approved. pdf (accessed 04.21.16.). FERC, 2004. Notice of Availability of Detailed Computations for the Consequence Assessment Methods for Incidents Involving Releases from Liquefied Natural Gas Carriers. Havens, J., Spicer, T., 2007. United States regulations for siting LNG terminals: problems and potential. J. Hazard. Mater 140, 439e443. Holster, E.R., Westwater, J.W., 1962. Film boiling on a horizontal plate. Am. Rocket Soc. J. 553e558. Johnson, D.W., Cornwell, J.B., 2007. Modeling the release, spreading, and burning of
LNG, LPG, and gasoline on water. J. Hazard. Mater 140, 535e540. http:// dx.doi.org/10.1016/j.jhazmat.2006.10.022. Juric, D., Tryggvason, G., 1998. Computations of boiling flows. Int. J. Multiph. Flow. 24, 387e410. Klimenko, V.V., 1981. Film boiling on a horizontal plate - new correlation. Int. J. Heat. Mass Transf. 24, 69e79. http://dx.doi.org/10.1016/0017-9310(81)90094-6. Klimenko, V.V., Shelepen, A.G., 1982. Film boiling on horizontal plate - a supplementary communication. Int. J. Heat. Mass Transf. 25, 1611e1613. chot, L.N., 2015. Modeling of a cryogenic liquid pool boiling by Liu, Y., Olewski, T., Ve CFD simulation. J. Loss Prev. Process Ind. 35, 125e134. NFPA 59A, 2013. Standard for the Production, Storage, and Handling of Liquefied Natural Gas (LNG). National Fire Protection Association. Reid, R.C., 1980. Boiling of LNG on Typical Dike Floor Materials. Rowley, R.L., Wilding, W.V., Oscarson, J.L., Zundel, J.L., Marshall, N.A., Daubert, T.L., Danner, R.P., 2003. DIPPR Data Compilation of Pure Compound Properties. Design Institute for Physical Property Data/AIChE. Son, G., Dhir, V.K., 1998. Numerical simulation of film boiling near critical pressures with a level set method. J. Heat. Transf. 120, 183e192. Son, G., Dhir, V.K., 1997. Numerical simulation of saturated film boiling on a horizontal surface. J. Heat. Transf. 119, 525e533. Tomar, G., Biswas, G., Sharma, A., Agrawal, A., 2005. Numerical simulation of bubble growth in film boiling using a coupled level-set and volume-of-fluid method. Phys. Fluids 17. http://dx.doi.org/10.1063/1.2136357. Tomar, G., Biswas, G., Sharma, A., Welch, S.W.J., 2008. Multimode analysis of bubble growth in saturated film boiling. Phys. Fluids 20, 1e8. http://dx.doi.org/10.1063/ 1.2976764. United States Government Accountability Office, 2007. Report GAO-07-316maritime Security: Public Safety Consequences of a Terrorist Attack on a Tanker Carrying Liquefied Natural Gas Need Clarification Highlights. chot, L., Olewski, T., Osorio, C., Basha, O., Liu, Y., Mannan, S., 2013. Laboratory scale Ve analysis of the influence of different heat transfer mechanisms on liquid nitrogen vaporization rate. J. Loss Prev. Process Ind. 26, 398e409. http:// dx.doi.org/10.1016/j.jlp.2012.07.019. Webber, D.M., Gant, S.E., Ivings, M.J., Jagger, S.F., 2010. LNG Source Term Models for Hazard Analysis: a Review of the State-of-the-art and an Approach to Model Assessment (No. Research Report: RR789). HSE Books, Buxton, Derbyshire, UK. Welch, S.W.J., Rachidi, T., 2002. Numerical computation of film boiling including conjugate heat transfer. Numer. Heat. Transf. B 42, 35e53. Welch, S.W.J., Wilson, J., 2000. A volume of fluid based method for fluid flows with phase change. J. Comput. Phys. 160, 662e682. http://dx.doi.org/10.1006/ jcph.2000.6481. Youngs, D.L., 1982. Time dependent multi-material flow with large fluid distortion. In: Morton, K.W., Baines, M.J. (Eds.), Numerical Methods in Fluid Dynamics. Academic Press, pp. 273e285. Yuan, M.H., Yang, Y.H., Li, T.S., Hu, Z.H., 2008. Numerical simulation of film boiling on a sphere with a volume of fluid interface tracking method. Int. J. Heat. Mass Transf. 51, 1646e1657. http://dx.doi.org/10.1016/ j.ijheatmasstransfer.2007.07.037. Zuber, N., 1959. Hydrodynamic Aspects of Boiling Heat Transfer. University of California, Los Angeles.
Contents lists available at ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
A CFD based model to predict film boiling heat transfer of cryogenic liquids chot b, Sam Mannan a, * Monir Ahammad a, Tomasz Olewski b, Luc N. Ve a b
Mary Kay O'Connor Process Safety Center, College Station, TX, USA Mary Kay O'Connor Process Safety Center Extension, Texas A&M University at Qatar, Doha, Qatar
a r t i c l e i n f o
a b s t r a c t
Article history: Received 7 July 2016 Received in revised form 8 September 2016 Accepted 24 September 2016 Available online 28 September 2016
Following an accidental spill of cryogenic liquid (e.g., LNG) on a solid substrate (e.g., concrete), the vapor generation corresponds to different boiling regimes i.e., film boiling, transition boiling, and nucleate boiling. As film boiling phenomena dictate the vapor generation in the early stage of the spill, it is considered as the most important boiling regime in the context of cryogenic (e.g., LNG) source-term estimations. This paper presents CFD simulations of cryogenic film boiling for liquid nitrogen (LN2) and LNG as pure methane. Different aspects of CFD modeling such as vapor-liquid interface morphology, the behavior of heat flux at the heated surface, the effect of wall superheats on bubbles generation frequency and bubbles departure diameter are presented. Based on the results of CFD simulations, a first principle model is applied to correlate the wall heat flux in the film boiling regime. This model can be used to enable a faster estimation of wall heat flux when CFD simulations and use of empirical correlations are not feasible. © 2016 Elsevier Ltd. All rights reserved.
Keywords: Computational fluid dynamics (CFD) Source-term modeling Film boiling Liquefied natural gas (LNG) Bubble generations
1. Introduction Due to the technological break-through of shale gas exploration, LNG industries in the USA and other countries are booming (FERC, 2016). Once an importing country, USA is turning to be LNG exporting nation. As a result, LNG manufacturing activities related to liquefaction, storage, transportation and re-gasification have increased. Experts (Havens and Spicer, 2007) and government agencies (United States Government Accountability Office, 2007) emphasized the need for better risk assessment procedures associated with these operations. Standards such as (NFPA 59A, 2013), therefore, recommend accurate estimation of LNG spill consequence via validated models. Consequences of a LNG release include asphyxiating vapor cloud, fire (flash and pool), vapor cloud explosion (in highly congested area), and rapid phase transition (RPT). It is to be noted that a RPT occurs due to very fast evaporation of LNG on water, causing a pressure wave, is unique to the spill on water and therefore will not occur during a spill on land (Bubbico and Salzano, 2009). A loss of containment of storage tanks, pipes or hoses carrying
* Corresponding author. E-mail address: [email protected] (S. Mannan). http://dx.doi.org/10.1016/j.jlp.2016.09.017 0950-4230/© 2016 Elsevier Ltd. All rights reserved.
cryogenic liquid e.g., LNG will result in the formation of the liquid pool on spilled substrate. Heat transfer from the surrounding to the pool results in vapor formation which disperses and may cause flammable and/or asphyxiating cloud in the dispersed area. Ignition condition and congestion level in the dispersed area dictates whether the dispersed cloud may escalate to consequence type such as vapor cloud explosion, flash fire or pool fire. The size of the cloud, or in other words, the amount of vapor formed from the boiling pool, determines the severity of the consequence phenomena. Thus, estimations of vapor formation from a boiling pool, also known as source-term, an input to dispersion model, is the key to accurate estimation of consequence severity of a loss of containment event (Webber et al., 2010). The sources of heat to an on-land cryogenic pool are conductive heat from the ground, convective heat from the atmosphere and solar radiation or radiative heat from the fire. Among the different heat sources, conductive heat contributes the most to vapor forchot et al., 2013). In mation during the early stage of the spill (Ve general, three different regimes of the conductive heat are recognized during a spill of a cryogenic liquid: film boiling at the beginning, nucleate boiling towards the end of the spill and transition boiling between these two. A pictorial depiction of boiling phenomena, owing to the conductive heat transfer from the substrate, is presented in Fig. 1. A large temperature gradient during
248
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Fig. 1. Cryogenic spill phenomena and associated boiling regimes.
the early stage of the spill (189 K for LNG assuming ground at 300 K), result in a persistent vapor film between the substrate and the boiling liquid, and is known as film boiling. As shown in the schematic diagram under curve AB, generated bubbles completely cover the substrate surface forming a continuous film between the liquid and solid, which restrict the heat flux. With time, the surface temperature falls, reaching a lowest heat transfer rate (Leidenfrost point) as shown by point B in Fig. 1. Further decrease in temperature gradient cannot sustain a continuous film, the film breaks off, and boiling liquid comes in contact with the heated substrate, which is known as transition boiling regime. The behavior of heat flux as a function of wall superheat and schematic diagram of transition boiling are presented by the curve BC and corresponding area below, as shown in Fig. 1. With the contraction of intermittent vapor film, more and more liquid participates in convective heat transfer from the heated surface, causing an increase in heat flux as a function of decreasing wall superheat. The duration of boiling in this regime is very short compared to film boiling and nucleate boiling (CD). As a result, the cumulative amount of vapor generation in transition boiling regime is quite small compared to other boiling regimes. With the further decrease of wall superheats, in the nucleate boiling regime, bubbles become isolated from each other and detach from the nucleation sites of the substrate. Geometrical properties of the substrate, such as surface roughness play a major role in estimating vapor generation due to nucleate boiling. Point D in Fig. 1 represents the onset of nucleate boiling (ONB), below which the heat transferred to the pool contributes in free convection. In the context of consequence analysis, consideration of film boiling regime is very important because of the fact that it corresponds to the highest temperature gradient at the beginning of spill. In 2004, U.S. Federal Energy Regulatory Commission (FERC) contracted ABS consulting to develop a case, to identify appropriate methods for performing consequence analysis of LNG release in water to estimate vapor dispersion and radiation hazards distance (FERC, 2004). The proposed methods have later become a “de facto” standard for consequence analysis of LNG release on water (Johnson and Cornwell, 2007). The detail computational methods have considered only the film boiling regime in estimating the vapor generation source-term. For the LNG release on land, 1-D heat transfer model is the most popular among the industries.
However, this model does not consider the effect of heat transfer resistance due to the presence of different boiling regimes in the liquid phase. Thus, the estimation of vapor source-terms results inaccurate estimation of LNG spill consequence severity. Therefore, this study focuses in film boiling regime to facilitate better consequence analysis.
2. Literature review The simplest of the vapor source-term model is based on 1-D heat conduction from the solid substrate. Reid utilized transient 1-D heat conduction equation to propose four different models as each model corresponds to different boundary conditions in estimating vapor generation rates (Reid, 1980). The most ideal of these models assumes a perfect contact between the liquid and solidchot et al., substrate. However (Reid, 1980), himself and others (Ve 2013) criticized this model for unrealistic assumptions and proposed to consider surface-thermal resistance due to the presence of film in the liquid solid interface. This is accounted via estimating heat transfer co-efficient from empirical film boiling correlations. The first notable correlation to estimate film boiling heat transfer co-efficient on a flat boiling surface is credited to (Berenson, 1961). His formulation was based on the original work of Rayleigh-Taylor (R-T) hydrodynamic instability model of (Zuber, 1959) and expressed as Equation (1).
Nu ¼ 0:425
rv ðrl rv Þghfg kv mv DT
1=4
s
3=8
g ðrl rv Þ
(1)
Experimental works of (Holster and Westwater, 1962) confirmed the coherency of R-T instability during the film boiling on a horizontal surface. Water and Freon 11 (CCl3F) systems were studied to validate the predictions of Berenson correlation. A generalization of film boiling correlations on flat surface were proposed by (Klimenko, 1981). This proposed correlation was validated against different liquids including cryogens (Klimenko and Shelepen, 1982). The previously mentioned FERC “de facto” consequence analysis method has considered this correlation to estimate vapor formation source-term. For an upward facing geometrical system, according to this correlation, Nusselt number is predicted as follows:
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Nu ¼ 3:02 102 Ar1=3 Pr 3 f 1 ðbÞ;
for
Nu ¼ 1:37 103 Ar1=2 Pr 3 f 2 ðbÞ;
for
1
1
Ar < 108 Ar > 108
(2) (3)
Where,
f1
¼ 1 for b > 0:71 ¼ 0:89 b1=3 for b < 0:71
f2
¼ 1 for b > 0:5 ¼ 0:71 b1=2 for b < 0:5
C DT Cp;v mv g L3 r ðr rv Þ b ¼ p; v ; Ar ¼ c v 2l ; Pr ¼ ; hfg kv mv rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h Lc sl and characteristic length; Lc ¼ Nu ¼ kv gðrl rv Þ In spite of the fact that both Berenson and Klimenko correlations are very useful in estimating vapor generation source-term, both have fundamental limitations. For example, it is unrealistic to consider a constant film thickness over the solid surface whereas in practical case it is continuously changing owing the bubble generations and the movement of boiling the liquid. Other limitations include assumptions of constant bubble diameter and height, periodic bubble generation from the heated surface. As a result, a temporal aspect of heat flux is lost when correlations were used. Furthermore, local physical properties of the liquid and vapor were replaced by the mean properties. Computational fluid dynamics (CFD) is a great tool for overcoming these limitations and therefore utilized in this research to capture the dynamic nature of film boiling and associated wall heat flux. At the end of the XX century, CFD simulation of film boiling on a horizontal surface was pioneered by (Son and Dhir, 1997). A moving mesh method is applied to study bubble and film dynamics for water (Dhir, 2001). Water at near critical conditions was also studied by the authors to provide steady-state bubble release pattern (Son and Dhir, 1998). A sub-cooled film boiling of water on a horizontal disk was studied by (Banerjee and Dhir, 2001). To overcome the computational difficulties in simulating both low and high density ratio fluids (Juric and Tryggvason, 1998), implemented interfacial source terms in the continuity equations on a Eulerian grid. This numerical recipe was further improved by (Esmaeeli and Tryggvason, 2004a, 2004b). Volume of fluid (VOF) method (Youngs, 1982) was used by (Agarwal et al., 2004; Welch and Rachidi, 2002; Welch and Wilson, 2000) to conjugate heat transfer system for studying saturated horizontal film boiling of water on steel. The film boiling of water was studied at saturation condition and near critical point by (Tomar et al., 2008, 2005) using coupled level-set and volume of fluid (CLSVOF) method. Film boiling of water on a sphere was studied by (Yuan et al., 2008) using a non-orthogonal body fitted coordinates. The first attempts to study boiling of cryogenic liquid were attempted by (Liu et al., 2015) and (Ahammad et al., 2016) using a commercial CFD package, ANSYSFLUENT. Liu proposed different boiling routes to model three boiling phases for a pool boiling system of liquid nitrogen. In spite of a decade progress in CFD application to model film boiling, no attempt has been made to study LNG system. This study is primarily focused on CFD applications to study film boiling cryogenic liquids i.e., LNG (as pure methane) and liquid nitrogen (LN2).
3. The film boiling model The Rayleigh- Taylor (R-T) instability approach (Zuber, 1959) is
249
adopted in modeling film boiling using a two-dimensional computational domain. In this approach, bubbles generated under the liquid phase tries to invert phases due to gravity. Bubble spacing, bubble size and bubble frequency can be used to estimate the vapor formations. These parameters are governed by the hydrodynamic factors of the boiling systems, therefore CFD method is utilized. A perturbation analysis of R-T instability (Carey, 1992) suggested that for a horizontal system, the highest amount of bubble generation would occur if the distance between two nearest bubble is equal to the “most dangerous wavelength” as given by the following Equation.
ld ¼ 2p
3s ðrl rv Þg
1=2 (4)
Fig. 2 graphically present the underlying concepts of R-T approach when adopted for CFD simulation of film boiling. On the left hand side, it shows that bubble generation from the sustained vapor film which rises against the gravity force. It is assumed that the vapor film thickness is large in compared to the substrate surface roughness. The location of the surface where a bubble evolves and detaches from the film is considered as a node point. Between two nearest bubble, the point where the film thickness is lowest is considered as an anti-node point. The dotted rectangle contains four nodes at its vertices and one in the center. Among these five nodes, there are four anti-nodes as shown on the right hand side of Fig. 2. Assuming uniformity of surface, i.e., bubble generation is consistent through-out the heated surface, a domain enclosed by a node point and anti-node point would suffice for this study. Thus, a domain of ld/2 3/2 ld is used where a constant wall temperature is applied as a boundary conditions in the heated surface. The initial vapor film is a sinusoidal perturbation to induce R-T instability and expressed as Equation (5). A linear temperature gradient (Equation (6)) is applied in the initial vapor film and saturation temperature is used elsewhere in the initial temperature field.
d¼
ld 64
Ty ¼
4 þ cos
2px
ld
Twall DT·y=d Tsat
for a ¼ 1 for a ¼ 0
(5)
(6)
3.1. Simulation conditions The computation of the film boiling model is performed by a commercial CFD package ANSYS FLUENT, version 14.0. A 2D computational domain of grid size of 5.79 mm 17.37 mm for LN2 and 9.78 mm 29.33 mm for LNG were used. Grid and time step size independence test were performed for both cases and it is found that a structural grid of 64 192 segments with a time step of Dt ¼ 104 s was sufficient to capture the bubble dynamics with reasonable accuracy. Temperature dependent physical properties of LN2 (Barron, 1999) and LNG as pure methane (Rowley et al., 2003) are used. Symmetric boundary conditions were applied on the both sides of the computational domain as shown in Fig. 2. The volume of Fluid (VOF) method is used to track the vapor-liquid interface. In this method, mass, momentum and energy equations are solved once with and the additional equation for computing the volume fraction (a) of each computational cell. Vapor is considered as the primary phase. Geometric reconstruction algorithm is used to better estimate the bubble interface. To account the surface tension force which stabilizes the interfacial instability, a continuum surface force (CSF) model is used (Brackbill et al., 1992). First order
250
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Fig. 2. Two-dimensional CFD set-up for studying film boiling of cryogenic liquids.
implicit discretization is used for transient time formulations and VOF equation. Energy and momentum equations were discretized using QUICK algorithm. PRESTO procedure is used to interpolate pressure values. PISO algorithm is used to solve the corresponding governing equations. 4. Results and discussions Liquid nitrogen (LN2) system with wall superheats of 32 K, 73 K, 103 K, 153 K and 196 K and LNG boiling system with wall superheats of 43 K, 83 K, 123 K and 163 K were simulated by following
aforementioned approach. The results of these simulations are discussed as follows: 4.1. Individual bubble growth Fig. 3(a) e (d) present evolution of bubble from the initial vapor film for a LN2 system of 103 K wall superheat. In these Figures, x ¼ 5.79 mm corresponds to node point. The initialized vapor film (by Equations (5) & (6)) stabilizes at the beginning and then starts growing. Vapor generates at the interface as the liquid evaporates. As vapor generates, it increases the film thickness and moves
Fig. 3. Evolution of bubbles for a wall superheat of 103 K for LN2 system.
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
4.2. Behavior of wall heat flux The heat flux from the solid wall during the film boiling strongly dependent on the thickness of the film as also observed in other studies (Agarwal et al., 2004; Tomar et al., 2005). The behavior of heat flux can be explained by the behavior of the sustained film thickness during the bubble evolution. As vapor phase presents greater thermal resistance than the liquid phase; and thus these locations where the film thickness is smaller have higher heat flux than that of a thicker film. From Fig. 4, the heat flux decreases as the film thickness increases steadily in the beginning of bubble generation. At 0.05 s the vapor reaches at the most stable uniformly distributed condition. Vapor start pushing the bubble at this time and slowly builds up the bubble. As vapor from the film moves towards the bubble, the thickness of the film decreases on the surface because of the vapor deficiency. As a result the area weighted average heat flux increases with time. At around 0.18s (see Fig. 4 (a)) the bubble grows to its maximum size and is about to leave the film. At this point, the average film thickness on the surface is minimum, resulting maximum heat flux. Immediately after bubble detachment, the retreated vapor from bubble stem takes some time to stabilize. Because of the presence of wake behind the bubble, the interface again necks slightly and generates a smaller peak as seen in the heat flux curve seen around time from 0.21s to 0.31s. This phenomenon repeats causing large peaks representing the bubble release and troughs as points of highest film thickness, (d) in Fig. 4. The CFD simulated heat flux of this work is
compared with corresponding estimations using Berenson and Klimenko correlations in Fig. 4. It is found that the time weighted average (TWA) simulated wall heat flux is bounded by the correlations for the case of wall superheat of 103 K. Fig. 5 depicts the trend of wall heat flux for LNG at a wall superheat of 83 K. Similar to LN2, the peaks represent the liberation of bubbles which follows a period of bubble growth. The average wall heat flux is obtained from the simulations is 23 kW/m2 which is very high compared to the estimations of Klimenko and Berenson correlations.
4.3. Effect of wall superheats on bubble generation frequency and bubble diameter Fast Fourier Transform (FFT) is performed on 16,384 data points of the simulated heat flux to determine the frequency of the bubble generation. Fig. 6 depicts the results of the FFT analysis for the simulated case of LNG with a wall superheat of 43 K. The dominant frequency of the bubble generation which corresponds to the x-axis value (Fig. 6) of the highest peak has been found to be 4.272 Hz. Similar analysis has been performed for LN2 and LNG for all simulated wall superheat conditions. Fig. 7 shows bubble generation frequency for various simulated wall superheats. It is clear that the bubble formation rate is strongly
50000 Simulation
45000
Berenson
Klimenko
TWA-Simulation
40000 35000
q (W/m2)
towards the node as shown in Fig. 3(b). Further generations of vapor in the film tend to move towards to the node, withdrawing a large quantity of vapor from the film, and results in a visible bubble at Fig. 3(c). Then, the average film thickness decreases as the bubble grows and starts departing. When the buoyancy force inside the bubble cannot be compensated by surface tension and gravity, bubble detach from the film Fig. 3(d). Due to surface tension force and capillary action, some portion of vapor from the bubble stem returns to the initial film thus increases the film thickness. The average minimum film thickness is usually occurring just before the liftoff of the bubble. Individual bubble growth mechanism for LN2 and LNG are similar and therefore is not described.
251
30000 25000 20000 15000 10000 5000 0
0
0.1
0.2
0.3
0.4
0.5
0.6
Time, s Fig. 5. Surface heat flux for the film boiling of LNG at a wall superheat DT ¼ 83 K.
Fig. 4. Trend of wall heat flux for LN2 film boiling at a wall superheat of 103 K.
252
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Bubble diameter, db (m)
0.01 LN2
0.009
LNG
0.008 0.007
db = 0.0025ΔT0.2465 R² = 0.8578
0.006 0.005
db = 0.001ΔT0.2766 R² = 0.7258
0.004 0.003 0.002 0.001 0 20
60
100
140
Wall Superheat, ΔT (K)
180
Fig. 8. Bubble diameters of the simulated film boiling cases at various wall superheats.
4.4. Proposed model to predict film boiling heat transfer
Bubble Frequency, f (Hz)
Fig. 6. FFT of heat flux variation for LNG at a wall superheat of DT ¼ 43 K.
12 LN2 10
LNG f = 0.4179ΔT0.6036 R² = 0.9514
8 6
f = 0.1233ΔT0.7863 R² = 0.9933
4 2 0 20
60
100
140
180
Wall Superheat, ΔT (K) Fig. 7. Bubble generation frequency of simulated film boiling cases at various wall superheats.
dependent on the wall superheats. Within the simulated window, a model is fitted to predict the bubble generation frequency as a function of wall superheats, and it is found that, bubble frequency of LN2 follows a power of 0.8 whereas LNG follows a power of 0.6. However, the co-efficient is higher for LNG in compared to LN2 in an order of 3.4. Fig. 8 shows the dependency of bubble diameters, estimated at the time of departure from the vapor film. For a particular wall superheat, the bubble size is larger for LNG in compared to LN2. A power model is fitted (shown on Fig. 8) for both LN2 and LNG to predict the bubble diameter as a function of wall superheats. For LNG the exponent of the model has found to be 0.25 and for LN2, its value is 0.28. However, the co-efficient for LNG is 2.5 times higher than that of LN2.
Fig. 9(a) compares the time weighted average (TWA) simulated heat fluxes with the heat fluxes estimated by using Berenson and Klimenko correlations, for a LN2 film boiling system. Similarly, Fig. 9(b) compares that of for LNG film boiling system. For the case of LN2, at a lower wall superheat, the simulated heat flux is closer to the estimation of Klimenko correlations and at higher wall superheat; the simulated heat flux is closer to estimation of Berenson correlations. The trend of the simulated results, depict that a logarithmic slope is much higher than the slope of both correlations. For the case of LNG, the TWA simulated heat fluxes are much higher than the estimations from Berenson and Klimenko correlations. It is important to mention that, the simulated results are not subjected to the assumptions such as constant film thickness, constant bubble diameters, absence of dynamic movement of vapor and liquid, usage of mean fluid properties etc. Thus, the simulated heat fluxes are more realistic than the aforementioned correlations. Film boiling simulation for a single wall superheat takes approximately 800 min of computer time (Processor: 3.33 GHz; RAM: 32 GB) to simulate two consecutive bubble detachments. Therefore, performing simulations for a large range of wall superheats might not be always possible when the time does not allow. Therefore, a first principle model based on simple heat balance is used to predict the wall heat flux for non-simulated wall superheats. The proposed model is based on the heat balance as the amount of heat transferred from the wall is approximately equal or proportional to the heat taken by the bubbles that were formed in the film boiling regime. This can be mathematically expressed as Equation (7).
qfVb $f $m$rv $hfg
(7)
For a two-dimensional CFD simulation, the above equation can be modified as follows:
q ¼ C·
pd2b 4
·f ·m·rv ·hfg
(8)
where, it is assumed that C is a dimensionless constant equals to 1.0. For some cases C might not be 1.0. For an example, if the local physical property data such as r, m and hfg strongly varies with the temperature field, C ¼ 1.0 may not be the proper representation. Furthermore, the bubble diameter db will also vary with the time of simulation and wall superheats, etc. The departure diameter of the first bubble for a particular wall superheat might not be equal to
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
253
q (W/m2)
50000
5000
q (W/m2)
5000
TWA Simulation Klimenko Berenson Proposed model
500 20
Wall superheat, ΔT (K)
Proposed model TWA Simulation Klimenko Berenson
500
200
20
(a) LN2
Wall superheat, ΔT (K) (b) LNG
200
Fig. 9. CFD based correlation to predict heat flux of (a) LN2 and (b) LNG film boiling.
that of the second bubble. The shape of the bubble can be different from spherical shape. Wake formation in the liquid and liquid depth may also affect the size and shape of the bubbles. Therefore, adjustment of the constant C might be needed to accurately predict the heat flux of the non-simulated cases. The parameters in Equation (8) can be estimated using the simulation results. Fourier transformation is employed to determine the bubble generation frequency f from the transient wall heat flux for different wall superheat as depicted in Figs. 6 and 7. The departure diameter of the simulated bubble is considered to estimate the volume of the bubble (for 3D simulation) or area of the bubble (for 2D simulation). Therefore, by replacing the term, f and db using the fitted model as shown in Figs. 7 and 8 and replacing the parameter m, i.e., number of nodes per unit area, with 1/L2c , the proposed model for LN2 and LNG are as follows:
q ¼ 8:02 109
q ¼ 1:7 107
rl rv rv Lv ·DT 1:34 ðfor LN2 Þ s
(9)
Acknowledgments This paper was made possible by a National Priority Research Project (NPRP) award [NPRP 6-425-2-172] from the Qatar National Research Fund (a member of The Qatar Foundation). Nomenclature
s r g
l
rl rv rv Lv ·DT 1:1 ðfor LNGÞ s
the heat flux in the film boiling regime. The proposed model can be used along with in 1-D standard heat transfer model to determine the liquid phase resistance, and thus, enabling more accurate estimation of vapor formation source-term.
(10)
Fig. 9 depicts the fitness of the correlation for the proposed correlation for the two-dimensional simulations performed in this study. It is shown that first principle based correlation predicts the simulated heat fluxes to an excellent agreement. 5. Conclusions Film boiling of cryogenic liquids, i.e., LN2 and LNG, is simulated by adopting Rayleigh-Taylor (R-T) approach in a commercial CFD environment (ANSYS Fluent). It is found that the dynamics of film thickness affects the wall heat flux and therefore the frequency of bubbles evolution and bubbles diameter. The bubble generation frequency is estimated from the transient wall heat flux by decomposing via FFT and the most dominant frequency is considered as the average bubble frequency for that particular wall superheat. In the simulated window, LN2 and LNG have shown different trends in bubble generation frequency and bubble departure diameter resulting from increased wall superheats. The simulated film boiling heat fluxes for LN2 were higher in comparison to the estimates by Klimenko correlations but it underestimates when compared to Berenson correlation. On the other hand, the LNG simulations over-estimates the wall heat fluxes when compared against both Berenson and Klimenko correlations. Finally, based on the trends of bubble generation frequency and bubble diameters, a first principle model is applied for correlating
hfg k
m p Cp T DT t Dt Lc
a h
d x, y q Vb f n db m
Surface tension, N/m Density, kg/m3 Acceleration of gravity, 9.81 m/s2 Wavelength, m Heat of vaporization, J/kg Thermal conductivity, W/m.K Dynamic viscosity, Pa.s Pressure, Pa Heat capacity, J/kg.K Temperature, K Wall superheat, K Time, s Time step, s Characteristic length, m Volume fraction; vapor volume fraction Average local heat transfer co-efficient, W/m2.K Thickness of the initial vapor film, m Axis of the simulated plane, m Wall heat flux, W/m2 Bubble volume, m3 Frequency of bubble generation, Hz Number of nodes per square area Diameter of the bubble just before detachment from the film, m Number of nodes per square meter
Subscripts v Vapor phase l Liquid phase wall, sat Indication of wall condition and saturation condition b Bubble y Y axis direction
254
M. Ahammad et al. / Journal of Loss Prevention in the Process Industries 44 (2016) 247e254
Dimensionless numbers f1, f2, b Dimensionless number Nu Nusselt number Ar Archimedes number Pr Prandtl number References Agarwal, D.K., Welch, S.W.J., Biswas, G., Durst, F., 2004. Planar simulation of bubble growth in film boiling in near-critical water using a variant of the VOF method. J. Heat. Transf. 126, 329. http://dx.doi.org/10.1115/1.1737779. Ahammad, M., Liu, Y., Olewski, T., Vechot, L., Mannan, M.S., 2016. Application of computational fluid dynamics (CFD) in simulating film boiling of cryogens. Ind. Eng. Chem. Res. 55, 7548e7557. http://dx.doi.org/10.1021/acs.iecr.6b01013. Banerjee, D., Dhir, V.K., 2001. Study of subcooled film boiling on a horizontal disc: Part IdAnalysis. J. Heat. Transf. 123, 271e284. http://dx.doi.org/10.1115/ 1.1345889. Barron, R.F., 1999. Cryogenic Heat Transfer. Taylor & Francis, Philadelphia, PA, USA; London, UK. Berenson, P.J., 1961. Film-boiling heat transfer from a horizontal surface. J. Heat. Transf. 351e356. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335e354. Bubbico, R., Salzano, E., 2009. Acoustic analysis of blast waves produced by rapid phase transition of LNG released on water. Saf. Sci. 47, 515e521. http:// dx.doi.org/10.1016/j.ssci.2008.07.033. Carey, V.P., 1992. Liquid-vapor Phase Change Phenomena: an Introduction to the Thermophysics of Vaporization and Condensation Process in Heat Transfer Equipment. Hemisphere Publishing Corp, Washington , D.C. Dhir, V.K., 2001. Numerical simulations of pool-boiling heat transfer. AIChE J. 47, 813e834. http://dx.doi.org/10.1002/aic.690470407. Esmaeeli, A., Tryggvason, G., 2004a. Computations of film boiling. Part I: numerical method. Int. J. Heat. Mass Transf. 47, 5451e5461. http://dx.doi.org/10.1016/ j.ijheatmasstransfer.2004.07.027. Esmaeeli, A., Tryggvason, G., 2004b. Computations of film boiling. Part II: multimode film boiling. Int. J. Heat. Mass Transf. 47, 5463e5476. http://dx.doi.org/ 10.1016/j.ijheatmasstransfer.2004.07.028. FERC, 2016. Approved North American LNG Import and Export Terminals [WWW Document]. http://www.ferc.gov/industries/gas/indus-act/lng/lng-approved. pdf (accessed 04.21.16.). FERC, 2004. Notice of Availability of Detailed Computations for the Consequence Assessment Methods for Incidents Involving Releases from Liquefied Natural Gas Carriers. Havens, J., Spicer, T., 2007. United States regulations for siting LNG terminals: problems and potential. J. Hazard. Mater 140, 439e443. Holster, E.R., Westwater, J.W., 1962. Film boiling on a horizontal plate. Am. Rocket Soc. J. 553e558. Johnson, D.W., Cornwell, J.B., 2007. Modeling the release, spreading, and burning of
LNG, LPG, and gasoline on water. J. Hazard. Mater 140, 535e540. http:// dx.doi.org/10.1016/j.jhazmat.2006.10.022. Juric, D., Tryggvason, G., 1998. Computations of boiling flows. Int. J. Multiph. Flow. 24, 387e410. Klimenko, V.V., 1981. Film boiling on a horizontal plate - new correlation. Int. J. Heat. Mass Transf. 24, 69e79. http://dx.doi.org/10.1016/0017-9310(81)90094-6. Klimenko, V.V., Shelepen, A.G., 1982. Film boiling on horizontal plate - a supplementary communication. Int. J. Heat. Mass Transf. 25, 1611e1613. chot, L.N., 2015. Modeling of a cryogenic liquid pool boiling by Liu, Y., Olewski, T., Ve CFD simulation. J. Loss Prev. Process Ind. 35, 125e134. NFPA 59A, 2013. Standard for the Production, Storage, and Handling of Liquefied Natural Gas (LNG). National Fire Protection Association. Reid, R.C., 1980. Boiling of LNG on Typical Dike Floor Materials. Rowley, R.L., Wilding, W.V., Oscarson, J.L., Zundel, J.L., Marshall, N.A., Daubert, T.L., Danner, R.P., 2003. DIPPR Data Compilation of Pure Compound Properties. Design Institute for Physical Property Data/AIChE. Son, G., Dhir, V.K., 1998. Numerical simulation of film boiling near critical pressures with a level set method. J. Heat. Transf. 120, 183e192. Son, G., Dhir, V.K., 1997. Numerical simulation of saturated film boiling on a horizontal surface. J. Heat. Transf. 119, 525e533. Tomar, G., Biswas, G., Sharma, A., Agrawal, A., 2005. Numerical simulation of bubble growth in film boiling using a coupled level-set and volume-of-fluid method. Phys. Fluids 17. http://dx.doi.org/10.1063/1.2136357. Tomar, G., Biswas, G., Sharma, A., Welch, S.W.J., 2008. Multimode analysis of bubble growth in saturated film boiling. Phys. Fluids 20, 1e8. http://dx.doi.org/10.1063/ 1.2976764. United States Government Accountability Office, 2007. Report GAO-07-316maritime Security: Public Safety Consequences of a Terrorist Attack on a Tanker Carrying Liquefied Natural Gas Need Clarification Highlights. chot, L., Olewski, T., Osorio, C., Basha, O., Liu, Y., Mannan, S., 2013. Laboratory scale Ve analysis of the influence of different heat transfer mechanisms on liquid nitrogen vaporization rate. J. Loss Prev. Process Ind. 26, 398e409. http:// dx.doi.org/10.1016/j.jlp.2012.07.019. Webber, D.M., Gant, S.E., Ivings, M.J., Jagger, S.F., 2010. LNG Source Term Models for Hazard Analysis: a Review of the State-of-the-art and an Approach to Model Assessment (No. Research Report: RR789). HSE Books, Buxton, Derbyshire, UK. Welch, S.W.J., Rachidi, T., 2002. Numerical computation of film boiling including conjugate heat transfer. Numer. Heat. Transf. B 42, 35e53. Welch, S.W.J., Wilson, J., 2000. A volume of fluid based method for fluid flows with phase change. J. Comput. Phys. 160, 662e682. http://dx.doi.org/10.1006/ jcph.2000.6481. Youngs, D.L., 1982. Time dependent multi-material flow with large fluid distortion. In: Morton, K.W., Baines, M.J. (Eds.), Numerical Methods in Fluid Dynamics. Academic Press, pp. 273e285. Yuan, M.H., Yang, Y.H., Li, T.S., Hu, Z.H., 2008. Numerical simulation of film boiling on a sphere with a volume of fluid interface tracking method. Int. J. Heat. Mass Transf. 51, 1646e1657. http://dx.doi.org/10.1016/ j.ijheatmasstransfer.2007.07.037. Zuber, N., 1959. Hydrodynamic Aspects of Boiling Heat Transfer. University of California, Los Angeles.