- Email: [email protected]
Simulation of shock-interface interaction using a lattice Boltzmann model
ELSEVIER
Nuclear Engineeringand Design 155 (1995) 67-71
Simulation of shock-interface interaction using a lattice Boltzmann model H. Ohashi, Y. Chert, M. A k i y a m a Department of Quantum Engineeringand Systems Science, Universityof Tokyo, Bunkyo-ku, Tokyo 113,Jat~
Abstract For simulation of fundamental processes of vapor explosion, a novel fluid dynamic simulation method is developed using the lattice Boltzmann model. We incorporate into the model capabilities to deal with comptess/b~ fluid flow and two-phase flow. By numerical simulation, we demonstrate fundanmntal c h a m t i c s of this medei and apply it to simulation of interaction betw~n shock wave and two-phase interfaces.
1. Introduction Simulation plays an important role in vapor explosion studies, in two possible ways: one is to reproduce the whole proczss of the vapor explosion to evaluate its effects, and the other is to simulate fundamental processes of the vapor explosion to understand their mechanisms (Akiyama, 1993). Formidable difficulties are encountered in the numerical work of the vapor explosion because it includes complex geometry, complex physics, phase changes and moving boundaries. In addition to that, conventional fluid descriptions require many constitutive models on which we hardly obtain information and validation from experimental results. So novel simulation methods based on mechanistic or physical backgrounds are need~a, m particular to elucidate the fundamental mechanisms by simulation. We have developed totally new algorithms to analyze fimdamental processes of the vapor explosion and applied them to simulation of pressure
shock, phase separation and interaction of pressure shock and two-pha.~ interface. These constitute key processes which caus¢ the coherent explosion. The lattice Boltzmann equation (LBE) mmthod is a recently developed numerical scheme to s i m ~ late fluid flow (McNamara, 1988). As a derivative of lattice gas automata (LGA) (Hasslacher, 1986), the LBE method solves macroscopk: c o m p ~ glow by dealing with the underlying micro-wotki. It describes the dynamics of partich) motion on a discrete lattice in terms of ensemble averaged particle distribution functions. Tim adv~mtages of adoption of the LBE model are two-fold in cornpar/son with the LGA model. First, the whit~ noise of the LGA model, which comes from the statistical nature of the automata, is completeiy removed. Second, the equations of ~ dynam/cs, the Navier-Stokes equations, are completely mcovered. Although the interua! stability ~ bib operation are no longer ensured in the LBE model, the simple implementation arid the ab~'ty
0029-5493/95/$09.50 © 1995 ElsevierScienceS.A. All rights reserved SSDI 0029-5493(94)00869-Z
68
H. Ohashi el at. / Nuclear Engineering and Design 155 0995) 67- 71
of accurate simulation, as well as the inheren.~ parallel nature, make the mt~el viable in the computational fluid dynamics (CFD) field. We have developed a class of multi-speed LBE model for compressible ~ow. We added one degree of freedom into the sound speed of modeled fltdd with the use of rest panicles and a particle reservoir. The sound speed is formalized to be dependent on thermodynamic variables in the way the perfect gas does in adiabatic cases. We also incorporated the capability of simulating multiphase immiscible fittids in the LBE model by introduction of multi-color particles and an interparticle potential. We will give results of shock wave simulation, phase separation and interaction of shock with two-phase interfaces.
2. LBE simdation The LBE simulation of fluid flow has become rapidly popular among the researchers of lattice hydrodynamics because it is noise-free and so flexible that the physics of real fluid can be fully recovered. Furthermore, the LBE models for fluid dynamics are extremely suitable for modern supereomputing because the calculation procedure is inherently parallel. This realization is due to some fundamental characteristics of the discrete microscopic model: the intrinsic stability, the easiness in deaiing with compficatexi geometric boundaries and the fast and local operations adequate for massively parallel processing. The lattice Boltzmann equation is expressed in the ensemble averaged form of the micro.dynamic equation for particles moving on a spatial lattice. It is given as follows Jv~,(x + c~, t + 1) -- N~,(x, t) = ~ ( N )
(1)
with the Boltzmann approximation applied to the collision term. Here Nm is the particle distribution and A~ is the collision term on the ith lattice node. Thus, the LBE exactly describes the twostep motion of particles in the micro-world. The first step is the streaming period during which the local particle distribution is transported to one of its nearc, t neighbors. The second step is the collision period during which the particle distribution
approaches its equilibrium value that depends locally on the macroscopic density and momentum. Recently a lattice BGK model was proposed to recover the full Navier-Stokes equations for isothermal fluid (Qian, 1992). The name comes f~om the adoption of the collision operator that had been employed in the kinetic BGK equation (Bhatnagar, 1954). The lattice BGK equation describes the propagation and collision of particles .occurring on the discrete spatial lattice at discrete time steps with a discrete velocity set. This equation appears to be Np,(x +cp,, t + l) - N.,~(x,0 = -( I/rXNp, - N ~ )
(2) where the collisionsof particlesare replaced by a relaxation process, in which the particle distribution is relaxed to itsequilibrium value over a time period, r. In order to reproduce visc~.usfluid flow which is governed by the Navier-Stokes equations, the equilibrium distribution should be Maxwellian and depend only on the local conserved densities of mass, momentum and energy. Unfortunately, the sound speed of the above model remains frozen so that it is difficult to increase the Mach number. It was claimed that the Maeh number must be less than x/5/2 I~hrough stability analysis on the ponitivity requirement of the particle density. More severe conditions would be expected if this analysis was done on the spatial gradient of hydrodynamic quantities. Therefore, some modifications of the model are needed to solve this difficulty. We have employed a particle reservoir to ensure mass conservation while the non-physical terms are removed by a modified Maxwell-Boltzmann type of equilibrium distribution. With such operations, we can freely set a certain number of density allocation factors on different sub-lattices. The state equation of the usual lattio: Boltzmann gas bears resemblance to that of perfect gas. Therefore, the phase transition phenomena could not occur in such a model. In regard to the state equation of gas in which the phase tramition is incorporated, the pressure is not always a monotonically increasing function of rmlss density. We have introduced into the LBE model a way to extend a recently suggested multi-speed
H. Ohashi et al./ Nuclear Engineering and L~sign 155 (1995) 67-71
LBE model to include phase transition at a prescribed temperature (Shan, 1993). We added a term in the original LBE with the consideration of inter-particle potential and analytically derived the critical magnitude of the inter-particle potential for the fixed critical density. The phase transition phenomenon is shown to occur on the multi-speed LBE model by carrying out a series of computer simulations.
69
C12
.........C21
3. Modeling tad calculation results Recently we extended a class of multi-speed LBE model into the compressible limit (Chen, 1993). The degree of freedom was introduced in the pressure term of the macroscopic momentum equation in the following way for a 2-D, threespeed model: c~p -- 0.[0.6( 1 -fo)P]
(3)
Here fo is the density allocating factor of the rest particles which is an arbitrary constant between 0 and 1. We realized that the free form of the pressure could actually be utilized to employ equations of state of real gas. In this model, the temperature has not been defined, so only the isothermal or isentropic cases could be considered. It is well known, however, that when the shock is weak it can be approximated as the isentropic process. This justifies the following simulation. Two-phase fluid can he simulated by considering the interaction potential between two fluid components in addition to the separate LBE description of the two components, This interaction potential becomes a connection term in the equation of state. We numerically simulated the propagation of a shock wave by using the aforementioned model. The lattice geometry and particle speeds we used are shown in Fig. 1. We performed the numerical simulation in the following way. The shock wave is formed by setting discontinuity of density in the initial conditions. The intensity of shock can be varied by adjusting the density ratio from 1.067 to 6.667. The Mach number and pressure ratios were measured when the steady stage was reached. The
L'14
c25
c2a
Fig. I. La'.ticegeometryand velocityvectors of nine-speed model. comparison between the results of our simulation and the theoretical prediction is shown in Fig. 2. We can sec that the deviation is negligible, in particular for cases of the small Mach number. Apart from the deviation analyzed above, we notice that the effect of high order terms, brought by the nature of the LBE method itself, becomes noticeable when the Mach numbs: becomes larger. These terms appear in the stre~g k,msor of the lattice gas and were derived recently for BGK models (Qian, 1993). Fig. 3 shows the pressure contour when an obstacle is put in the supersonic flow field. Compression shock can be seen in front of the obstacle with an expansion wave behind it. The inclusion of the inter-particle potential brings an additional term F~ into the RHS of the lattice Boltzmann equation. This additional collision term describes a part of the variation in particle distribution other than that by the local collision. As it is caused by the long-distance interaction of particles, we naturally relate it to the inter-particle potential V as follows F, -- - c , O ~ V .
(4)
Using this formulation we simudated the process of agglomeration. We employed a 100 x 100
70
H. Ohashi el al. I Nuclear Engineering anti Design 155 (1995) 67- 71
L~
4.5
I
o
t
.~ 4.0
~
O J
3.5
io" J
3.o o 2.5
¢~ 2.0
r'~
J
Y
J peril ct ga~ C, perf, c~, ga~ U
~-
1.0 1.0
1.1
1.2
1.3
1.4 1.5 1.6 Maeh Number
1.7
1.8
1.9
Fig. 2. LBE simulation of compressible shock.
(a)
t
= 100
(b)
t
= 300
(d)
t
= 2500
(£)
t
= 10000
Fig. 3. Pressure contour of supersonic flow with an obstacle.
lattice field in two space dimensions and set the critical density to be 0.693. The magnitude o f the inter-particle potential was chosen to be 0.27, well above the theoretically predicted critical value. In Fig. 4 the temporal evolution of the density field is shown, Gray levels in those figures are proportional to the local density; dark areas represent low-density regions while light color areas represent high-density regions. From this figure we clearly observe how the phase separation occurs and evolves. The scene of reconnection of interfaces of two bubbles can be recognized. This is a Iremendously difficult task to deal with for the conventional numerical techniques. It happened spontaneously in our simulation.
(a)
(e)
t
t
= 700
= 3500
Fig. 4. Temporal evolution of a tion process
density
field in the agglomera-
H, Ohashi et aLI Nuelea: Engineering and Design f55 (1995) 67-71
71
shock propagation and two-phase separation seem to be physically meaningful and it is coneluded that the method developed here becomes a powerful tool for the analysis of compt~ex fluid flow ~plmaring in fundamental processes of tim vapor explosion.
Appendix A: Nomenclatmre ~-~:,--': ~7-~/.W- .z
,.
~.. " ~:
~ ,,v..:2,/-:~
..
~.--., ,~2-
.....?,.,.
-..~.
.~
F interaction term between particles N particle distribution function V inter-particle potential particle velocity c density allocating factor fo pressure P A collision term P density relaxation time constant
--~.~ Ref~
Fig. 5. Interaction of plane shock with two-p~ase interfaces.
Lastly, employing the above simulation methods, we simulated the interaction between shock wave and two-phase interfaces. After the phase separation is built up, the plane shock is initiated at some place. As shown, for example, in Fig. 5, deformation and reflection of the shock wave could be clearly observed. 4. Conclusions We extended the range of applications of a multi-speed LBE model by introducing one degree of freedom in the sound speed of the fluid and inter-particle potential. The numerical results of
M. Ak~:ama (¢d.), Dynamics of the vapor explosion - - statt~ and perspective - - , Monbusho Gram-in-Aid Group Report. 1993. P. Bhatnagar, E.P. Gross and M.K. Krook, A re~xk'i for collision process~ in ga.7~s I, Phys. Rev. 94 (1954) 51t -51K Y. Chcal, H. Ohashi and M. Akiyama, The lattice Bottzm~rm equation method for the simulation of compre~bl¢ fluid flow, Proc. High Performara:¢ Computing in S/mulation MultiConference, March-April 1993, Arlington, 1993, pp. 27 - 32. B. Hasslacher, U. Frisch and Y. Pomcau, Lattice gas au:omata for the Navier-Stokcs equation, Phys. R ~ . Lett. 56 (1986) 1565-1567. G. McNamara and G. Zarmtti, Use of the Bo!tzmann equation to simulate lattice gas automata, Phys. Rev. Lctt. 61 { 1988) 2332-2334. Y.H. Qian, D. d'Humieres and P. Lallemaod, Lattice ~ K models for Navier-Stokes equations, Europhys. Left. 17 (1992) 479-484. Y.H. Qian and S.A. Orszag, Lattice BOK models for NavicrStokes equation; nonlinear deviation in compressible regimes, Europhys. 12tt. 21 (1993) 255-259. X. Shah and H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815-1819.
Nuclear Engineeringand Design 155 (1995) 67-71
Simulation of shock-interface interaction using a lattice Boltzmann model H. Ohashi, Y. Chert, M. A k i y a m a Department of Quantum Engineeringand Systems Science, Universityof Tokyo, Bunkyo-ku, Tokyo 113,Jat~
Abstract For simulation of fundamental processes of vapor explosion, a novel fluid dynamic simulation method is developed using the lattice Boltzmann model. We incorporate into the model capabilities to deal with comptess/b~ fluid flow and two-phase flow. By numerical simulation, we demonstrate fundanmntal c h a m t i c s of this medei and apply it to simulation of interaction betw~n shock wave and two-phase interfaces.
1. Introduction Simulation plays an important role in vapor explosion studies, in two possible ways: one is to reproduce the whole proczss of the vapor explosion to evaluate its effects, and the other is to simulate fundamental processes of the vapor explosion to understand their mechanisms (Akiyama, 1993). Formidable difficulties are encountered in the numerical work of the vapor explosion because it includes complex geometry, complex physics, phase changes and moving boundaries. In addition to that, conventional fluid descriptions require many constitutive models on which we hardly obtain information and validation from experimental results. So novel simulation methods based on mechanistic or physical backgrounds are need~a, m particular to elucidate the fundamental mechanisms by simulation. We have developed totally new algorithms to analyze fimdamental processes of the vapor explosion and applied them to simulation of pressure
shock, phase separation and interaction of pressure shock and two-pha.~ interface. These constitute key processes which caus¢ the coherent explosion. The lattice Boltzmann equation (LBE) mmthod is a recently developed numerical scheme to s i m ~ late fluid flow (McNamara, 1988). As a derivative of lattice gas automata (LGA) (Hasslacher, 1986), the LBE method solves macroscopk: c o m p ~ glow by dealing with the underlying micro-wotki. It describes the dynamics of partich) motion on a discrete lattice in terms of ensemble averaged particle distribution functions. Tim adv~mtages of adoption of the LBE model are two-fold in cornpar/son with the LGA model. First, the whit~ noise of the LGA model, which comes from the statistical nature of the automata, is completeiy removed. Second, the equations of ~ dynam/cs, the Navier-Stokes equations, are completely mcovered. Although the interua! stability ~ bib operation are no longer ensured in the LBE model, the simple implementation arid the ab~'ty
0029-5493/95/$09.50 © 1995 ElsevierScienceS.A. All rights reserved SSDI 0029-5493(94)00869-Z
68
H. Ohashi el at. / Nuclear Engineering and Design 155 0995) 67- 71
of accurate simulation, as well as the inheren.~ parallel nature, make the mt~el viable in the computational fluid dynamics (CFD) field. We have developed a class of multi-speed LBE model for compressible ~ow. We added one degree of freedom into the sound speed of modeled fltdd with the use of rest panicles and a particle reservoir. The sound speed is formalized to be dependent on thermodynamic variables in the way the perfect gas does in adiabatic cases. We also incorporated the capability of simulating multiphase immiscible fittids in the LBE model by introduction of multi-color particles and an interparticle potential. We will give results of shock wave simulation, phase separation and interaction of shock with two-phase interfaces.
2. LBE simdation The LBE simulation of fluid flow has become rapidly popular among the researchers of lattice hydrodynamics because it is noise-free and so flexible that the physics of real fluid can be fully recovered. Furthermore, the LBE models for fluid dynamics are extremely suitable for modern supereomputing because the calculation procedure is inherently parallel. This realization is due to some fundamental characteristics of the discrete microscopic model: the intrinsic stability, the easiness in deaiing with compficatexi geometric boundaries and the fast and local operations adequate for massively parallel processing. The lattice Boltzmann equation is expressed in the ensemble averaged form of the micro.dynamic equation for particles moving on a spatial lattice. It is given as follows Jv~,(x + c~, t + 1) -- N~,(x, t) = ~ ( N )
(1)
with the Boltzmann approximation applied to the collision term. Here Nm is the particle distribution and A~ is the collision term on the ith lattice node. Thus, the LBE exactly describes the twostep motion of particles in the micro-world. The first step is the streaming period during which the local particle distribution is transported to one of its nearc, t neighbors. The second step is the collision period during which the particle distribution
approaches its equilibrium value that depends locally on the macroscopic density and momentum. Recently a lattice BGK model was proposed to recover the full Navier-Stokes equations for isothermal fluid (Qian, 1992). The name comes f~om the adoption of the collision operator that had been employed in the kinetic BGK equation (Bhatnagar, 1954). The lattice BGK equation describes the propagation and collision of particles .occurring on the discrete spatial lattice at discrete time steps with a discrete velocity set. This equation appears to be Np,(x +cp,, t + l) - N.,~(x,0 = -( I/rXNp, - N ~ )
(2) where the collisionsof particlesare replaced by a relaxation process, in which the particle distribution is relaxed to itsequilibrium value over a time period, r. In order to reproduce visc~.usfluid flow which is governed by the Navier-Stokes equations, the equilibrium distribution should be Maxwellian and depend only on the local conserved densities of mass, momentum and energy. Unfortunately, the sound speed of the above model remains frozen so that it is difficult to increase the Mach number. It was claimed that the Maeh number must be less than x/5/2 I~hrough stability analysis on the ponitivity requirement of the particle density. More severe conditions would be expected if this analysis was done on the spatial gradient of hydrodynamic quantities. Therefore, some modifications of the model are needed to solve this difficulty. We have employed a particle reservoir to ensure mass conservation while the non-physical terms are removed by a modified Maxwell-Boltzmann type of equilibrium distribution. With such operations, we can freely set a certain number of density allocation factors on different sub-lattices. The state equation of the usual lattio: Boltzmann gas bears resemblance to that of perfect gas. Therefore, the phase transition phenomena could not occur in such a model. In regard to the state equation of gas in which the phase tramition is incorporated, the pressure is not always a monotonically increasing function of rmlss density. We have introduced into the LBE model a way to extend a recently suggested multi-speed
H. Ohashi et al./ Nuclear Engineering and L~sign 155 (1995) 67-71
LBE model to include phase transition at a prescribed temperature (Shan, 1993). We added a term in the original LBE with the consideration of inter-particle potential and analytically derived the critical magnitude of the inter-particle potential for the fixed critical density. The phase transition phenomenon is shown to occur on the multi-speed LBE model by carrying out a series of computer simulations.
69
C12
.........C21
3. Modeling tad calculation results Recently we extended a class of multi-speed LBE model into the compressible limit (Chen, 1993). The degree of freedom was introduced in the pressure term of the macroscopic momentum equation in the following way for a 2-D, threespeed model: c~p -- 0.[0.6( 1 -fo)P]
(3)
Here fo is the density allocating factor of the rest particles which is an arbitrary constant between 0 and 1. We realized that the free form of the pressure could actually be utilized to employ equations of state of real gas. In this model, the temperature has not been defined, so only the isothermal or isentropic cases could be considered. It is well known, however, that when the shock is weak it can be approximated as the isentropic process. This justifies the following simulation. Two-phase fluid can he simulated by considering the interaction potential between two fluid components in addition to the separate LBE description of the two components, This interaction potential becomes a connection term in the equation of state. We numerically simulated the propagation of a shock wave by using the aforementioned model. The lattice geometry and particle speeds we used are shown in Fig. 1. We performed the numerical simulation in the following way. The shock wave is formed by setting discontinuity of density in the initial conditions. The intensity of shock can be varied by adjusting the density ratio from 1.067 to 6.667. The Mach number and pressure ratios were measured when the steady stage was reached. The
L'14
c25
c2a
Fig. I. La'.ticegeometryand velocityvectors of nine-speed model. comparison between the results of our simulation and the theoretical prediction is shown in Fig. 2. We can sec that the deviation is negligible, in particular for cases of the small Mach number. Apart from the deviation analyzed above, we notice that the effect of high order terms, brought by the nature of the LBE method itself, becomes noticeable when the Mach numbs: becomes larger. These terms appear in the stre~g k,msor of the lattice gas and were derived recently for BGK models (Qian, 1993). Fig. 3 shows the pressure contour when an obstacle is put in the supersonic flow field. Compression shock can be seen in front of the obstacle with an expansion wave behind it. The inclusion of the inter-particle potential brings an additional term F~ into the RHS of the lattice Boltzmann equation. This additional collision term describes a part of the variation in particle distribution other than that by the local collision. As it is caused by the long-distance interaction of particles, we naturally relate it to the inter-particle potential V as follows F, -- - c , O ~ V .
(4)
Using this formulation we simudated the process of agglomeration. We employed a 100 x 100
70
H. Ohashi el al. I Nuclear Engineering anti Design 155 (1995) 67- 71
L~
4.5
I
o
t
.~ 4.0
~
O J
3.5
io" J
3.o o 2.5
¢~ 2.0
r'~
J
Y
J peril ct ga~ C, perf, c~, ga~ U
~-
1.0 1.0
1.1
1.2
1.3
1.4 1.5 1.6 Maeh Number
1.7
1.8
1.9
Fig. 2. LBE simulation of compressible shock.
(a)
t
= 100
(b)
t
= 300
(d)
t
= 2500
(£)
t
= 10000
Fig. 3. Pressure contour of supersonic flow with an obstacle.
lattice field in two space dimensions and set the critical density to be 0.693. The magnitude o f the inter-particle potential was chosen to be 0.27, well above the theoretically predicted critical value. In Fig. 4 the temporal evolution of the density field is shown, Gray levels in those figures are proportional to the local density; dark areas represent low-density regions while light color areas represent high-density regions. From this figure we clearly observe how the phase separation occurs and evolves. The scene of reconnection of interfaces of two bubbles can be recognized. This is a Iremendously difficult task to deal with for the conventional numerical techniques. It happened spontaneously in our simulation.
(a)
(e)
t
t
= 700
= 3500
Fig. 4. Temporal evolution of a tion process
density
field in the agglomera-
H, Ohashi et aLI Nuelea: Engineering and Design f55 (1995) 67-71
71
shock propagation and two-phase separation seem to be physically meaningful and it is coneluded that the method developed here becomes a powerful tool for the analysis of compt~ex fluid flow ~plmaring in fundamental processes of tim vapor explosion.
Appendix A: Nomenclatmre ~-~:,--': ~7-~/.W- .z
,.
~.. " ~:
~ ,,v..:2,/-:~
..
~.--., ,~2-
.....?,.,.
-..~.
.~
F interaction term between particles N particle distribution function V inter-particle potential particle velocity c density allocating factor fo pressure P A collision term P density relaxation time constant
--~.~ Ref~
Fig. 5. Interaction of plane shock with two-p~ase interfaces.
Lastly, employing the above simulation methods, we simulated the interaction between shock wave and two-phase interfaces. After the phase separation is built up, the plane shock is initiated at some place. As shown, for example, in Fig. 5, deformation and reflection of the shock wave could be clearly observed. 4. Conclusions We extended the range of applications of a multi-speed LBE model by introducing one degree of freedom in the sound speed of the fluid and inter-particle potential. The numerical results of
M. Ak~:ama (¢d.), Dynamics of the vapor explosion - - statt~ and perspective - - , Monbusho Gram-in-Aid Group Report. 1993. P. Bhatnagar, E.P. Gross and M.K. Krook, A re~xk'i for collision process~ in ga.7~s I, Phys. Rev. 94 (1954) 51t -51K Y. Chcal, H. Ohashi and M. Akiyama, The lattice Bottzm~rm equation method for the simulation of compre~bl¢ fluid flow, Proc. High Performara:¢ Computing in S/mulation MultiConference, March-April 1993, Arlington, 1993, pp. 27 - 32. B. Hasslacher, U. Frisch and Y. Pomcau, Lattice gas au:omata for the Navier-Stokcs equation, Phys. R ~ . Lett. 56 (1986) 1565-1567. G. McNamara and G. Zarmtti, Use of the Bo!tzmann equation to simulate lattice gas automata, Phys. Rev. Lctt. 61 { 1988) 2332-2334. Y.H. Qian, D. d'Humieres and P. Lallemaod, Lattice ~ K models for Navier-Stokes equations, Europhys. Left. 17 (1992) 479-484. Y.H. Qian and S.A. Orszag, Lattice BOK models for NavicrStokes equation; nonlinear deviation in compressible regimes, Europhys. 12tt. 21 (1993) 255-259. X. Shah and H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815-1819.