- Email: [email protected]
Simulation of two-phase flow and heat transfer phenomena in a boiling water reactor using the lattice Boltzmann method
ARTICLE IN PRESS
Physica A 362 (2006) 98–104 www.elsevier.com/locate/physa
Simulation of two-phase flow and heat transfer phenomena in a boiling water reactor using the lattice Boltzmann method Adrian Tentnera,, Hudong Chenb, Raoyang Zhangb a
Argonne National Laboratory, 9700 S. Cass Ave., Darien, IL 60439, USA Exa Corporation, 3 Burlington Woods Drive, Burlington, MA 01803, USA
b
Available online 17 October 2005
Abstract The paper describes a new computational tool based on lattice Boltzmann methods for the simulation of two-phase flow and heat transfer phenomena in boiling water reactor fuel bundles. r 2005 Elsevier B.V. All rights reserved. Keywords: Two-phase flow; Lattice Boltzmann method; Boiling water reactor
1. Introduction The work described has developed a new computational tool based on lattice Boltzmann methods (LBM) for the simulation of two-phase flow and heat transfer phenomena in boiling water reactor (BWR) fuel bundles. The designers and operators of BWRs are highly interested in improving their understanding of the detailed two-phase flow phenomena that occur inside a BWR fuel bundle. These phenomena include coolant phase changes and multiple flow regime transitions, which directly influence the coolant interaction with the fuel pins, spacers, and other elements of the fuel assembly [1]. These complex interactions between the coolant and the fuel assembly ultimately determine the reactor performance. 2. Background Traditionally, the best tools for the analysis of two-phase flow phenomena inside the BWR fuel assembly have been the sub-channel codes. However, the resolution of these codes is still too coarse for analyzing the detailed intra-assembly flow patterns, such as flow around a spacer element. Recent progress in computational fluid dynamics (CFD), coupled with the rapidly increasing computational power of massively parallel computers, shows promising potential for the fine-mesh, detailed simulation of fuel assembly two-phase flow phenomena. However, the phenomenological models available in the conventional CFD codes solving the Navier–Stokes equations are not as advanced as those currently being used in the sub-channel codes used in the nuclear industry. In particular, there are no models that are able to reliably predict the nature of the flow Corresponding author.
E-mail addresses: [email protected] (A. Tentner), [email protected] (H. Chen), [email protected] (R. Zhang). 0378-4371/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2005.09.028
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
99
regimes, and use the appropriate sub-models for those flow regimes are currently available. LBM-based codes which use a mesoscopic representation of the fluid have been used successfully in modeling the single-phase fluid flow in very complex geometries. They have the potential, which has not yet been fully exploited, to simulate complex two-phase flow phenomena characteristic for a BWR fuel assembly, including flow regime transitions, at a more fundamental level than possible with conventional CFD approaches. 3. Multiphase lattice Boltzmann algorithm with thermodynamics The lattice Boltzmann method has been shown to be ideally suited for the simulation of flows involving multiple immiscible phases [2–4]. However, all previously existing LBM codes are only applicable for isothermal flow situations, since no technique for dynamic coupling of temperature with the multiphase lattice Boltzmann algorithm existed in the published literature [4]. The reason for this limitation is that the lattice Boltzmann models involving temperature dynamics based on a self-consistently generated energy are known to be numerically unstable [5] and the computational cost for overcoming this problem is very high in practical applications [6]. The isothermal flow simplifying assumptions make it impossible to track the energy transport and conservation that are required to simulate the temperature dynamics needed to describe the water–vapor mixture flow typical for BWR applications. Another obstacle to using an LBM based code for BWR twophase flow simulations is the need to use non-local interaction terms based on the water equation of state (EOS), which introduce additional computational and numerical stability problems. In this work we have successfully addressed the LBM limitations described above. We developed a novel hybrid LBM approach that allows the simulation of two-phase fluid thermodynamics, including mass exchange between phases. This model describes the fluid thermodynamics by modeling the multiphase mass and momentum evolution via the lattice Boltzmann method, while the temperature dynamics is computed by solving a scalar transport energy equation [7]. The model also includes terms that describe the non-local interactions based on the actual water/ steam EOS, allowing the simulation of BWR two-phase phenomena. Modeling of turbulence and interfacial forces was not included in these calculations. Both will be explored in future work. 3.1. Implementation of lattice Boltzmann algorithm with non-ideal gas equation of state We have implemented a lattice Boltzmann algorithm that simulates non-ideal gas fluid flow. This approach has demonstrated its ability for capturing the intrinsic phase separation physics in boiling processes involving a mixture of water and vapor. The basic lattice Boltzmann equation is given as f i ðx þ ci ; t þ 1Þ f i ðx; tÞ ¼ C i ðx; tÞ þ df i .
(1)
The collision term C i ðx; tÞ is given by the so called BGK model that determines the viscosity in terms of the relaxation time: 1 C ¼ ðf f eq Þ. (2) t It has been shown that lattice Boltzmann models of this kind gives rise to the Navier–Stokes fluid flow [8–10]: qt r þ r ðr~ uÞ ¼ 0,
(3)
rðqt~ uþ~ u r~ uÞ ¼ rp þ r ðrn r~ uÞ,
(4)
where r and ~ u are, respectively, the fluid density and velocity. Without the additional term df i the pressure p in (4) obeys an ideal gas-like equation of state, p ¼ RT 0 r, where T 0 is a prescribed temperature value and R is the universal gas constant. The resulting kinematic viscosity is determined by the relaxation parameter t: n ¼ ðt 1=2ÞRT 0 .
(5)
The local fluid viscosity can be changed by modifying the relaxation time t. However, viscosity effects were not studied in these initial calculations, which were conducted using a constant relaxation time t ¼ 1.
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
100
One of the key developments of our work was to use the additional term df i in Eq. (1) in order to produce desired non-ideal gas properties. Generally speaking, this additional term represents a body force. The first key observation is that a non-ideal gas equation of state is produced if this body force term is self-consistently generated from the gradient of the local density and temperature values. Specifically, the momentum due to the additional term df i in (1) is explicitly constructed such that X ci df i ¼ rpðr; TÞ þ rðRrT 0 Þ. (6) i
Here pðr; TÞ can be in principle any arbitrary function of local density and temperature, for example the one corresponding to the classical Van der Waals equation of state or a more realistic one. The second term in (6) is used to cancel the existing ideal gas pressure in the original lattice model. As a consequence, an algorithm describing the thermodynamics of a fluid with a non-ideal gas equation of state is obtained when the evolution of the temperature distribution or an energy equation is provided. Although the surface tension treatment was not included in these calculations, it is noted that the new algorithm is an extension of the two-phase LBM model of Shan and Chen [2,3] and implicitly contains a surface tension mechanism [11] which depends on the formulation of Eq. (6). Matching of the actual water–vapor surface tension is possible through the addition of local phase interface terms in (6) and will be explored in future work. Fig. 1 shows a set of graphs describing fluid density distributions at different times in a square domain with heated walls, resembling the cross-section of a BWR fuel assembly. The results shown in Fig. 1 were obtained for a water/steam system with a temperature T ¼ 0:9T c , where T c is the critical temperature, and initial density r chosen to be within the phase co-existence domain. Spontaneous phase separation between water and vapor phases is predicted by the new lattice Boltzmann code. Vapor bubbles (red color) and water (blue color) emerge out of the initial two-phase mixture as a result of infinitesimally small random density fluctuations. In time, the small vapor bubbles coalesce into ever larger size bubbles which directly describes the expected behavior of realistic fluid undergoing phase transition. To demonstrate the physics of this process, we used the Van der Waals non-ideal gas equation of state for the fluid equation of state (EOS). In this case, pressure is related to the local specific volume V (i.e. the inverse number density V ¼ 1/r) and temperature T according to a p þ 2 ðV bÞ ¼ RT, (7) V where a and b are constant parameters which depend on the fluid being analyzed. The Van der Waals EOS has the advantage of having a simple analytical form and its thermodynamic properties are well understood and defined. This makes it suitable for demonstrating the fundamental capability of the current approach for multiphase flow simulations. Fig. 2 plots the pressure p isotherms given by the Van der Waals equation (7) versus the specific volume V at three representative temperatures T. The key difference between a fluid with a non-ideal gas EOS (such as Van der Waals) and an ideal gas is that at temperatures lower than T c the pressure of such a fluid is no longer a monotonically decreasing function of the specific volume so that a given pressure corresponds to more than one density value. Consequently, a substance with an initial density in the two-phase range will undergo a spontaneous phase separation into two immiscible phases. As shown in Fig. 2, the resulting density values at a given temperature are located at the intersections of a given pressure isotherm (solid line with symbols) and the
Fig. 1. Calculated spontaneous water/vapor separation process. Three consecutive times are shown.
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
101
Fig. 2. Van der Waals isotherms for a water/steam system.
phase coexistence curve (the solid blue line). This fundamental physical mechanism is robust across equations of state and when implemented in the algorithm developed during our work realistically describes formation of two phases typical of water–vapor mixture. Also in Fig. 2, we insert an isotherm (the purple line with symbols) of the water–vapor equation of state obtained from the experimental water–steam tables. The results of calculations performed using the EOS based on the water–steam tables are presented in Section 4. 3.2. Implementation of a scalar solver for the temperature dynamics As discussed above, the simulation of two-phase flow including phase changes requires specification of the temperature dynamics. For this purpose we have developed and incorporated in the new LBM model a finite difference scheme for the solution of the energy equation, which describes the scalar transport of energy in two-phase flow. The specific form of the energy equation uses the fluid enthalpy h: rðqt h þ ~ u rhÞ ¼ ðqt p þ ~ u rpÞ þ r ðrkrTÞ.
(8)
One motivation for this choice is that the pressure distribution (and its functional dependencies) is continuous leading to smooth changes of the first term on the right-hand side of (8), unlike the corresponding term in the internal energy equation which is highly discontinuous across interfaces due to density jumps. Furthermore, the enthalpy equation in the above differential form more accurately describes the physics of multiphase flows than the internal energy equation. Once the enthalpy h is computed, the temperature T is obtained using a function T ¼ Tðh; pÞ based on the water/steam tables which covers the entire phase spectrum, including liquid, two-phase mixture, and vapor. The use of this function naturally and accurately accounts for the latent heat effects present during phase changes. In a typical BWR channel during steady-state operation the pressure gradients are small and the pressure terms on the right-hand side of (8) are negligible when compared with the conductivity term. Therefore, the effect of the pressure terms on the enthalpy behavior was neglected in the initial calculations reported in Section 4. Future work will require additional consideration with respect to the treatment of the pressure time and space derivatives in the enthalpy equation during transients or in systems with large pressure gradients. The liquid and vapor conductivities used in the enthalpy calculations were obtained from the water/steam tables. The two-phase mixture conductivity was calculated through an interpolation between the liquid and vapor conductivities based on the local void fraction.
ARTICLE IN PRESS 102
A. Tentner et al. / Physica A 362 (2006) 98–104
In solving enthalpy equation (8) we separated the changes in enthalpy into two categories: (a) enthalpy changes due to fluid convection, and (b) enthalpy changes due to conduction. The enthalpy solution was therefore implemented in two steps. In the first step the enthalpy changes due to fluid convection were calculated, following the mass and enthalpy transport along the lattice directions. This approach was essential in ensuring the exact energy conservation. In the second step we calculated the enthalpy changes due to conduction along the Cartesian x–y coordinate directions, using a temperature (T) based formulation of the last term in (8). We completed the basic software protocol of the new hybrid approach by coupling this scalar energy equation to the above multiphase lattice Boltzmann algorithm through the equation of state. 4. Results of BWR channel analysis To demonstrate our novel approach for LBM simulation of two-phase flow we computed water–vapor phase transition in a vertical heated channel with parameters similar to a BWR fuel assembly sub-channel. In these initial calculations we have studied a 2-D geometry using 32 256 grid points which allows us to obtain prototypical results with a high computational efficiency. To reduce the computational effort the length/width ratio of the channel is smaller than that of the actual BWR channel, and the wall heat flux has been increased accordingly to allow the energy transferred to fluid to match the actual BWR channel value. Sub-cooled single-phase liquid enters the bottom of the pipe (T ¼ 262:7 C, r ¼ 790:64 Kg=m3 Þ with a uniform fixed velocity V ¼ 1:5 m=s, is heated by the heat flux received from the walls as it moves upwards, and begins to vaporize. The vapor leaves the vertical heated channel at the upper end, where a fixed pressure boundary condition is used (P ¼ 68:25e5 Pa). The equation of state used in the simulation is based on a routine that reproduces the water/steam tables for the range of interest. At the start of the simulation the channel is filled with water at T ¼ 262:7 C. The coolant density distributions in the channel after 10,000 and 11,000 computational time steps are shown in Figs. 3A and B, respectively. The number of steps required to traverse the channel is 12,800. The sub-cooled fluid enters at the bottom of the channel, is heated by the channel walls while flowing upwards and exits at the top of the channel. We can clearly see the formation of the low density vapor-phase region at the outlet of the
Fig. 3. Density (Kg/m3) distribution in the BWR channel at t ¼ 10 000 cycles (A) and t ¼ 11 000 cycles (B).
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
103
Fig. 4. Radial vapor volume fraction distribution at t ¼ 10 000 cycles.
channel in both figures. Two lower density regions can be seen below the continuous vapor region in Fig. 3A, indicating the formation of vapor bubbles near the walls. In Fig. 3B two irregular lower density regions, resembling streams of coalescing bubbles can be observed near the walls, below the continuous vapor region. The transition from these regions to the continuous vapor region is similar to the transition from the continuous liquid to continuous vapor flow regimes observed in BWR channels. Both figures show the formation of liquid films on the heated channel walls, an important feature of the two-phase flow in BWR fuel bundles. It is noted however that the phase separation in the bubbly flow regime is not as clearly defined as the separation illustrated in Fig. 1, where the Van der Waals EOS was used. The lower definition of the phase interfaces is due in part to the coarser computational mesh used in the BWR calculations. This effect may also be due to the use of the water/steam table EOS which does not display the two spinoidal points and the positive slope region of the pressure isotherms typical for the Van der Waals EOS. The effects of the mesh resolution and the EOS shape on phase separation will be studied in future work. As we approach steady state, the high density (water) phase and low density (vapor) phase show sharper interfaces and a configuration typical of a BWR channel. Fig. 4 illustrates the radial volume fraction distribution at selected axial locations at 10 000 cycles. The flattening profiles at the exit are indicative of the transition from fully developed bubbly or slug flow to the annular flow regime. 5. Conclusion A new lattice Boltzmann model was developed that includes thermodynamics and allows the modeling of non-ideal gas fluids. Preliminary calculations have demonstrated the capability of the model to simulate twophase flow phenomena typical for a BWR fuel assembly. This new model opens the way for the simulation abinitio of many two-phase flow phenomena of interest to BWR fuel assembly design, including flow regime transitions. References [1] R.T. Lahey, F.J. Moody, The Thermal-Hydraulics of a Boiling Water Nuclear Reactor, American Nuclear Society, 1993. [2] X. Shan, H. Chen, Phys. Rev. E 47 (1993) 1815. [3] X. Shan, H. Chen, Phys. Rev. E 49 (1994) 2941.
ARTICLE IN PRESS 104 [4] [5] [6] [7] [8] [9] [10] [11]
A. Tentner et al. / Physica A 362 (2006) 98–104 S. Chen, G. Doolen, Ann. Rev. Fluid Mech. 30 (1998) 329. H. Chen, C. Teixeira, K. Molvig, Int. J. Mod. Phys. C 8 (4) (1997) 675. C. Teixeira, H. Chen, D. Freed, Comput. Phys. Commun. 129 (2000) 207. R. Zhang, H. Chen, Phys. Rev. E 67 (2003) 066711. S. Chen, H. Chen, W. Matthaeus, D. Martinez, Phys. Rev. Lett. 67 (1991) 3776. H. Chen, S. Chen, W. Matthaeus, Phys. Rev. E 45 (1992) 5339. Y. Qian, D. d’Humieres, P. Lallemand, Europhys. Lett. 17 (1992) 479. X. He, S. Chen, R. Zhang, J. Comput. Phys. 152 (1999) 642–663.
Physica A 362 (2006) 98–104 www.elsevier.com/locate/physa
Simulation of two-phase flow and heat transfer phenomena in a boiling water reactor using the lattice Boltzmann method Adrian Tentnera,, Hudong Chenb, Raoyang Zhangb a
Argonne National Laboratory, 9700 S. Cass Ave., Darien, IL 60439, USA Exa Corporation, 3 Burlington Woods Drive, Burlington, MA 01803, USA
b
Available online 17 October 2005
Abstract The paper describes a new computational tool based on lattice Boltzmann methods for the simulation of two-phase flow and heat transfer phenomena in boiling water reactor fuel bundles. r 2005 Elsevier B.V. All rights reserved. Keywords: Two-phase flow; Lattice Boltzmann method; Boiling water reactor
1. Introduction The work described has developed a new computational tool based on lattice Boltzmann methods (LBM) for the simulation of two-phase flow and heat transfer phenomena in boiling water reactor (BWR) fuel bundles. The designers and operators of BWRs are highly interested in improving their understanding of the detailed two-phase flow phenomena that occur inside a BWR fuel bundle. These phenomena include coolant phase changes and multiple flow regime transitions, which directly influence the coolant interaction with the fuel pins, spacers, and other elements of the fuel assembly [1]. These complex interactions between the coolant and the fuel assembly ultimately determine the reactor performance. 2. Background Traditionally, the best tools for the analysis of two-phase flow phenomena inside the BWR fuel assembly have been the sub-channel codes. However, the resolution of these codes is still too coarse for analyzing the detailed intra-assembly flow patterns, such as flow around a spacer element. Recent progress in computational fluid dynamics (CFD), coupled with the rapidly increasing computational power of massively parallel computers, shows promising potential for the fine-mesh, detailed simulation of fuel assembly two-phase flow phenomena. However, the phenomenological models available in the conventional CFD codes solving the Navier–Stokes equations are not as advanced as those currently being used in the sub-channel codes used in the nuclear industry. In particular, there are no models that are able to reliably predict the nature of the flow Corresponding author.
E-mail addresses: [email protected] (A. Tentner), [email protected] (H. Chen), [email protected] (R. Zhang). 0378-4371/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2005.09.028
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
99
regimes, and use the appropriate sub-models for those flow regimes are currently available. LBM-based codes which use a mesoscopic representation of the fluid have been used successfully in modeling the single-phase fluid flow in very complex geometries. They have the potential, which has not yet been fully exploited, to simulate complex two-phase flow phenomena characteristic for a BWR fuel assembly, including flow regime transitions, at a more fundamental level than possible with conventional CFD approaches. 3. Multiphase lattice Boltzmann algorithm with thermodynamics The lattice Boltzmann method has been shown to be ideally suited for the simulation of flows involving multiple immiscible phases [2–4]. However, all previously existing LBM codes are only applicable for isothermal flow situations, since no technique for dynamic coupling of temperature with the multiphase lattice Boltzmann algorithm existed in the published literature [4]. The reason for this limitation is that the lattice Boltzmann models involving temperature dynamics based on a self-consistently generated energy are known to be numerically unstable [5] and the computational cost for overcoming this problem is very high in practical applications [6]. The isothermal flow simplifying assumptions make it impossible to track the energy transport and conservation that are required to simulate the temperature dynamics needed to describe the water–vapor mixture flow typical for BWR applications. Another obstacle to using an LBM based code for BWR twophase flow simulations is the need to use non-local interaction terms based on the water equation of state (EOS), which introduce additional computational and numerical stability problems. In this work we have successfully addressed the LBM limitations described above. We developed a novel hybrid LBM approach that allows the simulation of two-phase fluid thermodynamics, including mass exchange between phases. This model describes the fluid thermodynamics by modeling the multiphase mass and momentum evolution via the lattice Boltzmann method, while the temperature dynamics is computed by solving a scalar transport energy equation [7]. The model also includes terms that describe the non-local interactions based on the actual water/ steam EOS, allowing the simulation of BWR two-phase phenomena. Modeling of turbulence and interfacial forces was not included in these calculations. Both will be explored in future work. 3.1. Implementation of lattice Boltzmann algorithm with non-ideal gas equation of state We have implemented a lattice Boltzmann algorithm that simulates non-ideal gas fluid flow. This approach has demonstrated its ability for capturing the intrinsic phase separation physics in boiling processes involving a mixture of water and vapor. The basic lattice Boltzmann equation is given as f i ðx þ ci ; t þ 1Þ f i ðx; tÞ ¼ C i ðx; tÞ þ df i .
(1)
The collision term C i ðx; tÞ is given by the so called BGK model that determines the viscosity in terms of the relaxation time: 1 C ¼ ðf f eq Þ. (2) t It has been shown that lattice Boltzmann models of this kind gives rise to the Navier–Stokes fluid flow [8–10]: qt r þ r ðr~ uÞ ¼ 0,
(3)
rðqt~ uþ~ u r~ uÞ ¼ rp þ r ðrn r~ uÞ,
(4)
where r and ~ u are, respectively, the fluid density and velocity. Without the additional term df i the pressure p in (4) obeys an ideal gas-like equation of state, p ¼ RT 0 r, where T 0 is a prescribed temperature value and R is the universal gas constant. The resulting kinematic viscosity is determined by the relaxation parameter t: n ¼ ðt 1=2ÞRT 0 .
(5)
The local fluid viscosity can be changed by modifying the relaxation time t. However, viscosity effects were not studied in these initial calculations, which were conducted using a constant relaxation time t ¼ 1.
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
100
One of the key developments of our work was to use the additional term df i in Eq. (1) in order to produce desired non-ideal gas properties. Generally speaking, this additional term represents a body force. The first key observation is that a non-ideal gas equation of state is produced if this body force term is self-consistently generated from the gradient of the local density and temperature values. Specifically, the momentum due to the additional term df i in (1) is explicitly constructed such that X ci df i ¼ rpðr; TÞ þ rðRrT 0 Þ. (6) i
Here pðr; TÞ can be in principle any arbitrary function of local density and temperature, for example the one corresponding to the classical Van der Waals equation of state or a more realistic one. The second term in (6) is used to cancel the existing ideal gas pressure in the original lattice model. As a consequence, an algorithm describing the thermodynamics of a fluid with a non-ideal gas equation of state is obtained when the evolution of the temperature distribution or an energy equation is provided. Although the surface tension treatment was not included in these calculations, it is noted that the new algorithm is an extension of the two-phase LBM model of Shan and Chen [2,3] and implicitly contains a surface tension mechanism [11] which depends on the formulation of Eq. (6). Matching of the actual water–vapor surface tension is possible through the addition of local phase interface terms in (6) and will be explored in future work. Fig. 1 shows a set of graphs describing fluid density distributions at different times in a square domain with heated walls, resembling the cross-section of a BWR fuel assembly. The results shown in Fig. 1 were obtained for a water/steam system with a temperature T ¼ 0:9T c , where T c is the critical temperature, and initial density r chosen to be within the phase co-existence domain. Spontaneous phase separation between water and vapor phases is predicted by the new lattice Boltzmann code. Vapor bubbles (red color) and water (blue color) emerge out of the initial two-phase mixture as a result of infinitesimally small random density fluctuations. In time, the small vapor bubbles coalesce into ever larger size bubbles which directly describes the expected behavior of realistic fluid undergoing phase transition. To demonstrate the physics of this process, we used the Van der Waals non-ideal gas equation of state for the fluid equation of state (EOS). In this case, pressure is related to the local specific volume V (i.e. the inverse number density V ¼ 1/r) and temperature T according to a p þ 2 ðV bÞ ¼ RT, (7) V where a and b are constant parameters which depend on the fluid being analyzed. The Van der Waals EOS has the advantage of having a simple analytical form and its thermodynamic properties are well understood and defined. This makes it suitable for demonstrating the fundamental capability of the current approach for multiphase flow simulations. Fig. 2 plots the pressure p isotherms given by the Van der Waals equation (7) versus the specific volume V at three representative temperatures T. The key difference between a fluid with a non-ideal gas EOS (such as Van der Waals) and an ideal gas is that at temperatures lower than T c the pressure of such a fluid is no longer a monotonically decreasing function of the specific volume so that a given pressure corresponds to more than one density value. Consequently, a substance with an initial density in the two-phase range will undergo a spontaneous phase separation into two immiscible phases. As shown in Fig. 2, the resulting density values at a given temperature are located at the intersections of a given pressure isotherm (solid line with symbols) and the
Fig. 1. Calculated spontaneous water/vapor separation process. Three consecutive times are shown.
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
101
Fig. 2. Van der Waals isotherms for a water/steam system.
phase coexistence curve (the solid blue line). This fundamental physical mechanism is robust across equations of state and when implemented in the algorithm developed during our work realistically describes formation of two phases typical of water–vapor mixture. Also in Fig. 2, we insert an isotherm (the purple line with symbols) of the water–vapor equation of state obtained from the experimental water–steam tables. The results of calculations performed using the EOS based on the water–steam tables are presented in Section 4. 3.2. Implementation of a scalar solver for the temperature dynamics As discussed above, the simulation of two-phase flow including phase changes requires specification of the temperature dynamics. For this purpose we have developed and incorporated in the new LBM model a finite difference scheme for the solution of the energy equation, which describes the scalar transport of energy in two-phase flow. The specific form of the energy equation uses the fluid enthalpy h: rðqt h þ ~ u rhÞ ¼ ðqt p þ ~ u rpÞ þ r ðrkrTÞ.
(8)
One motivation for this choice is that the pressure distribution (and its functional dependencies) is continuous leading to smooth changes of the first term on the right-hand side of (8), unlike the corresponding term in the internal energy equation which is highly discontinuous across interfaces due to density jumps. Furthermore, the enthalpy equation in the above differential form more accurately describes the physics of multiphase flows than the internal energy equation. Once the enthalpy h is computed, the temperature T is obtained using a function T ¼ Tðh; pÞ based on the water/steam tables which covers the entire phase spectrum, including liquid, two-phase mixture, and vapor. The use of this function naturally and accurately accounts for the latent heat effects present during phase changes. In a typical BWR channel during steady-state operation the pressure gradients are small and the pressure terms on the right-hand side of (8) are negligible when compared with the conductivity term. Therefore, the effect of the pressure terms on the enthalpy behavior was neglected in the initial calculations reported in Section 4. Future work will require additional consideration with respect to the treatment of the pressure time and space derivatives in the enthalpy equation during transients or in systems with large pressure gradients. The liquid and vapor conductivities used in the enthalpy calculations were obtained from the water/steam tables. The two-phase mixture conductivity was calculated through an interpolation between the liquid and vapor conductivities based on the local void fraction.
ARTICLE IN PRESS 102
A. Tentner et al. / Physica A 362 (2006) 98–104
In solving enthalpy equation (8) we separated the changes in enthalpy into two categories: (a) enthalpy changes due to fluid convection, and (b) enthalpy changes due to conduction. The enthalpy solution was therefore implemented in two steps. In the first step the enthalpy changes due to fluid convection were calculated, following the mass and enthalpy transport along the lattice directions. This approach was essential in ensuring the exact energy conservation. In the second step we calculated the enthalpy changes due to conduction along the Cartesian x–y coordinate directions, using a temperature (T) based formulation of the last term in (8). We completed the basic software protocol of the new hybrid approach by coupling this scalar energy equation to the above multiphase lattice Boltzmann algorithm through the equation of state. 4. Results of BWR channel analysis To demonstrate our novel approach for LBM simulation of two-phase flow we computed water–vapor phase transition in a vertical heated channel with parameters similar to a BWR fuel assembly sub-channel. In these initial calculations we have studied a 2-D geometry using 32 256 grid points which allows us to obtain prototypical results with a high computational efficiency. To reduce the computational effort the length/width ratio of the channel is smaller than that of the actual BWR channel, and the wall heat flux has been increased accordingly to allow the energy transferred to fluid to match the actual BWR channel value. Sub-cooled single-phase liquid enters the bottom of the pipe (T ¼ 262:7 C, r ¼ 790:64 Kg=m3 Þ with a uniform fixed velocity V ¼ 1:5 m=s, is heated by the heat flux received from the walls as it moves upwards, and begins to vaporize. The vapor leaves the vertical heated channel at the upper end, where a fixed pressure boundary condition is used (P ¼ 68:25e5 Pa). The equation of state used in the simulation is based on a routine that reproduces the water/steam tables for the range of interest. At the start of the simulation the channel is filled with water at T ¼ 262:7 C. The coolant density distributions in the channel after 10,000 and 11,000 computational time steps are shown in Figs. 3A and B, respectively. The number of steps required to traverse the channel is 12,800. The sub-cooled fluid enters at the bottom of the channel, is heated by the channel walls while flowing upwards and exits at the top of the channel. We can clearly see the formation of the low density vapor-phase region at the outlet of the
Fig. 3. Density (Kg/m3) distribution in the BWR channel at t ¼ 10 000 cycles (A) and t ¼ 11 000 cycles (B).
ARTICLE IN PRESS A. Tentner et al. / Physica A 362 (2006) 98–104
103
Fig. 4. Radial vapor volume fraction distribution at t ¼ 10 000 cycles.
channel in both figures. Two lower density regions can be seen below the continuous vapor region in Fig. 3A, indicating the formation of vapor bubbles near the walls. In Fig. 3B two irregular lower density regions, resembling streams of coalescing bubbles can be observed near the walls, below the continuous vapor region. The transition from these regions to the continuous vapor region is similar to the transition from the continuous liquid to continuous vapor flow regimes observed in BWR channels. Both figures show the formation of liquid films on the heated channel walls, an important feature of the two-phase flow in BWR fuel bundles. It is noted however that the phase separation in the bubbly flow regime is not as clearly defined as the separation illustrated in Fig. 1, where the Van der Waals EOS was used. The lower definition of the phase interfaces is due in part to the coarser computational mesh used in the BWR calculations. This effect may also be due to the use of the water/steam table EOS which does not display the two spinoidal points and the positive slope region of the pressure isotherms typical for the Van der Waals EOS. The effects of the mesh resolution and the EOS shape on phase separation will be studied in future work. As we approach steady state, the high density (water) phase and low density (vapor) phase show sharper interfaces and a configuration typical of a BWR channel. Fig. 4 illustrates the radial volume fraction distribution at selected axial locations at 10 000 cycles. The flattening profiles at the exit are indicative of the transition from fully developed bubbly or slug flow to the annular flow regime. 5. Conclusion A new lattice Boltzmann model was developed that includes thermodynamics and allows the modeling of non-ideal gas fluids. Preliminary calculations have demonstrated the capability of the model to simulate twophase flow phenomena typical for a BWR fuel assembly. This new model opens the way for the simulation abinitio of many two-phase flow phenomena of interest to BWR fuel assembly design, including flow regime transitions. References [1] R.T. Lahey, F.J. Moody, The Thermal-Hydraulics of a Boiling Water Nuclear Reactor, American Nuclear Society, 1993. [2] X. Shan, H. Chen, Phys. Rev. E 47 (1993) 1815. [3] X. Shan, H. Chen, Phys. Rev. E 49 (1994) 2941.
ARTICLE IN PRESS 104 [4] [5] [6] [7] [8] [9] [10] [11]
A. Tentner et al. / Physica A 362 (2006) 98–104 S. Chen, G. Doolen, Ann. Rev. Fluid Mech. 30 (1998) 329. H. Chen, C. Teixeira, K. Molvig, Int. J. Mod. Phys. C 8 (4) (1997) 675. C. Teixeira, H. Chen, D. Freed, Comput. Phys. Commun. 129 (2000) 207. R. Zhang, H. Chen, Phys. Rev. E 67 (2003) 066711. S. Chen, H. Chen, W. Matthaeus, D. Martinez, Phys. Rev. Lett. 67 (1991) 3776. H. Chen, S. Chen, W. Matthaeus, Phys. Rev. E 45 (1992) 5339. Y. Qian, D. d’Humieres, P. Lallemand, Europhys. Lett. 17 (1992) 479. X. He, S. Chen, R. Zhang, J. Comput. Phys. 152 (1999) 642–663.