Application of Computational Fluid Dynamics in ... - ACS Publications

Article pubs.acs.org/IECR

Application of Computational Fluid Dynamics in Simulating Film Boiling of Cryogens Monir Ahammad,† Yi Liu,† Tomasz Olewski,‡ Luc N. Véchot,‡ and M. Sam Mannan*,† †

Mary Kay O’Connor Process Safety Center, Artie McFerrin Department of Chemical Engineering, Texas A&M University, College Station, Texas 77843-3122, United States ‡ Mary Kay O’Connor Process Safety CenterQatar, Chemical Engineering Program, Texas A&M University at Qatar, Doha, Qatar ABSTRACT: This paper presents a computational fluid dynamic (CFD) model, simulating film boiling based on Rayleigh−Taylor (R-T) instability, using the volume of fluid (VOF) method to track the liquid/vapor interface. Film boiling of cryogenic liquids (e.g., LNG and liquid nitrogen) is simulated to estimate the vapor generation rate during an accidental spill. The simulated heat fluxes were compared with heat fluxes obtained from Berenson and Klimenko correlations. The effects of wall superheats on the bubble generation frequency were studied. This study helps researchers to understand the physics of film boiling that are useful during the risk assessment of a cryogenic spill scenario. For example, it was found that the bubble released from the node and the antinode points between the consecutive bubble generations cycles do not follow the alternating nature under the realistic film boiling conditions. Therefore, empirical expressions assuming alternating bubble generation might be unsuitable for cryogenic vaporization source term estimation. regime is called film boiling. As the solid surface cools down, the vapor cannot be sustained as a continuous film, thus breaking to intermittent film, and liquid comes in contact with solid intermittently. This regime is called the transient boiling regime. For cryogenic boiling, the transient regime is very short compared to the film boiling. Therefore, cumulative vapor generation in this regime is smaller than that for the other boiling regimes. Further cooling of the heated surface, i.e., the ground, will cause nucleate boiling. In this regime, bubbles are generated from the nucleation sites of the substrate. Substrate geometrical properties, such as surface roughness, play a significant role in determining the vapor generation in this regime. Nucleate boiling occurs at the later stage of the spill, when the pool is developed and the temperature gradient is small. For an accurate estimation of LNG spill consequence, the change of boiling regimes during vaporization should be considered. However, as a part of a wider study to incorporate the changes in boiling regimes during vaporization, this study only focuses on the film boiling part. It is expected that the highest boiling rate would occur in the film boiling regime owing to the largest temperature gradient in the substrate. Zuber6 introduced the concept of Rayleigh−Taylor (R-T) hydrodynamic instability to model film boiling on a horizontal surface. He proposed that the nearest distance of evolving vapor bubbles should be bounded by Taylor “critical” and “most dangerous” wavelengths as defined in eqs 1 and 2,

1. INTRODUCTION With tremendous recent growth of LNG (liquified natural gas) liquefaction, storage, transportation, and regasification all over the world, the need for better risk assessment procedures associated with these operations was emphasized by the experts and government agencies.1,2 Standards such as NFPA 59A, therefore, recommend the calculation of an LNG spill consequence, to be a realistic estimation via validated models.3 Upon a loss-of-containment, LNG will form a cryogenic liquid pool on land or water depending on the source of release. The severity of the consequences such as fire and explosion depends on pool spreading and dispersion modeling, which in turns depend on the source terms, such as the vaporization rate of the pool. For cryogenic liquids, the vaporization rate is governed by the heat flux provided to the pool. The overall heat flux is dominated, at least for the beginning of the spill, by conductive heat transfer from the substrate.4 Most of the existing models utilize a 1-D heat conduction ideal model to estimate LNG evaporation rate when the pool is fully developed. However, in the context of consequence estimation due to LNG release, NFPA 59A prescribes the analysis of the first 10 min developing scenario to determine the worst-case condition.3 Moreover, the assumption of perfect contact between the liquid and solid substrate in 1-D ideal models is unrealistic5 and does not address the presence of different boiling regimes in estimating the heat transfer rate. In the early stage of the spill, high temperature gradient is present between the boiling liquid and the substrate, resulting in a continuous sustained vapor film. The boiling liquid does not come in contact with the solid substrate, and rather takes heat from the vapor film present in-between. This boiling © XXXX American Chemical Society

Received: March 13, 2016 Revised: June 18, 2016 Accepted: June 23, 2016

A

DOI: 10.1021/acs.iecr.6b01013 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research respectively. He also predicted the minimum heat flux of film boiling (Leidenfrost point) by assuming that two bubbles evolve per cycle from a square cell where the base of the cell is the distance between the bubbles. In a horizontal cocurrent system, where the dense phase lays over the less dense phase, the interface may become unstable in the presence of disturbance. This phenomenon is referred to as the R-T instability. In order to grow a bubble, Carey’s7 perturbation analysis suggests that the initial perturbation wavelength should be greater than a corresponding “critical” wavelength, as defined by ⎡ ⎤1/2 σ ⎥ λc = 2π ⎢ ⎢⎣ (ρl − ρv )g ⎥⎦

β=

(2)

8

Berenson formulated an empirical expression (eq 3) to calculate the heat transfer coefficient, for saturated film boiling on horizontal surfaces, based on the same analysis. In deriving this expression, it was assumed that bubbles are spaced on a square grid of thin vapor film by a distance of “most dangerous” wavelength. As a result, two bubbles can generate per λD2 area of heated surface at any particular moment. This twodimensional study was improved for three dimensions by Sernas et al.9 They showed that the three-dimensional Taylor wavelength, λD3, is √2 times larger than λD and the release was four bubbles per cycle from an λD32 area. ⎡ ρ (ρ − ρ )ghfg ⎤1/4 ⎡ ⎤3/8 σ v l v ⎥ ⎢ ⎥ Nu = 0.425⎢ ⎢⎣ ⎥⎦ ⎢⎣ g (ρl − ρv ) ⎥⎦ k vμv ΔT

(3)

10

Holster and Westwater experimentally confirmed that the film boiling from a horizontal surface follows R-T instability. Their experimental studies have found that the film boiling heat fluxes for water and Freon-11 (CCl3F) are in agreement with the prediction of the Berenson correlation (eq 3). Klimenko11 attempted to generalize the film boiling correlations on horizontal flat plates for different liquids including cryogens, therefore extending the experimental database for film boiling.12 His correlation predicts the Nusselt number during film boiling, in a geometrical system consisting of an upward facing horizontal surface, as Nu = 3.02 × 10−2Ar1/3Pr1/3f1 (β )

for Ar < 108

(4)

Nu = 1.37 × 10−3Ar1/2Pr1/3f2 (β )

for Ar > 108

(5)

, Ar =

gLc 3ρv (ρl − ρv ) μv

2

, Pr =

Cp,vμv

h ̅ Lc and characteristic length, Lc = kv

kv

,

σl g (ρl − ρv )

The fundamental limitations of the existing empirical expressions, e.g., the Berenson and the Klimenko correlations, are unrealistic assumptions such as a constant film thickness, proportional relationship between the bubbles diameter and height, periodic bubble generations, and periodic and alternating nature of bubbles liberating from the node and antinode points of the heated surface. Hence, such correlations cannot predict the temporal variation of the heat flux. Additionally, these expressions (i.e., eqs 3−5) do not use local physical properties of the liquid and vapor but instead are replaced by the mean properties. To overcome these limitations, the computational fluid dynamics (CFD) tool can be employed to further clarify the physics, and to accurately estimate the dynamic nature of film boiling and the associated heat flux reliably. Numerical simulation of horizontal film boiling was pioneered by Son and Dhir.13 The authors studied bubble and film dynamics for water boiling using the moving-mesh method and also presented a combined scheme for nucleate and film boiling.14 Further considerations for near critical conditions of water in an axisymmetric horizontal film boiling have been studied to provide a steady-state bubble release pattern.15 Panzarella et al.16 modeled film boiling of water by using a lubrication approximation and thereby solving a strongly nonlinear evolution equation. Banerjee17 simulated subcooled film boiling of water on a horizontal disk. Juric18 used added interfacial source terms in the continuity equations on a Eulerian grid to simulate horizontal film boiling of low density ratio fluid to high density ratio fluid. This numerical method is further improved by Esmaeeli19,20 by elimination of iterative algorithm. Welch21 used Youngs’22 volume of fluid (VOF) method to simulate saturated horizontal film boiling and conjugate heat transfer. Using this method, Welch and Rachidi23 simulated film boiling of water in contact with steel. Yuan et al.24 simulated the film boiling of water on a sphere on a nonorthogonal body fitted coordinates. Agarwal et al.25 simulated film boiling of water at 373 °C, 219 bar using a variant of the VOF method to study the unsteady bubble release patterns, transport coefficients, and influence of fluid properties. Tomar et al.26,27 and Hens et al.28 studied water and refrigerant R134a at near and far critical pressures using coupled level-set and volume of fluid (CLSVOF) method. Welch and Biswas29 and Tomar et al.30 investigated the effect electrical potential on heat transfer by performing direct simulation of film boiling. Liu et al.31 studied pool boiling of liquid nitrogen using commercial computational fluid dynamics package ANSYS Fluent. As mentioned above, many studies address film boiling numerical approaches, bubble generation dynamics, and associated heat transfer, but no notable attempts have been taken to simulate cryogenic fluid boiling particularly LNG for the application of vaporization source term estimation. The most studied boiling systems are water and refrigerants at near-critical pressures. Hence, this study addresses film boiling of cryogenic systems (i.e., LNG and LN2) at atmospheric conditions, and is particularly useful for reliable estimation of cryogenic boiling, e.g., LNG source term modeling.

(1)

3 λc

hfg

Nu =

If the length of interface in the horizontal direction is less than λc, the interface is stable because a perturbation of wavelength less than λc cannot arise; therefore, no film boiling will be observed. Carey’s study further suggests that the initial disturbance wavelength that corresponds to maximum vapor rise from the vapor film is referred as the “most dangerous” wavelength λD. λD =

Cp,v ΔT

where f1 = 1 for β > 0.71 = 0.89β −1/3 for β < 0.71 f2 = 1 for β > 0.5 = 0.71β −1/2 for β < 0.5 B


Article

Industrial & Engineering Chemistry Research

Figure 1. (a) Actual vapor−liquid interface. (b) Interface reconstructed using PLIC. (c) Interface using donor−acceptor scheme.

∂α + v ⃗ · ∇α = 0 ∂t

The remainder of this paper is organized as follows: Section 2 describes the mathematical formulations of the CFD model that is used to simulate a two-dimensional film boiling system. Section 3 describes the details of simulation setup in ANSYS Fluent 14.0, a commercial CFD software program, to solve the mathematical formulations presented in section 2. Grid and time step sensitivity analysis is also presented in this section. Film boiling of liquid nitrogen (LN2) and LNG (as 100% liquid methane) is simulated using this setup. Section 4 discusses the results of the film boiling simulations which mainly focus on interface morphology, behavior of heat flux, effects of wall superheat, and the assumptions of alternating bubble generations that were considered in other CFD studies.13 In this section, the simulated film boiling heat fluxes for LN2 and LNG were also compared with the estimated heat flux by using the Klimenko and the Berenson correlations.

(6)

Several researchers applied this method for simulating film boiling. Hardt and Wondra33 proposed a method for applying VOF to perform film boiling simulations and droplet evaporation. Kunkelmann34 implemented VOF solver in the open-source CFD package “OpenFoam” to solve incompressible two-phase problems. Kunugi35 performed a comprehensive review of the latest simulations. Further information on VOF use can be found in Kunugi’s work.35 In ANSYS Fluent, the tracking is accomplished by the solution of continuity equation for the volume fraction of one phase (page 47536). For the qth phase, this equation is ⎤ 1⎡∂ ⎢⎣ (αqρq ) + ∇·(αqρq vq⃗ )⎥⎦ = Sαq + ρq ∂t

n

∑ (ṁ pq − ṁqp) p=1

(7)

where ṁ qp is the mass transfer rate from phase q to phase p, and Sαq is the mass source term. The primary fluid in this study is considered as the gas phase. Thus, the above-mentioned equation is solved for the liquid phase only. The vapor phase volume of fraction is calculated using the constraint

2. MATHEMATICAL FORMULATIONS Commercial CFD software ANSYS Fluent was used to implement and solve the formulation of the film boiling model as described in the following. 2a. Interface Tracking Using VOF Method and Calculation of Curvature. Volume of fluid (VOF) method is used to capture the vapor−liquid interface in a fixed Eulerian mesh. A single set of momentum equation as shown in next sections is solved to determine the volume fraction (α) in each computational cell. In a control volume (i.e., a single cell of the solution domain), the summation of the volume fractions of liquid and vapor phase is equal to unity. For example, boiling of one component two-phase system, each cell in the solution domain is either filled with liquid, vapor, or a mixture of liquid and vapor phases (say 50% liquid and 50% vapor); where in any cases, the summation of volume fraction will be 1. Thus, for a multiphase cell, the cell or node property value or field variable value represents the volume-averaged value. If qth fluid’s volume fraction in a cell is denoted by αq, then the following apply: if αq = 0, the cell is empty of qth fluid; if αq = 1, the cell is full of qth fluid; if 0 < αq < 1, there is an interface in the cell between qth fluid and at least one other fluid. The VOF method is based on the conservation of α with respect to time and space as expressed in the following equation. This fundamental idea was originated by Hirt and Nichols.32

n

∑ αq = 1 q=1

(8)

On the basis of the need for interface reconstruction, VOF methods are categorized into two categories, i.e., those that require vapor/liquid interface reconstruction and those that do not. SOLA-VOF (based on donor−acceptor method),37 FCTVOF,38 and CICSAM (compressive interface capturing scheme for arbitrary meshes)39 do not require interface reconstruction. SLIC (simple line interface calculation)40 and PLIC (piecewise linear interface calculation)22 methods require reconstruction of interface. Because the PLIC method provides a more accurate solution than the other methods, a more realistic representation of the actual interface is possible via this method (Figure 1). Thus, it was chosen for use as the interface reconstruction method in this study. 2b. Material Properties. The material properties in the transport equations are computed as follows.

C

ρ = (1 − α)ρv + αρl

(9)

μ = (1 − α)μv + αμ l

(10) DOI: 10.1021/acs.iecr.6b01013 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Figure 2. R-T instability approach for film boiling simulation using CFD.

k = (1 − α)k v + αkl

momentum and energy discretization, the QUICK algorithm is used.36 2f. Discretization of VOF Equation. First order implicit scheme is used for time discretization, therefore, ANSYS Fluent’s standard finite difference interpolation scheme; QUICK is used to obtain the face fluxes for all cells.

(11)

2c. Governing Equations. A single momentum equation is solved for the entire computational domain. The calculated velocity field is shared among the vapor and liquid phases. The momentum equation shown in eq 12 is dependent on the volume fractions of both phases via average material properties.

αqn + 1ρqn + 1 − αqnρqn

∂(ρv ⃗) + ∇·(ρvv⃗ ⃗) ∂t

Δt T

= −∇p + ∇·[μ(∇v ⃗ + ∇v ⃗ )] + ρg ⃗ + F ⃗

(12)

= [Sαq +

(14)

n

∑q = 1 αqρq Eq n

(15)

Here, Eq is based on specific heat of the qth phase and the temperature of the computational cell. 2d. Continuum Surface Force Model. Surface tension in the interface creates a jump in density and energy across the interface. The continuum surface force (CSF) model developed by Brackbill et al.41 is used to capture this jump conditions via addition of surface force as source term in the momentum equation. Fvol = σqp 1 2

ρκq∇αq (ρp + ρq )

(17)

3. SIMULATION SETUP During film boiling, a sustained film of vapor is always present between the solid substrate and boiling liquid. Phase change (liquid to vapor) occurs at the vapor−liquid interface. It is assumed that the surface roughness of the substrate is much smaller than the film thickness; hence, there is no significant effect on film boiling. It is also assumed that the bubble generation from the vapor film follows a regular pattern. Therefore, heat transfer from a large area can be estimated by repeating the simulated domain. From Figure 2, one bubble evolves in each cycle from a square cell of area λd22, where λd2 is the “most dangerous” Taylor wavelength. The point at which bubbles are growing is referred to as the node, and the valley of two adjacent bubbles is called the antinode. To consider the symmetry of the bubbles over the entire hot surface, a horizontal length of λd2/2 needs to be simulated as shown in Figure 3. To capture the bubble dynamics properly, the height of the two-dimensional simulation domain is considered as 3 times the width of the domain.42 The bottom of the domain is considered as a hot surface at constant temperature whereas the top of the domain is considered as vapor outlet. The simulation is initialized with a linear temperature profile in the vapor film at the bottom of the computational domain as shown in the following equations. The initialized vapor film takes the form of a sinusoidal perturbation that induces R-T instability in the computational domain, thus enabling the simulation of film boiling phenomena. The simulation setup is

Enthalpy (E) and temperature (T) of each cell is considered as a mass-averaged quantity given by the following equation:

∑q = 1 αqρq

∑ (ṁ pq − ṁqp)]V

Iterative solution of a standard scalar transport equation for the secondary-phase volume fractions at each time step determines the volume fraction values at the current time step.

(13)

The energy equation is also dependent on the volume fractions of both phases via the material properties, and is shared among the phases within a computational cell.

E=

f

p=1

is

∂ (ρE) + ∇·(v ⃗(ρE + p)) = ∇·(k∇T ) ∂t

∑ (ρqn+ 1vnf + 1αqn,+f 1)

n

For the incompressible flow, the mass conservation equation ∇·v ⃗ = 0

V+

(16)

2e. Discretization. First order implicit discretization is used for transient formulation of time. PISO (pressure implicit with splitting of operator) algorithm is used to solve the governing equations. PRESTO, pressure discretization scheme, is selected because of its effectiveness in multiphase system. For D


Article


7.21 mm. Therefore, the deviation in height for these two grids is about 100%. However, the height of the bubble for 96 × 288 grids is almost same as 64 × 192 grids, 7.27 mm; therefore, the deviation in height is less than 1%. From Figure 4, a deviation in the bubble radius is also observed. Usually, higher mesh resolution results in greater bubble radius. However, the 32 × 96 mesh shows a greater diameter because of the fact that at this iteration level the bubble was still growing whereas, for other grid resolutions, the bubbles were about to leave the film. Comparison of the interface of the 32 × 96 grid with the 64 × 192 grid at a height of 1 mm indicates that the bubble radius 32 × 96 grid is 155% that of 64 × 192. Similarly, in a comparison of the results of 64 × 192 grids with 96 × 288 grids at a height of 6 mm, it is found that the bubble radius of 96 × 288 is 18% greater than that of 64 × 192. Although there is a difference between the interface evolution for different grids, considering the benefits and the cost of computer run time, 64 × 192 grid resolution is used as a working grid in this study. 3b. Time Step Sensitivity Analysis. The time step has been chosen to satisfy the Courant−Friedrichs−Lewy (CFL) condition for the convergence of the simulations. Figure 5

Figure 3. Setup of 2D film boiling simulations.

benchmarked by reproducing the film boiling case described in section 4.4 of Gibou et al.43 δ=

⎛ 2πx ⎞⎞ λd2 ⎛ ⎜⎜4 + cos⎜ ⎟⎟⎟ 64 ⎝ ⎝ λd2 ⎠⎠

⎧Twall − ΔT ·y/δ for α = 1 Ty = ⎨ Tsat for α = 0 ⎩

(18)

⎪

⎪

(19)

3a. Grid Sensitivity Analysis. Three mesh resolutions of size 32 × 96, 64 × 192, and 96 × 288 were used to assess the sensitivity of the film boiling model to the mesh resolution. Figure 4 shows the bubble interfaces for three different mesh

Figure 5. Time convergence study showing the bubble interface at 0.67 s using 64 × 192 meshes for different time steps.

depicts the difference between the bubble evolution for a time step of Δt = 0.0001 s and Δt = 0.00001 s. The maximum difference between the heights is less than 0.1%. Therefore, a working time step (Δt = 0.0001 s) is used for all the simulations. It takes about 800 min to simulate 2 bubble generations in the ANSYS Fluent package installed in a terminal server with Intel I Xeon I CPU, dual 3.33 GHz processors and 32.0 GB of installed RAM.

4. RESULTS AND DISCUSSION The temperature dependent physical properties of both liquid and vapor phases of nitrogen are estimated on the basis of the correlations provided in the DIPPR44 database and the physical properties of LNG, as pure methane is collected from Barron.45 4a. Interface Morphology and Behavior of Heat Flux During LN2 Film Boiling. Figure 6 presents the evolution of a bubble from the initial sustained film at a wall superheat of 32K. In the beginning, the node point has the highest thickness of the vapor film. The average film thickness is also highest. Therefore, the heat flux as shown in Figure 7 is the lowest. As it grows, the generated vapors move toward the node, and

Figure 4. Grid sensitivity analysis of film boiling simulation of LCH4.

resolutions at a time of 0.18 s which corresponds to 90 000 iterations for ΔT = 43 K. It is observed from the figure that the bubble evolution speed is greater for the denser grid. The difference between the height and diameter of the bubbles for the grid size of 32 × 96 with the grid size of 64 × 192 is not very significant. The bubble heights at 0.18 s of flow time for 32 × 96 grids and 64 × 192 grids are correspondingly 3.6 mm and E


Article


Figure 6. Evolution of bubbles during LN2 film boiling at wall superheat of 32 K.

therefore, it also draws the vapor from its sides resulting in decrease of the average film thickness. As shown in Figure 6c, when the bubble is about to detach, the average vapor film thickness is the lowest; therefore, the wall heat flux reaches its peak which is shown in Figure 7. After the detachment of the bubble, some vapor from the vertical stem returns to the film in contact with the wall. The average film thickness therefore increases again, resulting in a drop in the wall heat flux. This process is repeated in between the node and antinode. Figure 7 depicts the wall heat flux for bubble generation over 7 cycles. The time weighted average wall heat flux due to film boiling simulation is compared with the Berenson and Klimenko correlations as depicted in Figure 7. The Berenson correlation was validated for pentane and other high boiling point fluids46 whereas the Klimenko correlation is validated against the liquid nitrogen.11 It is clear that the average heat flux from film boiling

Figure 7. Simulated wall heat flux for film boiling of LN2, ΔT = 32 K.

Figure 8. Velocity vectors for LN2 a wall superheat of 73 K: (a) before bubble release and (b) after bubble release. F


Article


Furthermore, Figure 9 shows an interesting pattern on the wall. A wavy temperature profile exists at the boundary condition. This might be connected with the movement of a minimum thickness point. The local highest temperature gradient can be observed at the locations of the local minimum thickness point. Furthermore, it is clear that the isotherm follows the shape of interface. 4b. LNG Film Boiling. Figure 10a−c shows the bubble evolution and streamlines of vortex formation. Similar to LN2 film boiling, LNG vapor in the film moves toward the node and forms the bubble. When the bubble leaves the film, due to wake formation in the liquid phase, the bottom interface of the bubble inverts to take the shape of the cup. The streamlines clearly depict the locations of heat transfer in the vapor−liquid interface. The maximum amount of heat transfer occurs at the minimum thickness point. Generated vapor moves toward the bubble, and therefore, the bubble continuously pushes the liquid. As a result there is wake formation in the liquid phase. It is clear from Figure 10 that the vortex forms at the location inside the bubble with a lowest curvature. After the bubble is released, the formation of vortex is at the location of the retracted film. Figure 11 depicts the area weighted average wall heat flux as a function of time. Due to the heat transfer, vapor is generated at the vapor liquid interface, and thus, the film thickness increases. Therefore, the average heat flux decreases. As the vapor production continues, vapor from the film moves toward the node. As a result, the bubble starts growing at the node, and because of the vapor deficiency in the film, the average film thickness decreases, and the heat flux increases. When the bubble leaves the film, a portion of vapor retracts to the film causing an increase of average film thickness; therefore, heat flux at the wall suddenly drops. Each peak in Figure 11 represents the formation of one bubble. At the beginning of the simulation, the liquid was quiescent. However, after one or two bubble formations, there is a significant amount of wake or churning motion in the liquid phase, which further affects the film thickness. The complicated interactions of the liquid phase, film thickness, and interface movements represent a realistic film boiling scenario. Therefore, the heat flux peaks as shown in Figure 11 are not entirely periodic. 4c. Effect of Wall Superheat. Figure 12 depicts the effect of wall superheat on the bubble generation frequency during the film boiling of LN2. For the wall superheat of 32 K (Figure 12a−g), the first bubble generates at 0.61 s, and the subsequent bubble generates at 1.19, 1.76, 2.24, 2.53, and 3.17 s. As the wall superheat is increased to 73 K (Figure 12h−m), the bubble generation frequency increased. In this case, the first bubble was generated at 0.25 s, and subsequent bubbles were generated at 0.55, 0.72, 1.0, and 1.25 s. Further increase of the wall superheat to 103 K (Figure 12n−t) depicts that the bubble generation frequency was also increased. In this case, the first bubbles were generated at 0.19 s, and subsequently, the bubble generated at 0.37, 0.62, 0.83, 0.99, and 1.16 s. Therefore, the increase of wall superheat increases the amount of vapor generation as the frequency of bubble generation is higher for greater wall superheat. 4d. Nonalternating Bubble Generation. Many previous studies assumed alternating bubble generations from the node and antinode points.24,8,12,14 However, from this study, it is found that the alternating character of bubble generations depends on the steadiness of the liquid pool and also the depth of the pool. Figure 12 shows the bubble generation from node

simulation is slightly higher than that derived by the Klimenko correlation. Figure 7 shows the heat flux distribution of the simulation performed for 4 s. Seven bubbles were released in this period of time. It is clear from this figure that the heat flux varies dynamically with bubble release. The time weighted average heat flux for this case is about 4200 W/m2, whereas the Berenson correlation estimates about 11000 W/m2 and Klimenko estimates about 2000 W/m2. It is to be noted that the Klimenko correlation was validated for LN2 and is more reliable than the Berenson correlation for simulating LN2 film boiling. Figure 8a, the velocity vector plot, shows that velocity is small in most parts of the liquid phase away from the bubble. The highest velocity exists in the vapor film. Vapor is moving toward the center of the bubble, i.e., toward the node. In the center line of the bubble, close to the necking point, the velocity is higher. It is also observed that high velocity gradients exist near the necking point as the liquid moves in to fill the gap. The minimum wall film thickness is observed close to the necking point. Thus, a maximum amount of vapor is produced in that zone. After the bubble was released, as shown in Figure 8b, a portion of vapor from the bubble stem returned to the film due to capillary forces. Vapor in the film travels toward the antinode. However, the vapor generated from the minimum thickness point is still moving toward the node. When the bubble is released, wake formation in the liquid below the bottom interface of the bubble pushes the bottom interface up, forming an inverted cup formation. As the bubble rises, it tries to adjust the spherical shape due to the surface tension. However, the inverted cup formation becomes ellipsoidal (or more like a cap shape) as it grows upward. Meanwhile, the disturbance travels toward the antinode as shown in Figure 8b, and subsequently, the bubble forms at the antinode. From Figure 8, it is seen that vapor velocity is higher than the liquid velocity. This can also be correlated with the temperature distribution profile. The temperature distribution as shown in Figure 9 indicates that the liquid is mostly uniform at its boiling point, whereas the vapor phase is superheated. From a comparison of both figures, it can be reasonably concluded that the velocity vector is directly related with the temperature.

Figure 9. Temperature distribution during the film boiling for LN2 at 73 K wall superheat of solution time = 0.25 s. G


Article


Figure 10. Streamlines showing the vortex formations for LNG at a wall superheat of ΔT = 43 K.

which assume the alternating nature of bubble release locations may cause erroneous estimation of vapor generation.

5. CONCLUSIONS Inaccurate estimations of LNG vapor generations due to accidental spill may cause underestimation of consequence severity during the risk assessment of the facility. To provide a more realistic estimation of the vapor dispersion source term, e.g., the vapor generation rate from a LNG pool, CFD simulations of cryogenic film boiling using volume of fluid (VOF) method in ANSYS Fluent have been studied. The simulations provide insights into the physical processes of vapor formation that are useful for estimating the LNG vapor generation. The key conclusions of this study are the following: (A) The simulated wall heat flux of LN2 film boiling is slightly greater than that estimated by using the Klimenko correlation and much lower than the Berenson correlation. It is important to note that the Klimenko correlation was validated for liquid nitrogen while the Berenson correlation was validated using noncryogenic liquid. (B) The simulated wall heat fluxes for LNG, as pure methane, were found to be significantly higher than the estimates of both the Klimenko and the Berenson correlations. As a result, the vapor generation estimation during the risk assessment would be underpredicted if such correlations were used. A conservative approach to determine the consequence severity of an accidental spill would be to use the estimated heat flux via CFD simulations. (C) The frequency of bubble generation is dependent on the degree of wall superheats. Increase of wall superheats increases the bubble generation rate; i.e., bubbles form faster. It also enhances the instabilities in the vapor film. Thus, during the early stage of an accidental spill of LNG, the vapor generation rate will be higher in comparison to the later stages of a spill. (D) It is observed that the bubble release locations depend on the dynamics of the vapor film movement, the motion of the liquid, and the movement of the liquid surface when the pool depth is not significant. As a result, the alternating nature of bubble release locations, from the node point to the antinode point in the consecutive bubble cycles, will be an unrealistic assumption. Thus, use of empirical expressions, which assume the alternating nature of bubble generation from the node and antinode points, might be unsuitable for cryogenic vaporization source term estimation.

Figure 11. Surface heat flux for the film boiling of LNG at a wall superheat ΔT = 43 K.

and antinode points for different wall superheat during the film boiling of liquid nitrogen. In Figure 12, if the bubble is released from the right side of the contour, it is said to be released from the node point whereas if it is released from the left side of the contour it is said to be released from antinode points. As depicted in Figure 12, for a wall superheat of 32 K, the first five bubbles were released from the node point, and the sixth bubble was released from the antinode point. Similarly, for ΔT = 73 K, the first and second bubbles were released from node points, and subsequently, three bubbles were released from the antinode points. For ΔT = 103 K, the first bubble was released from the node point; the second and third bubbles were released from the antinode point, and the subsequent three bubbles were released from node points. The dependency of the bubble release location depends on the motion of the vapor which is further influenced by force created at the vapor−liquid interface by the liquid velocity. When the depth of the pool is significantly small, the top surface movement also influences the motion of the liquid and, therefore, the release location of the bubble generation. During the accidental spill of cryogenic liquid, e.g., LNG, the film boiling will occur mainly in the spreading pool front where the depth of the pool is significantly low. Therefore, bubbles generating from the spreading front may not follow the switching nature of bubble release locations from node point to antinode point. Thus, empirical expressions H


Article


Figure 12. Effect of wall superheat on bubble generation frequency during LN2 film boiling.

■

AUTHOR INFORMATION

ṁ qp = mass transfer rate from the qth phase to the pth phase, kg/s δ = thickness of the initial vapor film, m n = index; number of variables; iteration number; total number of phases V = volume of a cell, m3 F⃗ = body force, N/m3 Fvol = surface tension equivalent body force, N/m3 E = total enthalpy, J/m3

Corresponding Author

*E-mail: [email protected]. Tel: +1 (979) 862-3985. Notes

The authors declare no competing financial interest.

■

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).

Dimensionless Numbers

■

f1, f 2, β = dimensionless number Nu = Nusselt number Ar = Archimedes number Pr = Prandtl number

NOMENCLATURE σ = surface tension, N/m ρ = density, kg/m3 g = acceleration of gravity, 9.81 m/s2 λ = wavelength, m hfg = heat of vaporization, J/kg k = thermal conductivity, W/m K μ = dynamic viscosity, Pa.S p = pressure, Pa Cp = heat capacity, J/kg K T = temperature, K ΔT = wall superheat, K t = time, s Δt = time step, s Lc = characteristic length, m α = volume fraction; vapor volume fraction h̅ = average local heat transfer coefficient, W/m2 K v ⃗ = fluid velocity vector, m/s Sαq = mass source term for qth phase, kg/m3

Subscripts

l = liquid phase v = vapor phase c = critical D; d2 = most dangerous; two-dimensional most dangerous D3 = three-dimensional most dangerous sat = saturated condition p, q = pth and qth phase wall = at wall condition x, y = in x and y direction b = at boiling point f = dimension index, 1 or 2 Superscripts

n or n+1 = value at nth or (n + 1)th iteration I


Article


■

Interface Tracking Method. Int. J. Heat Mass Transfer 2008, 51, 1646− 1657. (25) Agarwal, D. K.; Welch, S. W. J.; Biswas, G.; Durst, F. Planar Simulation of Bubble Growth in Film Boiling in Near-Critical Water Using a Variant of the VOF Method. J. Heat Transfer 2004, 126, 329. (26) Tomar, G.; Biswas, G.; Sharma, A.; Agrawal, A. Numerical Simulation of Bubble Growth in Film Boiling Using a Coupled LevelSet and Volume-of-Fluid Method. Phys. Fluids 2005, 17, 112103. (27) Tomar, G.; Biswas, G.; Sharma, A.; Welch, S. W. J. Multimode Analysis of Bubble Growth in Saturated Film Boiling. Phys. Fluids 2008, 20, 092101. (28) Hens, A.; Biswas, G.; De, S. Analysis of Interfacial Instability and Multimode Bubble Formation in Saturated Boiling Using Coupled Level Set and Volume-of-Fluid Approach. Phys. Fluids 2014, 26, 012105. (29) Welch, S. W. J.; Biswas, G. Direct Simulation of Film Boiling Including Electrohydrodynamic Forces. Phys. Fluids 2007, 19, 012106. (30) Tomar, G.; Gerlach, D.; Biswas, G.; Alleborn, N.; Sharma, A.; Durst, F.; Welch, S. W. J.; Delgado, A. Two-Phase Electrohydrodynamic Simulations Using a Volume-of-Fluid Approach. J. Comput. Phys. 2007, 227, 1267−1285. (31) Liu, Y.; Olewski, T.; Véchot, L. N. Modeling of a Cryogenic Liquid Pool Boiling by CFD Simulation. J. Loss Prev. Process Ind. 2015, 35, 125−134. (32) Hirt, C.; Nichols, B. Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. J. Comput. Phys. 1981, 39, 201−225. (33) Hardt, S.; Wondra, F. Evaporation Model for Interfacial Flows Based on a Continuum-Field Representation of the Source Terms. J. Comput. Phys. 2008, 227, 5871−5895. (34) Kunkelmann, C.; Stephan, P. CFD Simulation of Boiling Flows Using the Volume-of-Fluid Method within OpenFOAM. Numer. Heat Transfer, Part A 2009, 56, 631−646. (35) Kunugi, T. Brief Review of Latest Direct Numerical Simulation on Pool and Film Boiling. Nucl. Eng. Technol. 2012, 44, 847−854. (36) Ansys Inc. ANSYS Fluent Theory Guide; 2013, No. November; 591. (37) Nichols, B. D.; Hirt, C. W.; Hotchkiss, R. SOLA-VOF: A Solution Algorithm of Transient Flow Flow with Multiple Free Boundaries; Los Alamos Scientific Lab: Los Alamos, NM, 1980. (38) Rudman, M. VOLUME-TRACKING METHODS FOR INTERFACIAL FLOW CALCULATIONS. Int. J. Numer. Methods Fluids 1997, 24, 671−691. (39) Ubbink, O.; Issa, R. I. A Method for Capturing Sharp Fluid Interfaces on Arbitrary Meshes. J. Comput. Phys. 1999, 153, 26−50. (40) Noh, W. F.; Woodward, P. SLIC (Simple Line Interface Calculation). In The Fifth Internatioal Conference on Numerical Methods in Fluid Dynamics; Springer and Twente University: Berlin and Enschede, 1976. (41) Brackbill, J. U.; Kothe, D. B.; Zemach, C. A Continuum Method for Modeling Surface Tension. J. Comput. Phys. 1992, 100, 335−354. (42) Ansys Inc. Horizontal Film Boiling (Tutorial); 2011; pp 1−20. (43) Gibou, F.; Chen, L.; Nguyen, D.; Banerjee, S. A Level Set Based Sharp Interface Method for the Multiphase Incompressible Navier– Stokes Equations with Phase Change. J. Comput. Phys. 2007, 222, 536−555. (44) Rowley, R. L.; Wilding, W. V.; Oscarson, J. L.; Zundel, J. L.; Marshall, N. A.; Daubert, T. L.; Danner, R. P. DIPPR Data Compilation of Pure Compund Properties; Design Institute for Physical Property Data/AIChE, 2003. (45) Barron, R. F. Cryogenic Heat Transfer; Taylor & Francis: Philadelphia, PA, 1999. (46) Berenson, P. J. Experiments on Pool-Boiling Heat Transfer. Int. J. Heat Mass Transfer 1962, 5, 985−999.

REFERENCES

(1) United States Government Accountability Office. Report GAO07-316-Maritime Security: Public Safety Consequences of a Terrorist Attack on a Tanker Carrying Liquefied Natural Gas Need Clarification Highlights; 2007. (2) Havens, J.; Spicer, T. United States Regulations for Siting LNG Terminals: Problems and Potential. J. Hazard. Mater. 2007, 140, 439− 443. (3) NFPA 59A. Standard for the Production, Storage, and Handling of Liquefied Natural Gas (LNG); National Fire Protection Association, 2013. (4) Véchot, L.; Olewski, T.; Osorio, C.; Basha, O.; Liu, Y.; Mannan, S. Laboratory Scale Analysis of the Influence of Different Heat Transfer Mechanisms on Liquid Nitrogen Vaporization Rate. J. Loss Prev. Process Ind. 2013, 26, 398−409. (5) Webber, D. M.; Gant, S. E.; Ivings, M. J.; Jagger, S. F. LNG Source Term Models for Hazard Analysis: A Review of the State-of-the-Art and an Approach to Model Assessment; Research Report: RR789; Buxton, Derbyshire, UK, 2010. (6) Zuber, N. Hydrodynamic Aspects of Boiling Heat Transfer; University of California: Los Angeles, 1959. (7) Carey, V. P. Liquid-Vapor Phase Change Phenomena: An Introduction to the Thermophysics of Vaporization and Condensation Process in Heat Transfer Equipment; Hemisphere Publishing Corp: Washington, DC, 1992. (8) Berenson, P. J. Film-Boiling Heat Transfer From a Horizontal Surface. J. Heat Transfer 1961, 83, 351−356. (9) Sernas, V.; Lienhard, J. H.; Dhir, V. K. The Taylor Wave Configuration During Boiling from a Falt Plate. Int. J. Heat Mass Transfer 1973, 16, 1820−1821. (10) Hoster, E. R.; Westwater, J. W. Film Boiling on a Horizontal Plate. Am. Rocket Soc. J. 1962, 32, 553−558. (11) Klimenko, V. V. Film Boiling on a Horizontal Plate - New Correlation. Int. J. Heat Mass Transfer 1981, 24, 69−79. (12) Klimenko, V. V.; Shelepen, A. G. Film Boiling on Horizontal Plate - a Supplementary Communication. Int. J. Heat Mass Transfer 1982, 25, 1611−1613. (13) Son, G.; Dhir, V. K. Numerical Simulation of Saturated Film Boiling on a Horizontal Surface. J. Heat Transfer 1997, 119, 525−533. (14) Dhir, V. K. Numerical Simulations of Pool-Boiling Heat Transfer. AIChE J. 2001, 47, 813−834. (15) Son, G.; Dhir, V. K. Numerical Simulation of Film Boiling Near Critical Pressures With a Level Set Method. J. Heat Transfer 1998, 120, 183−192. (16) Panzarella, C. H.; Davis, S. H.; Bankoff, S. G. Nonlinear Dynamics in Horizontal Film Boiling. J. Fluid Mech. 2000, 402, 163− 194. (17) Banerjee, D.; Dhir, V. K. Study of Subcooled Film Boiling on a Horizontal Disc: Part IAnalysis. J. Heat Transfer 2001, 123, 271− 284. (18) Juric, D.; Tryggvason, G. Computations of Boiling Flows. Int. J. Multiphase Flow 1998, 24, 387−410. (19) Esmaeeli, A.; Tryggvason, G. Computations of Film Boiling. Part I: Numerical Method. Int. J. Heat Mass Transfer 2004, 47, 5451− 5461. (20) Esmaeeli, A.; Tryggvason, G. Computations of Film Boiling. Part II: Multi-Mode Film Boiling. Int. J. Heat Mass Transfer 2004, 47, 5463−5476. (21) Welch, S. W. J.; Wilson, J. A Volume of Fluid Based Method for Fluid Flows with Phase Change. J. Comput. Phys. 2000, 160, 662−682. (22) Youngs, D. L. Time Dependent Multi-Material Flow with Large Fluid Distortion. In Numerical Methods in Fluid Dynamics; Morton, K. W., Baines, M. J., Eds.; Academic Press, 1982; pp 273−285. (23) Welch, S. W. J.; Rachidi, T. Numerical Computation of Film Boiling Including Conjugate Heat Transfer. Numer. Heat Transfer, Part B 2002, 42, 35−53. (24) Yuan, M. H.; Yang, Y. H.; Li, T. S.; Hu, Z. H. Numerical Simulation of Film Boiling on a Sphere with a Volume of Fluid J


Application of Computational Fluid Dynamics in ... - ACS Publications

Recommend Documents