LBM simulation of bubble behavior and CHF prediction of pool boiling in narrow channels under different body forces

(1)

ADVANCED LEVEL, 30 ECTS

-STOCKHOLM BEIJING, 2018

(2)

(3)

LBM simulation of bubble behavior

and CHF prediction of pool boiling

in narrow channels under different body forces

Jing Yuan

Thesis Submitted to

Tsinghua University KTH Royal Institute of Technology In partial fulfilment of the requirement

for the degree of Master of Nuclear Energy and Nuclear Technology Engineering

In

Nuclear Energy Engineering

In partial fulfilment of the requirement for the degree of

Master of Science In

Engineering Physics

Co-supervisor: Associate Prof. Zhu Hongye

Co-supervisor: Associate Prof. Ma Weimin

TSINGHUA UNIVERSITY /INSTITUTE OF NUCLEAR AND NEW

ENERGY TECHONOLOGY

KTH-ROYAL INSTITUTE OF TECHNOLOGY /DEPARTMENT OF PHYSICS

UNDER THE COOPERATION AGREEMENT ON DUAL MASTER’S DEGREE PROGRAM IN NUCLEAR ENERGY RELATED DISCIPLINES

(4)

(5)

LBM simulering av bubbla beteende och CHF

förutsägelse av poolen kokar i smala kanaler under

olika kroppskrafter

Abstrakt

Med utvecklingen av stora fartyg och flytande kraftverk används kärnreaktorer alltmer under marina förhållanden. Som ett allvarligt hot mot reaktorns säkerhet är mekanismerna för värmeöverföringskrisen mer komplicerade i fartygsbaserade reaktorer än i land baserade reaktorer. Syftet med detta dokument är att undersöka beteendet (kärnbildning, tillväxt och avgång) hos bubblor under poolkokning i vertikala smala kanaler och förutse det kritiska värmeflödet (CHF).

(6)

riktning desto mindre är det längsgående längdavståndet för bubblans avgångspunkt. När tyngdriktningen var horisontell var längdriktningsavståndet för bubbelavgångspunkten noll. När kroppskraften var nära den horisontella riktningen ( 70 ) uppträdde den laterala förskjutningen av bubbelavståndspunkten. Baserat på resultaten av LBM-numerisk simulering, kombinerat med den nära väggen bubbla-trängningsmodellen, förutspådades CHF av poolen som kokar i slitsen under verkan av olika kroppskrafter. Relationerna mellan CHF och gravitation och kanalgeometri analyserades och optimeringsscheman föreslogs för rymdkonfigurationen av kärnreaktorbuntar. För en given geometri av kanalen ökade värdet av CHF-förutsägelse när kroppskraften ökade. När kroppskraften var konstant ökade värdet av CHF-prediktionen, eftersom slitsens bredd ökade. Under samma gravitation och slitsbredden var CHF för den rektangulära kanalen större än CHF i den cirkulära kanalen, så den rektangulära kanalen hade bättre värmeväxlingsförmåga.

I detta dokument uppnåddes reglerna för bubblaxpansion och avskiljning under olika kroppskrafter, och sedan uppskattades CHF i reaktorkärnan, vilken var lärande för reaktorens termiska utformning under marina förhållanden. Den använda forskningsmetoden har ett referensvärde för ytterligare undersökning av reaktorns kärnvärmeöverföringsförstörningsfenomen under dynamiska marina förhållanden.

(7)

Abstract

With the development of large ships and floating power stations, nuclear reactors are increasingly applied in marine conditions. As a serious threat for the safety of the reactors, the mechanisms of heat transfer crisis is more complicated in the ship-based reactors than that in land-based reactors. The objective of this paper is to investigate the behavior (nucleation, growth and departure) of bubbles during pool boiling in vertical narrow channels and predict the critical heat flux (CHF).

(8)

the lateral displacement of the bubble detachment point appeared. Based on the results of LBM numerical simulation, combined with the near-wall bubble crowding model, the CHF of the pool boiling in the slit under the action of different body forces was further predicted. The relationships between CHF and gravity and channel geometry were analyzed, and optimization schemes were proposed for the space configuration of nuclear reactor bundles. For a given geometry of the channel, the value of CHF prediction increased as the body force increased. When the body force remained constant, the value of CHF prediction increased as the width of the slit increased. Under the same gravity and the slit width, the CHF of the rectangular channel was larger than the CHF of the circular channel, so the rectangular channel had better heat exchange capability.

In this paper, the rules of bubble growth and detachment under different body forces were obtained, and then the CHF of the reactor core was estimated, which was instructive for the reactor thermal design under marine conditions. The research method used has a reference value for further exploration of reactor core heat transfer deterioration phenomenon under dynamic marine conditions.

Key words: Lattice Boltzmann Method; Marine Condition; Pool Boiling;

(9)

Chapter 1 Introduction

1.1 Background

Vapor-liquid two-phase flow is widely found in some power plants, such as cooling systems for boilers and nuclear reactor cores. According to the flow characteristics of liquid working fluids, boiling heat transfer can be divi ded into pool boiling and flow boiling. In pool boiling, the temperature difference and the turbulence caused by the nucleation, growth, and shedding of air bubbles on the heating surface are the main factors that cause the liquid working fluid to move. Since the 1920s and 1930s, the boiling heat transfer pioneers represented by Jacob[1], Fritz[2], and Nukiyama[3] have systematically studied pool boiling heat transfer. Using horizontally placed platinum wire as the heat source, Nukiyama [3] conducted an experimental study of pool boiling heat transfer and expressed the experimental result s as the relationship between heat flux and wall superheat, thus obtaining the well -known pool boiling heat transfer curve (Nukiyama boiling Curve). As the superheat increasing, there are a natural convection regime, a nucleate boiling regime, a transiting boiling regime, and a film boiling regime in turn. The nucleate boiling regime has the characteristics of small temperature difference and great heat transfer capability. Industrial applications are generally designed in this regime. In the pool boiling heat transfer curve, the peak of heat flux corresponding to the end of the nucleate boiling regime (ie, the starting point of the transition boiling regime) is called the critical heat flux (CHF). It is of great practical significance in industry practice. For a heating device that controls heat flow, once the heat flux exceeds CHF, the wall superheat will escalate sharply and may burn the heating surface. Therefore, the critical heat flux, also known as the burn point, is the highest limit of heat flux for safety.

(13)

etc. Since the launch of the United States "Naut ilus" nuclear submarine in 1954 more than 300 nuclear-powered ships have been cruising in the oceans. Using nuclear power for ships can greatly increase their capacities of endurance and concealment, and reduce the dependence on oxygen supply[4]. At the same time, offshore floating island nuclear power stations have attracted increasing attention in recent years. This small, por table nuclear power plant is a scaled-down version of the onshore nuclear power plant . It can be installed on a ship to provide energy supply safe and effective for remote islands, as well as provide electricity, energy and fresh water resources for oil and gas exploration platforms. Floating reactors are located in the sea and are not easily affected by earthquakes and tsunami. Even if an earthquake occurs, the seismic waves cannot be transmitted by seawater. Moreover, the ocean itself can also be used as an emergency cooler. Once extreme accidents occur, the seawater can be introduced into the floating reactor to prevent the core melting process and ensure reactor safety. Due to their small size, floating platforms can be towed to a special place for maintenance and handling.

However, due to marine conditions such as breeze and waves over the sea, additional movements such as swing, heaving, and inclining will occur on ships and offshore floating island nuclear power plants. The coolant will be subjected to inertial forces, and the effective body force will change periodically resulting in changes of flow characteristics, physical properties, and stability, and even affect the safety of nuclear reactors[4]. Therefore, in order to ensure the normal operation of nuclear power equipment, it is necessary to carry out the research on thermal hydraulic characteristics of nuclear reactors under ocean conditions.

(14)

capacity of the coolant change drastically due to the vaporization and even a boiling crisis occurs.

In recent years, the mesoscopic-scale lattice Boltzmann method (LBM) is receiving more and more attention. It shows great potential in the simulation of complex fluid systems. Conventional macroscopic numerical simulation methods are based on solving discrete control equations, while the lattice Boltzmann method can directly introduce forces among particles based on the theory of mesoscopic dynamics. Therefore, it has unique advantages of solving problems involving interface dynamics such as multiphase flow, multicomponent flow, etc, and complex boundaries such as flow in porous media. Based on the above background, this thesis will use the lattice Boltzmann method vapor phase model to study the two-phase flow characteristics under the non-inertial system.

1.2 Literature review

1.2.1 Non-inertial thermal hydraulics

The conversion process from a single-phase to two-phase exists in boiling flow. For steam generators and other thermal equipment, vap or bubbles will change the heat transfer characteristics and resistance characteristics. The oscillation in two-phase flow, especially the slug flow and plug flow will accelerate the fatigue damage of the heat transfer pipe, resulting in shortened service life of the heat transfer pipe[5], and the inertia forces in non-inertial systems are different for the gas-phase and liquid-phase forces because of the density difference between vapor and liquid phases[5]. The interface structure, flow characteristics and heat transfer characteristics are different from that in the traditional inertial system.

(15)

research is mature but mainly focuses on the overall parameters. The mechanism of thermal hydraulic characteristics the influence under ocean conditions does not be revealed.

John, et al.[11] find that the change of the gravitational field caused by the undulating motion will cause the periodic pulsation of the flow when the system is initially steady. The pulsation amplitude of the flow is linear with

/

max

a g . When the system is initially instable, the flow in the system is more

complicated due to undulating motivation. When the frequency of the flow instability and the change of the gravitational field are close to each other, the flow pulsation will be more severe due to resonance . The critical heat flux density will decrease under that condition.

(16)

and propose for the first time that the resonance of fluctuations in flow and flow instability caused by undulation should be avoided in the engineering design.

Zhang Jinling[17], Pang Fengge[18], Su Guanghui, et al.[19] conduct experimental and theoretical studies on the thermal-hydraulic characteristics of natural circulation under undulating conditions. It is shown that the gravitational field of the coolant system changes due to the undulation. As a result, different from forced circulation, significant periodic oscillation of the natural circulation flow rate appears with the undulation. Ishida, et al. [6] find that synchronous oscillations of the core flow and reactor power caused by undulations appear. The density wave oscillation in the system overlaps with the flow rate oscillation as the power increased. In addition, Ishida believes that the non-condensing gas in the upper space of the pressurizer helps to suppress the flow pulsation. Pendyala, et al.[20, 21] study the flow and heat transfer characteristics of forced circulation under undulating conditions. The studies show that change in the body force field causes a significant flow pulsation in the laminar region, which in turn increase the heat transfer. For the turbulent region, the flow rate does not change significantly under undulating conditions. The transfer characteristics in the turbulent flow are the same as in the stable state. Hong, et al. conduct experiments on the frequency of bubble detachment and the bubble detachment point under fluctuating conditions. The results indicate the change of bubble behavior mainly depends on whether the flow rate changes.

Aslam [22], Ma[23], Wu[24], Chen[25], Hirano[26], and Watanabe[27], et al. study the effects of the changes in high-frequency external force field caused by ground motions under seismic conditions on the thermal-hydraulic characteristic. The research content mainly focus on the effects of the free liquid surface sloshing caused by the high-frequency external force field on fluid-solid coupling, characteristics of boiling heat transfer and changes of reactivity.

(17)

relevant theoretical studies. Ma, Wu, and Chen, et al.[23-25] study the free surface sloshing of the sodium pool in the fast reactor in the high-frequency external force field caused by the earthquake, and obtain the dynamic response relationship between the free-surface sloshing of the sodium pool and the fluid-solid coupling of the equipment in the tank. Suppression program of liquid sloshing is proposed. Hirano and Tamakoshi, et al. [26] propose to equate the high-frequency external force field under seismic conditions with the internal body force of the fluid, and analyze by the modified TRAC-BF1 program the nuclear-thermal coupling phenomenon when the resonance of fluid surface sloshing in boiling water reactors under seismic conditions and thermal-hydraulic instabilities occurs. Watanabe, et al. [27] simulate the effect of free liquid surface sloshing in the boiling water reactor on the bubble behavior under seismic conditions using the dynamic mesh technology. It indicates that severe deformation of the gas-liquid interface due to the liquid surface sloshing.

Pang Fengge[18], Su Guanghui[19], Gao Yuzhen[28, 29], et al. establish a mathematics-physics model of the natural circulation capacity when the reactor is affected by swinging motions. It is demonstrated that the mechanism of the effect of the swing motion on the natural circulation. The additional inertial force applied to the coolant promotes or hinders the flow, affecting the original pressure drop. Yang Lan[30] and Jiang Chunlin[31], et al. analyze the effect of the swing motion on the operating condition of the coolant system with the constant coolant average temperature. It is found that the core power oscillates due to the swing motion and the total natural circulation capacity is reduced. Murata, et al.[32, 33] study the natural circulation heat transfer characteristics under swing conditions. The conclusions show that the heat transfer under swing conditions is affected by both the ship motion and the coolant flow.

1.2.2 Bubble dynamics and boiling mechanism

(18)

sub-cooled boiling regime and the fully-developed sub-cooled boiling regime. It is also called the demarcation point between the higher sub-cooled boiling regime and the lower sub-cooled boiling regime. At this point, the void fraction begins to increase. The characteristics of the two-phase flow need to be considered.

There are many achievements in the research on ONB. Since the 1960s, scholars from various countries have conducted in-depth research on this topic. The researches of Bergles & Rohsenow, Davis & Anderson et al are representative. Bergles & Rohsenow[34] study the starting point of sub-cooled boiling in a pressure range of 0.1 to 13.6 MPa within an annular channel heated by inner tube with an inner tube diameter of 1.6383 mm and an outer tube diameter of 5 mm quartz tube. The fitting of the empirical relationship between sub-cooled boiling starting point and heat flux density in the steam-water system is given as:

0.0234

1,156 2.41/

5.3 [1.8( ) ] p

ONB w s ONB

q  p T T (1-1)

Davis and Anderson[35] derive the relationship between forced convection starting point and heat flux density in large-diameter pipes using analytical methods. 2 1 ( / 8 )( ) ONB sf fg sg s w s q  K H  C T T T (1-2) where, H is the latent heat of vaporization, _fg K is the saturated liquid _sf

thermal conductivity, _sg is the saturated vapor density, and  is the surface tension coefficient, C is the ratio of the height of the bubble(i.e. the ₁

distance between the top of bubble and the wall) to the radius of curvature of the bubble, andT is the temperature.

Bowring R.W[36] modifies the above equation using Pr formula in terms

of various fluids based on formula (1-2).

2 1

( / 8 )[( ) / Pr ]

ONB sf fg sg s w s ONB l

(19)

Saha[37] believes that the point of net vapor generation must meet both the thermal and fluid dynamic limits and be separated by Peclet number Pe . Equation (1-4) is derived to predict the thermodynamic equilibrium dryness at the point of sub-cooled boiling net vapor generation in Saha model.

0.0022( ) Pe 70000 153.8( ) Pe 70000 e pf FDB f fg FDB fg qD C X k i q X Gi         _{ } _  (1-4)

where q is the heating load. G is the mass flow rate. De is the equivalent diameter of the pipe. Cpf, kf, and ifg are the specific heat, thermal conductivity,

and latent heat of vaporization, respectively. Pe is the Peclet number. The

Saha model is considered to be a good model for predicting the net vapor

generation point and is suitable for working fluid such as water and refrigerants.

Levy[38] proposes that the force balance of the bubble at the nucleation point and the temperature distribution of the single-phase liquid away from the wall need to be considered for determining the bubble detachment point. The balance of forces is constituted by buoyancy, surface tension, and wall shear force. The bubble is not hemispherical but rather spherical . The thickness of the boundary layer is assumed as the height of the monolayer bubble. The temperature of the bubble detachment point can be determined by considering the force on the bubble and the temperature distribution of the single -phase liquid away from the wall. However, a large error is verified by many researchers when applying Levy model at low flow rates.

Sun Qi[39] modifies the levy model and points out that buoyancy cannot be neglected at low flow rates. Sun Qi [39]'s sub-cooled boiling NVG calculation model is given as follows.

(20)





0.5 _-1.2





0.41 ' Re Prf f i fg g f B k h r r D _  (1-6) f i B fg f u D Re      (1-7) 2 3 ( ) 4 B f g e f D G g D      _   (1-8) 2 ( ) 1.41[ g f g ] gj f u       (1-9)

where, kf ( kW/m2﹒K) is the thermal conductivity of liquid phase. μf( kg/m﹒

s) is the viscosity of liquid phase. σ(N/m) is the surface tension. g( m/s2) is the gravitational acceleration. Prf is the liquid Prandtl number.

1.3 Contents and methods

1.3.1 Contents

(21)

Figure 1.1 The narrow channel in the rob bundle space

Therefore, the purpose of this thesis is to study the bubble dynamics in the space of rob bundle under the non-inertial frame and the heat transfer mechanism of the pool boiling. The nucleation, growth, detachment of the bubbles in pool boiling under oceanic and seismic conditions were discussed. The transform between nucleate boiling and the membrane boiling, heat transfer coefficient and critical heat flux were studied. Specifically, the following contents were included: a numerical model for the pool boiling bubble dynamics in the rob bundle space under non-inertial frame was established; the effects of hull motions such as inclining, swing and undulations on the bubble motion and heat transfer were studied; CHF for pool boiling in the rob bundle space was predicted; The effect of the rob bundles space configuration and the position of the partition on bubble motion and heat transfer was investigated and an optimization scheme was proposed.

1.3.2 Methods

(22)

Among the current interface processing methods, the most important methods are VOF [40] (Volume of Fluid) method and Level Set method. VOF method was first proposed by Hirt and Nichols in 1981[41]. The basic principle of VOF method is to determine the free surface by studying the ratio function F of fluid volume to the cell volume instead of the particle motion in the free surface. Then the change of fluid is tracked. VOF method can deal with strong non-linear phenomena such as reentry of free surfaces. It requires short calculation time and low storage, but it is a bit tedious when dealing with changes in F, and there are certain human factors. The VOF method constructs and tracks the free surface based on the volume ratio function F. If F=1, the unit is all occupied by the specified phase fluid. If F=0, the unit is no specified phase fluid unit. When 0 <F <1, the unit is called the interface unit. Assuming any point in the flow field(x, y), the function f(x , y , t) is defined as follows.

[42-44]

1 ( , )

( , , )

0 ( , )

fluid particles of this phase exist at x y f x y t

fluid particles of this phase do not exist at x y

  

 (1-10)

The conservative form of the transmission equation is expressed as:

0 f uf f t x y   _ _ _    (1-11)

Using the staggered grid shown in Figure 1.2, the shadows represent the liquid portion, and the time is in the first-order difference format. The transmission equation differential is shown as follows.

1 1 1 1 1 , , , , , , 2 2 2 2 =0 n n i j i j i j i j i j i j i i F F F F F F t x y      _ _  _ _  _      (1-12) Where ₁ , 2 i j F   and , 1 2 i j F 

 are the fluxes through the grid cell boundaries

(23)

1 1 , , 2 2 1 1 , , 2 2 1 ( ) 1 ( ) j i y i j i j j x i j i j i F dt uf dy t y F dt f dx t x           _  _{ }    _    

 

(1-13)

Figure 1.2 Variables on grid cells

Level Set method proposed by Osher et al.[45] in 1988 defines the material interface that moves with time as the zero isosurface of the distance function

( , )x t

 . ( , )x t satisfies certain equations. At each time t, as long as the value of the function ( , )x t is obtained, the location of its isosurface, that is, the location of the material interface can be determined. Level Set method has advantages in the interface description of the merger, crossover, and fragmentation phenomenon. However, it has a rounding effect on the interface.

[46]

The constructor ( , )x t is such that at any moment, the free interface( )t

(24)

At any moment t , for any point x on the free interface ( )t , ( , )x t 0   . · 0, d d dx V V dt dt dt  _ _{  }_ _ (1-14) This equation is called Level Set equation.

The lattice Boltzmann method (LBM, Lattice Boltzmann Method) is a new numerical method for computational fluid dynamics that evolved from lattice gas automata[47, 48]. In 1970's Hardy et al. propose the first lattice gas automata model in history. In the lattice gas automata, the fluid particles exist on discrete grid nodes and move along the grid lines. All the particles perform collision and migration movements in accordance with certain rules. At the boundary, the lattice gas automata simply uses the rebound format to facilitate the calculation. At the same time, the presence or absence of particles on the grid line is respectively described by 1 or 0. Thus, the lattice gas automaton can be unconditionally stable. Therefore, this model has many advantages such as stability, simple boundary treatment, and parallelism.[49-51] The lattice Boltzmann method gestated on the basis of it retains its advantages while overcoming many deficiencies, such as the inability to satisfy the Galileo invariance, the dependence of the fluid state equation on the macro vel ocity, the local noise, and high computational capabilities required by the collision operator.

(25)

Galilei invariance. The fluid state equation is decoupled from the flow velocity. After 20 years of development, the lattice Boltzmann method has been widely used in many fields such as micro/nano -scale flow[55], porous media flow[56], multi-phase multi-mass flow[57, 58], and non-Newtonian flow[59].

Compared with traditional computational fluid dynamics methods (such as finite element method, finite difference method, etc.), LBM has the following advantages[60]:

(1) Simple algorithm. Complex nonlinear macroscopic phenomena are simulated by a simple linear operation and a relaxation process.

(2) Ability to handle complex boundary conditions.

(3) The pressure in the lattice Boltzmann method can be directly solved by the state equation.

(4) Easy to program and the calculation before and after processing is also simple.

(5) Massive parallelism.

(26)

Chapter 2 LBM modelling of two-phase flow and boiling

The numerical study of vapor-liquid two-phase flow has been one of the hot issues in computational fluid dynamics (CFD). The lattice Boltzmann method has unique advantages in solving problems involving interface dynamics due to the mesoscopic scale and particle properties. Therefore, the multi -phase flow numerical study based on the lattice Boltzmann method has received extensive attention from scholars in the past two decades. .

At the macro level, fluids are assumed to be continuous. The motion of fluid satisfies the laws of conservation of mass, momentum, and energy. The motion law is described by Euler's equations, Navier-Stokes equations, and so on. The existing simulation methods (such as finite difference method, finite element method, finite analysis method, etc.) are mostly taken to discrete the nonlinear partial differential Euler equations, Navier-Stokes equations. Then calculate the discrete equations by computer. Due to the large scale of these CFD methods, it is often not suitable for problems tracking the phase interface of multi-phase or multi-component fluids and obtaining information on bubble nucleation, growth, and detachment. Considering that at the microscopic level, the fluid is composed of a large number of discrete molecules, and the motion characteristics of the molecules are influenced by the intermolecular interaction forces and external forces. All macroscopic thermodynamic properties of the fluid and the motion law are realized as the irregular molecular thermal motion at the microscopic level. Based on this theoretical basis, a very direct idea is to simulate the motion law of each discrete molecule and derive the flow characteristics of the system through statistical rules. Due to the limited computing capacity, the simulations performed by this method are extremely limited on both space and time scales. At the same time, the macroscopic properties of the fluid are not sensitive to the details of the irregular thermal motion of each specific molecule.

(27)

is to replace complicated and variable macroscopic phenomena with simple rules of microscopic particle motion. The movement of particles at each time step consists of two steps: flow and collision. When multiple particles reach the same grid node, they interact and redistribute according to the collision rules. Each particle moves in its speed direction to the nearest grid node according to the flow rules. [61]

The colored model, pseudo-potential model, and free energy model develop to describe multi-phase flow. The colored model is proposed by Gunstensen et al.[62] in 1991, based on a two-component lattice gas automata model, and improved by Grunau et al.[63] in 1993. Different component fluids are distinguished in different colors and the interactions are described by color gradients. In this model, the surface tension is anisotropic . So false flow phenomenon occurs near the phase interface. Large amount of calculation is required. In 1995, Swift et al.[64] proceed from the theory of free energy, using the inverse function of free energy to construct an equilibrium distribution function, and introduced the thermodynamic pressure tensor of nonideal fluid to satisfy energy conservation. The main drawback of the free energy model is that Galilei invariance is not satisfied.

In 1993, Shan et al.[47] propose the pseudopotential model and construct an interaction potential function to describe non-local interaction forces between particles to control the state equation. Yang [65] proves through quantitative experiments that the simulation results of bubble growth of in pseudo-potential model are very close to the experimental observations. Therefore, the two-phase flow simulation of the lattice Boltzmann method in this paper is based on the pseudo-potential model. The two-phase interactions are described by introducing the lattice Boltzmann particle distribution function evolution equation with the interparticle force.

2.1 Lattice Boltzmann equation

(28)

proposed the two-phase flow model of the pseudo-lattice Boltzmann method in 1993, the model has achieved vigorous development and wide application.

In the lattice Boltzmann method, the fluid state is characterized and described by the fluid particle distribution function[66].

2.1.1 Boltzmann Equation and Maxwell Distribution The Boltzmann equation can be expressed as follows:

 

· f f f J f t         r a  (2-1)

where, f( , , )r  t is the particle distribution function at r with the velocity  at time t . J f is the collision term and represents the variation of the ( ) number of molecules generated by the collision in integral form.

In a vector field, a particle at r with the velocity  at time t , if there is

no collision in the time interval dt , will move to rdr at time tdtand the velocity is adt . The actual particle distribution function with the velocity adt at rdr is f(rdr,adt t d d, ) r  . The change of the distribution function is caused by a collision. So the following equation is given.



, ,





, ,



 

f rdr  adt tdt d dr  f r  t drd J f drd (2-2)

The Maxwell distribution is an equilibrium distribution. The Maxwell distribution is the most basic physical law in the theory of gas dynamics. It specifies the most probable velocity distribution of the gas velocity in a stationary equilibrium state. It has been proved by measurement using molecular beam techniques. For D-dimensional space, the Maxwell distribution can be expressed as:









2 2 2 2 1 eq g D g f n e R xp R T T   _          u (2-4)

where , n is the molecular number density, u is the macroscopic flow velocity of the fluid, R is the gas constant, and T is the gas kinetic g

temperature.

(29)

2

3 1 1

2R Tg n



2c fd (2-5) where , c is the thermal movement speed of the molecule , c  u。 2.1.2 BGK approximation

The collision term J f is a non-linear term in the distribution ( ) function and is related to the intermolecular force. The existence of the collision term brings great difficulties to solving the Boltzmann equation. The numerical integra tion of each particle by the computer will introduce a huge amount of calculation. The macroscopic properties of the fluid, such as the density field and the velocity field, are insensitive to the details of the collision betwe en the single fluid microscopic particl e. In 1954 , Bhatnagar , et al.[ 6 7 ] propose to replace the collision term J f with a simple operat or ( ) _f to simplify the calculation of the collision function .

Assume the collision effect is to chan ge the particle distribution function f so that it tends to the Maxwell equilibrium distribution

eq

f . u_{and T in the} eq

f equations are the macroscopic and macroscopic temperatures of the gas, which are obtained from the velocity distribution function f . Because the equilibrium distribution function changes with time and space, it is called a local equilibrium distribution function. Assume that the rate of change is proportional to the difference between feqand f and the proportionality factor is  . Then we can introduce operators as follows.



, ,





, ,



eq

f f t f t 

  _ r   r  _ (2-6) Substitut e this operator into the Boltzmann equation and i ntroduce the relaxation time 0 which represents the average time interval

between two collisions.

0

1 



 (2-7)

(30)









0 1 · , , eq , , f f f f t f t t            _ _ _ a r r _(2-8)

Multiplying the distribution function by the mass of the particle m , it represents the density distribution of the particle. The corresponding Maxwell equilibrium distribution becomes :









2 2 2 2 1 g eq D g f e R T R p T x    _          u (2-9)

2.1.3 Lattice Boltzmann equation

Lattice Boltzmann equation is a special discrete form of the Boltzmann -BGK equation. Th e discretizatio n include s velocity discretization, time discretization and spatial discretization. The particle velocity  is simplified to a finite dimensional velocity space



e e0, ,...,1 eN



, and the continuous distribution function f is also

discretized into



f0,f1,..., fN



, where fi  fi( , , )r ei t , i0,1,... ,N .The

speed-discrete Boltzmann equation can be expressed as:

0 1 · eq i i i i i i f e f f f F t          (2-10)

where , fieq is the local equilibrium distribution function in the discrete

velocity space . F is the external force term in the discrete velocity i

space, which is the projection of the external force termaf in the discrete velocity space .





i _i

F  a_ f (2-11) A second-order Taylor expansion of the local equilibrium distribution function fieq in the discrete velocity space with low flow velocity is

(31)

where, _i is the weight coefficient.

Discretize equation (2-10) both in space and time. Then integrate along its characteristic line. The integral term in the right is approximated by a first-order rectangular method. The evolution equation of the fluid particle distribution function using the BGK collision operator can be obtained:





 

1

 

, , , eq , i i i i i i f t t t f t f t f t F t         _  _  r e r r r (2-13) where, 0 dt 

  is the dimensionless relaxation time. Fit is the volume

force term.

The local macro-density  and local macro-velocity u in the local equilibrium distribution function can be expressed as:

 

, i i f t 



r _(2-14)

 

, i i i f t u



e r _(2-15)

In the lattice Boltzma nn method, the temporal and spatial discretizations are related by discrete speeds and are not independent of each other, such as in the direction x , xeixt. In physical space, this

means that the particle mov es from one node to another in a time step, and collid es with another particle at this node. Based on this characteristic, the lattice Boltzmann method divides the particle motion into two simple processes of migration and collision, so it has good parallel characteristics and capabilities to deal with complex boundary.

Using the Chapman -Enskog expansion[ 6 8 ], the corresponding macroscopic equations for the lattice Boltzmann basic model can be derived. The derivation results show that the lattice Boltzmann method can be restored to the standard Navier -Stokes equations for low flow (low Mach number) fluids.

2.1.4 Energy equation and source term

(32)

considered as a passive scalar following the fluid moti on and satisfy the convection-diffusion equation. So the temperature can be regarded as an additional component of the fluid system. In 1997, Shan propose a two-component model to simulate the heat flow problem based on the pseudopotential model[69]. The component I model is a flow model and the component II is a temperature field model. This model is called a passive scalar model. In the passive scalar model, the temperature field is given by the following equation.

 

· ·( ) T T T t    _{ } _{   }  u (2-16) where,  is the thermal diffusion coefficient defined by formula (2-17).

2 1 ( ) 2 s T t c      (2-17)

where, Tis the relaxation time of the temperature distribution function.  is

the source term of the energy equation that is related to the vapor -liquid phase transition. The temperature distribution function evolution equation is given as follows. 1 ( , ) ( , ) ( , ) eq( , ) i i t t i i i t i T T  t  T t T t T t           _  _ r e r r r _(2-18)

The equilibrium temperature distribution function is given by equation (2-19).

 

2 ₂ 2 4 2 · · ( , ) 1 2 2 i eq i i i s s s T t T c c c              e u e u u r (2-19)

The definition of the local macro temperature T is as follows.

i i

T 



T _(2-20)

Excluding viscous heat dissipation, the source term  is derived from the entropy balance equation as follows. The entropy balance equation can be written as:

·( )

dS

T T

dt

    (2-21)

(33)

2 1 1 p p p TdS c T d c dT T d c dT T d T d T T T        _ _               _ _   _ _ _{ }  _ _            (2-22)

Substitute the derivation result of equation (2-22) into equation (2-21). Then the following equation is derived.

2 · dT T p d T dt c_ c_ T _ dt          _  _ _ _      (2-23) Due to

 

· dT T T dt t      u (2-24) the following equation is obtained from equations (2-23) and (2-24).

 

·





2 T T p d T T T t c_ T _ dt     _{ } _{   }   _{ }    u   u (2-25)

Compare expressions (2-16) and (2-25), and the expression of the source term can be derived. 2 T p d T c_ T _ dt        _ _      u (2-26)

In equation (2-26), the full-derivative term for density to time appears. In numerical calculations, this term introduces a large amount of calculations. So according to the continuity equation

 

· u 0 t  _  _{ } _  (2-27) A further simplified energy source term expression can be obtained.

1 1 p T c_ T _    __ _     _ _         u (2-28)

2.2 DdQm discrete model

(34)

Boltzmann model can be restored to the macroscopic equation that need to be solved . Therefore, in the lattice Boltzmann method, the construct ing of discrete velocity model is very important. Some physical quantit ies will not be able to satisfy the corresponding conservation laws if the number of discrete speeds is insufficient . Unnecessary calculations will be introduced if the number of discrete speeds is excessive , resulting in a decrease in computational effici ency. In 1992, Qian , et al.[ 5 4 ] propose the DdQm series model, where d is the spatial dimension and m is the number of discrete spee ds. This series of models has become a classic basic model in the development and application of the lattice Boltzmann method.

For 2D and 3D simulations, the most common discrete velocity models are the D2Q9 and D3Q19 models. The discrete speeds are shown in Figure 2.1 and Figure 2.2, respectively .

(35)

e1 e2 e3 e₄ e5 e6 e7 e8 e9 e10 e11 e12 _e 13 e14 e15 e16 e17 e18 Figure 2.2 D3Q19 model

The DdQm model uses the following equilibrium distribution functions.

 

2 2 2 2 4 · · 1 2 2 i s i s eq i i s c c f c              e u e u u (2-29)

where , c is the lattice speed of sound. s

The speed matrix in the D2Q9 model is defined as follows.

0 1 0 1 0 1 1 1 1 0 0 1 0 1 1 1 1 1 E    c       _ _ (2-30) The speed matrix in the D3Q19 model is defined as follows.

0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 E c                         (2-31)

where , c x/ t,x and t are respectively the grid step and the time

step. For the convenience of calculation, the grid step len gth is usually the same and is equal to the time step, ie. x y t. In the D2Q9

model and the D3Q19 model, c2_s c2 3. The kinematic viscosity coefficient is

(36)

For the D2Q9 and D3Q19 models, the weight coefficients i are

given in Table 2.1.

Table 2.1 Selection of Weight Coefficients of Equilibrium Distribution Functions for D2Q9 and D3Q19 Models i  ei D2Q9 D3Q19 0  0 4/9 1/3 1  1 1/9 1/18 2  ₂ _1/36 _1/36

In t his paper simulat ion is processed in two-dimensional space using the D2Q9 model.

2.3 The conversion of grid units and physical units

When performing the numerical simulation of the actual physical model, the unit of the actual physical quantity n eeds to be c onsidered. Usually, dimensionless method for physical parameters is implemented in order to facilitate the program code . The heat transfer criteria number of flowing need to be consistent before and after the non-dimensionalization. In the lattice Boltzman n method, we also make similar considerations to convert the actual physical units into grid units.

(37)

density, and u the reference speed. Then conversion of paramete rs r

under the grid unit and the actual physical unit is defined in Table 2.2. Table 2.2 Conversion of basic physical parameters under grid unit and actual unit The parameters under the

grid unit

The parameters under the

physical unit Conversion relationship

L L r L L L   _ r     t  _t t r t r u L    s c c_s s s r c c u   _ r r L u   

For specific problems, the L,  ,  andc used in the simulation are s

known. The actual physical density , velocity cs, and kinematic viscosity

coefficient  are also available. Therefore, the conversion between the grid unit and the actual physical unit is completed through the relationship shown in Table 2.2.

2.4 Boundary conditions in LBM

(38)

determine the value of the distribution function at the boundary node through the known macroscopic boundary conditions, such as density, velocity, pressure, temperature, and force conditions of the flow, so as to perform the next calculation.

The temperature distribution function in D2Q9 model is shown in Figure 2.3.

T1 T5 T2 T6 T3 T7 T4 T8 T0

Figure 2.3 Temperature Distribution Function in D2Q9 Model

The five boundary conditions involved in the lattice Boltzmann method in this paper are described below.

2.4.1 Fixed wall temperature boundary condition

In this paper, the left boundary heating point of the calculation area is defined as the fixed wall temperature boundary condition. Under the fixed wall temperature boundary conditions, T1, T5, and T8 are unknown temperature

distribution functions. Assume that they are equal to the equilibrium distribution function. The unknown local macro temperature is represented by a variableT . s

 

2

 

4 2 22





· · , 1 , 1,5,8 2 2 i eq i i i i s s s s T T t T i c c c                e u e u u r (2-33)

Add the three distribution functions of equation (2-33). Due to the symmetry of the discrete speeds the following equation can be obtained.





1 8 2 5 1 1 3 3 6 s x x T T T T     u  u (2-34)

where, u is the velocity component perpendicular to the left boundary. Since x

(39)

according to the definition of the temperature distribution function and the local macro temperature the following equation is given.

h i

i

T 



T _(2-35)

The expressions of unknown value T can be derived through simultaneous s

equations (2-34) and (2-35).



0 2 3 4 6 7



2 6 1 3 _x 3 _x s h T  T      T T T T T T  u  u (2-36)

The unknown temperature distribution functions T1, T5, and T8 can be

determined.

2.4.2 Adiabatic temperature boundary condition

The adiabatic temperature boundary condition is used on the left boundary in the calculation area except the heating point.

0, 0 j T x    (2-37)

The second-order precision finite difference scheme is ued:



1, 2, 0,



0, 4 j 3 j / 2 j x j T T T T x       (2-38)

Since the internal temperature field is known, the wall tempe rature under adiabatic conditions can be determined according to equations (2 -37) and (2-38).

2.4.3 Periodic boundary condition

(40)

Layers of virtual grids



mN_y1;n1, 2, ,N_x



and



m0;n1, 2, ,N_x



are added to the upper and lower boundaries of the grid in the calculation area



n1, 2, ,N_x; m 1, 2, ,N_y



. According to the symmetry of the grid and the discrete velocity, the particle distribution function under the periodic boundary conditions can be expressed as:

 









 

, 0 , 2,5, 6 , 1 ,1 4, 7,8 i i i y y i f n f n N i f n N f n i          (2-39)

2.4.4 Outflow boundary condition

The right boundary of the calculation area adopts outflow boundary conditions. The equation that satisfies the right outflow boundary condition can be written as:

0, _x f f U x N t x        (2-40) where, U is the characteristic speed of the fluid. Use the first -order accuracy difference format for equation (2-40).



_x,





x, , t





x, ,



t f N m t f N m t f N m t        (2-41)



x,





x, , t





x 1, , t



x f N m t f N m t f N m x           (2-42)

Equations (2-41) and (2-42) can be used to give the distribution function evolution equation on the outflow boundary:



, ,





, ,





1, ,



1 i t i t t x i x x f N m t f N m t m t f N             (2-43)

(41)

2.4.5 Speed boundary condition

For the left boundary particle distribution function evolution equation, the velocity boundary conditions must also be set. The bottom boundary velocity is set to 0. f1, f5, and f8 are unknown particle distribution functions. According

to ( , ) i i f t 



r _(2-44) ( , ) 0 i i i f t u



e r  _(2-45)

the following equations are given.







2



6 8 1 5 8 0 3 4 7 5 7 5 2 4 6 1 8 3 6 7 f f f f f f f f f f f f f f f f f f f f f                           (2-46)

Apply the bounce format to the discrete function that is vertical to the wall to determine the unknown.

1 3

f  f (2-47)

Unknown particle distribution functions f2, f5, f6 and density  can be

obtained from equations (2-46) and (2-47):













1 3 2 4 8 6 2 4 5 7 0 2 4 3 6 7 1 2 1 2 2 f f f f f f f f f f f f f f f f                        (2-48)

2.5 Interparticle force

(42)

phases depends on the introduced interparticle force terms. The expression of force between particles is as follows.

 

0

( )

int  c g x

F x x (2-49) It can also be rewritten as:

 

2 0 ( ) / 2 int  c g x F x (2-50) where, ( )x is called effective mass and g is the coefficient of interaction between particles. c depends on the lattice structure. In the D2Q9 model and 0

the D3Q19 model, c0 6.0.

The fluid state equation corresponding to the e xpression of the interaction force between the particles is represented by the equation (2 -51). The state equation combines with the two-phase flow model of the lattice Boltzmann method through the effective mass.

 

2 2 0 2 s c pc  g_  _ (2-51)

Equation (2-50) can also be rewritten as:

 



2



0 c 2 p s c g      (2-52)

In this paper, the Peng-Robinson state equation (P-R EOS) is selected to determine the pressure p. The P-R state equation is given as tfollows.

2 2 2 ( ) 1 1 2 RT a T p b b b            (2-53)

In this formula, ( )T  [1 (0.37464 1.54226 0.269922)(1 T T/ _c)]2, where  is the acentric factor, which is related to the working medium. Critical parameters T and c p are obtained by equations (2-54) and (2-55). c

2 2

0.45724 c / c

a R T p (2-54)

0.0778 _c/ _c

(43)

Parameters related to the P-R state equation are shown in Table 2.3. Table 2.3 Parameter Values in P-R State Equation

a b R 

2/49 2/21 1 0.344

Combining equations (2-49) and (2-50), the new interparticle force forms are used in this paper:

 

2

 

0 0

( ) (1 ) / 2

int   c g x   c g x

F x x (2-56) where  is the weight factor, and the corresponding value of P -R equation is 1.16. The numerical discrete form of equation (2-56) is as follows.

 





 





1







2

 



( ) 2 , , int x x G  G              



 



 F x x x x x x x x x x x x _(2-57)

where, G x x( , ) is the Green function, which reflects the interaction among fluid particles. Generally only the influence of neighboring nodes is considered. The negative value of Green's function represents that the particles are attracted to each other. According to experience, the selection of the Green's function is shown in equation (2-58).





1 2 , 1 , , 2 0, otherwise G g g   _      _{ }    x x x x x x (2-58)

In the D2Q9 model, g12 ,g g2 g/ 2 . In the D3Q19 model,

1 , 2 / 2

g g g g .

2.6 Method to implement forces

(44)

The idea of introducing force in this paper is to add the discrete force expression directly to the right side of the particle distribution function evolution equation.

In equation (2-13), the force term is given as follows.

 

,

 

, ?

V i t

f t F t 

 r  r (2-59) where, Fi

 

r,t is a numerical discrete form of force. In this paper, the

discrete force format proposed by Guo [70] in 2002 is implemented. Considering the discrete lattice effect and the contribution of body force to momentum, the discrete format is given as follows.





2 4 1 1 2 i i i i i s s c -c       _   _       e U e U F e F (2-60)

where, F is the force, U is the real speed of the fluid. According to the conservation of momentum before and after the collision, the real velocity expression of the fluid is given as:

2

t



U u F (2-61)

The speed in the equilibrium distribution function should also be replaced by the following equation:

2 t i i i f  u 



e  _{F (2-62)} Comparing formula (2-61) and (2-62) it is found that in Guo's discrete format of force u is equal to the real fluid velocity U.

In the two-phase flow model, the force is given as: ( ) ( )

int g

 

F F x _{F x (2-63)}

where Fint( )x is the interaction force between particles. F x is body force g( )

given as:

( ) ( ( ) )

g x   x ave g

(45)

Where, g is the gravitational acceleration. ave is the average fluid density

within the entire fluid calculation area.

2.7 Model Validation

2.7.1 Equation of State

As described in Section 2.5, many researchers use the Peng-Robinson equation to describe the thermodynamic state of two -phase fluids in the LBM calculation. In order to evaluate the state equation and the simulation accuracy of the LBM calculation model for the thermodynamic state of the two -phase fluid, the saturated vapor and liquid saturation densities at different saturation temperatures were numerically simulated. The results were compeared with the IAPWS-IF-97 water properties. In the calculation, a two -dimensional uniform grid with 150×150 nodes was adopted, and the boundary was no-slip and adiabatic conditions. At the initial time, superheated steam with a radius of 10 to 40 grids was set in the center of the calculation domain, and the rest was sub-cooled water. The fluid temperature was evenly distributed. After the calculations began, under the effect of interparticle forces, the two-phase density and temperature constantly changed and reached a steady state after more than 50,000 time steps, forming a clear bubble as shown in Figure 2.4. At this point, steam and liquid were in a stable state of saturation.

By adjusting the initial temperature of the fluid and the vapor and liquid densities, the saturation temperature and saturation density under different pressures are calculated to form the saturation curve of water-steam (T-ρ) as shown in Figure 2.5. The critical temperature T of the Peng-Robinson c

equation was 0.07292, and the critical density c was 2.995. In general the

(46)

Figure 2.4 Bubble and Surrounding Density Distributions in Saturation

Figure 2.5 Numerical simulation of water-saturation saturation curve 2.7.2 Horizontal heated wall bubble departure diameter

The pool boiling problem of horizontally heated wall has been studied by many scholars. To further evaluate the simulation accuracy of the boiling process by the LBM model used in this paper, numerically simulation of the pool boiling on horizontally heated wall was conducted to observe the nucleation, growth, and detachment processes of the bubble, and analyze the effect of body force and contact angle on the detached diameter of the bubble. Compare the results with experimental correlation.