Numerical Investigation on Physical Foaming and Decay Process Using Multicomponent Thermal Lattice Boltzmann Model

Foam materials produced by physical process are common in the field of material engineering such as foamed asphalt. Considering heat and phase transfer, the complicated nonlinear phenomenon of foaming and decay process is difficult to be described numerically. In this paper, a multicomponent thermal lattice Boltzmann model for simulation of physical foaming and decay process is proposed. The model combines a thermal single-component phase change model with a multicomponent model, and the two models are verified separately. The physical foaming and decay process is solved numerically based on the proposed model. Simulation results show the influence of the distribution and content of the phase-change component on the expansion rate and decay time during physical foaming and decay process.


Introduction
Physical foaming and decay process caused by bubbles generation and crushing of one component (the phase-change component) in multicomponent fluids, such as foamed asphalt foaming process [1] is widespread in material engineering and other field. The macro dynamic behaviour of such foamflow is very complex because of the interaction between different fluids and the process of mass and heat transfer. To investigate the foam-flow, multicomponent multiphase fluids accompanied by phase change with heat transfer need to be described numerically. However, conventional computational fluid dynamics (CFD) methods face difficulties in predicting such complex flow.
Different from conventional CFD methods, based on microscopic nature and mesoscopic characteristics, the lattice Boltzmann method (LBM) performs well in numerical simulations of many kinds of complex flow. Many thermal lattice Boltzmann models have been developed in the last few years [2][3][4]. However, those models were only used to study the phase change in single component multiphase fluids. In this work, we focus on thermal multicomponent multiphase fluids, of which one component is accompanied by the phase change process. The single-component phase change model presented by Gong and Cheng [5] is combined with a multicomponent model to simulate the phase change process of one component in multicomponent fluids with heat transfer, which can be used to investigate the foaming and decay process of foam-flow. The rest of the paper is organized as follows. Section 2 reviews the presented thermal phase change model and multicomponent model, and describes the combination of the two models. In Section 3, the phase change model and the multicomponent model are validated separately. In Section 4, foaming and decay process is simulated and the influence of the distribution and content of the phase-change component on the two important indicators (expansion rate and decay time) for studying foam-flow is discussed. Conclusions are made in Section 5.

Thermal single-component phase change model
In LBM, the fluid is described by the evolution of density distribution functions. The evolution equation of the density distribution function with the Bhatnagar-Gross-Krook collision operator [6] is written as: is the density distribution function along α direction at the lattice site r and time t, δ t is the time step, e а is the discrete velocities, τ is the dimensionless relaxation time, r is the body force term, N is the number of the discrete velocities(In this work, we use two dimensions with nine directions (D2Q9) lattice configuration), ( , ) is its corresponding equilibrium distribution function, which has the following form: where c s is the lattice sound speed which is equal to 1/ 3 , w D is the weight coefficient, and 0 1 4 5 8 4 / 9, 1 / 9, 1 / 36 w w w . ρ and u are the density and the macroscopic velocity, respectively. They are defined as: The exact-difference-method force scheme [7] is adopted and ( , ) f t ' r is given by: is the velocity change caused by body force during time step δ t , F is the body force including the interparticle interaction force F int , the gravitational force F g and so on. F int and F g can be obtained as follows: where c 0 = 6.0 [8] and g is a Green's function, which depends on the strength of the interaction force between particles, and represents attractive (repulsive) force when g < 0 (g > 0). G is the acceleration of gravity and ave U is the average density of the whole computation domain. ( ) \ r is the effective mass at lattice site r, and it is a function of ρ , which can be given by: where P(r) is the local pressure. Peng-Robinson (P-R) equation of state (EOS) is selected to calculate local pressure and introduce the influence of temperature field on density distribution. It is given by:
By averaging the moment before and after the collision, the actual physical velocity U is given by To solve energy equation, another distribution function is introduced. The evolution equation of the temperature distribution function is written as: ( where T W is the relaxation time for temperature with thermal diffusivity, λ is the thermal conductivity and c υ is the specific heat capacity. ( , ) eq h t D r is the equilibrium temperature distribution function which has the following form: T is the macroscopic temperature which is obtained by T h D D ¦ , φ in Equation (9) is the source term describing the energy change due to phase change. Gong and Cheng [5] derived the source term from the thermodynamic relation of entropy. And it can be written as:

Multicomponent model
The single-component pseudopotential model can be easily extended to a multicomponent model by treating each component using a separate density distribution function. Each set of density distribution functions can be evolved separately following the Equation (1). In order to separate different components, the interparticle force between different components is introduced into the multicomponent model. Summing over all the components and interacting sites, the interparticle force between different components acting on the particles of the σth component at site r, F intmc , has the following form like Equation (5): And the macroscopic density of multicomponent fluids is defined as: The macroscopic velocity and the actual physical velocity can be obtained as follows: where V F is the total body force of the σth component.

Combined model
To simulate foaming and decay process, the two models mentioned above are combined. For simplicity, we consider fluids consisting of two components. We set water as the first component and the effective mass 1 \ is calculated by Equation (6) with P-R EOS. As the temperature rises, a liquidvapour phase change process of the first component will appear. The second component is set to be a kind of component without phase change in our simulation, and the effective mass 2 \ is given by 2 2 \ U . The EOS of the multicomponent fluids [9] is calculated by:  (11), the source term in the evolution equation of the temperature can be written as: 1 1 212 It should also be mentioned that density of the first component is different compared with which that in the single-component thermal phase change model at the same temperature. That is caused by the interparticle force between different components. In this work, we set ρ 2 = 1.0, g 11 = 0.1(or -0.1), g 12 = 0.005, and g 22 = 0. For the single-component flow, the specific heat capacity (c υ ) is considered constant. However, in a multicomponent system, c υ of the mixture flow is no longer a constant. It is related to the distribution of each component. We define it as a parameter, which changes along with the mass fraction of each component. c υ can be expressed as:

Model verification
As the combined model is proposed by integrating the thermal phase change model and the multicomponent model, we verify two models separately. The phase change model is validated through the simulation of bubble growth and departure from a horizontal surface. Meanwhile we verify the multicomponent model for surface tension by binary static droplet simulations. It should be noted that variables in this paper are in lattice units, which can be converted into real physical properties based on the principle of reduced properties.

Thermal phase change model validation
Bubble growth and departure from a heated horizontal isothermal surface is simulated to study the relationship between departure diameter of a bubble and gravity acceleration. And the results are compared with analytical solutions and the simulation results in Gong and Cheng's [5] work. A 150 × 450 lattice structure is chosen to simulate the bubble growth and departure from a heated horizontal isothermal surface. Initially, the computation domain is occupied by the liquid phase. At the centre of the bottom wall, a high temperature core is specified to ensure that nucleation starts there. And simulation conditions are set as Gong and Cheng's [5] work. As time goes by, a bubble will grow and depart from the bottom wall. The process of bubble growth and departure under the condition of G = 0.000 04 is shown in Figure 1.   The relationship between departure diameter of the bubble and gravity acceleration [10] is given as follows: the simulation results of bubble departure diameter D are similar with the reference results in Gong and Cheng's [5] work, and can be fitted into regression formula D = 0.3061|G| -0.4902 . The exponent of G is -0.4902, which agrees well with the analytical result given by Equation (18). And simulation results showed above verify the correctness of the thermal phase change model.  Expansion rate and decay time are two important indicators for studying foam-flow. The liquid level of the second component, which can be used to investigate the expansion rate and decay time of foaming and decay process, is measured during bubble generation and crushing. The influence of distribution and content of the first component on the expansion rate and decay time is investigated using the control variable method. Figure 5(a) shows the average liquid level of the second component during foaming and decay process under different first component distributions (the initial number of the droplets is 2×3, 2×4, 2×5 separately at the same content of the first component). We can see that droplets distribution does not affect the expansion rate, but the more uniform the distribution, the slower the decay is. And Figure 5(b) shows the average liquid level of the second component during foaming and decay process under different first component contents (the initial radius of droplets is 4, 5, 7 separately at the same number of the droplets). We can see that as the content of the first component increases, the expansion rate increases and decay becomes faster. And step decay occurs happens when bubbles collapse as shown in Figure 5. It is because the bubbles generated in this simulation are large and the number of the bubbles is small, just as Jenkins mentioned [11].

Conclusion
In this paper, we propose a multicomponent lattice Boltzmann model for simulation of physical foaming and decay process caused by bubbles generation and crushing of one component (the phasechange component) in multicomponent fluids. The model combines a single-component thermal phase change model and a multicomponent model. Both models are verified. The combined model is used to simulate physical foaming and decay process. It is found that the distribution of the phase-change component does not affect the expansion rate, but the more uniform the distribution, the slower the decay is. And simulation results also show that as the content of the phase-change component increases, the expansion rate increases and decay becomes faster.