Numerical investigation of two-phase flow through tube bundles based on the lattice Boltzmann method

ABSTRACT In pursuit of stable operation in thermal hydraulic engineering, this paper researches two-phase flow through tube bundles by the lattice Boltzmann method. The single-relaxation-time lattice Boltzmann model and the multicomponent, multiphase pseudo-potential lattice Boltzmann model are adopted. Specifically, the two-phase flow past double-tandem circular cylinders, double-parallel circular cylinders, in-line tube bundles and staggered tube bundles with various spacing ratios is investigated. The vortex variable distribution and vortex shedding law under different working conditions are further compared to obtain the time-averaged drag, oscillating drag and Strouhal number. The variation law of lift coefficient and lift power spectrum with different spacing ratios is revealed. With increasing spacing ratio, the mutual interference between the double-tandem circular cylinders is weakened, the amplitude of the time-averaged drag on the double-parallel circular cylinders is reduced and the lift power spectrum changes from a double-peak distribution to a single-peak structure. The drag and lift of the center column of in-line tube bundles are lower than those on both sides, while the oscillating drag of the odd-row columns of staggered tube bundles is higher than that of even-row columns. These results contribute toward laying a solid foundation for future research on multiphase tube bundle flow issues in nuclear engineering.


Introduction
As a clean and low-carbon energy source, nuclear energy can reduce greenhouse gas emissions and mitigate the environmental effects of energy consumption while meeting the growing demand for electricity. The active and steady development of nuclear power is an important technological choice for China to achieve its carbon peak and carbon neutrality goals. As several major accidents in the nuclear industry have had a serious impact on social development and environmental protection, scholars have begun to search for reactor types with inherent safety. The high-temperature gas-cooled reactor, one of the fourth-generation advanced reactor types, is expected to have promising prospects for development Wu et al., 2022).
The two-phase flow issue is a ubiquitous and complicated phenomenon with a wide range of applications in industry and engineering, such as applications in nuclear reactors (Amini & Schleiss, 2009;Li et al., 2020;Mosavi et al., 2019). For instance, there is gas-liquid two-phase flow between fuel rod bundles in a boiling water reactor, and the scouring of the bundles by the fluid generates CONTACT Nan Gui guinan@mail.tsinghua.edu.cn; Shengyao Jiang shengyaojiang@sina.com alternating shedding vortices. At the same time, periodic Karman vortex streets are formed in the wake flow, triggering oscillations in the flow and heat transfer characteristics (Poullikkas, 2003). The long-term oscillations around the columns can lead to material fatigue, which may cause deformation and wear of the pipeline. It can also jeopardize the safety and reliability of nuclear reactor equipment if resonance phenomena occur (Q. Lu et al., 2019). Furthermore, research on the bubble generation mechanism and prediction method (Ban et al., 2015) has aroused scholars' interest, but there is no uniform perception. Therefore, an in-depth study on the thermodynamic properties (Lahey & Drew, 2001) of the two-phase flow process is needed to reveal the intrinsic mechanism of the engineering phenomenon and to ensure operational safety. Two-phase flow past circular cylinders or tube bundles is usually found in industrial equipment with heat exchange, such as steam generators and condensers. Extensive numerical and experimental studies have been conducted on the thermal hydraulics. Some scholars have studied the distribution of the void fraction, fluid velocity and hydrostatic pressure around the cylinders with vertical upward bubble flow (Inoue et al., 1986), and analyzed the variation of vortex, pressure and velocity fields in the wake flow, as well as the impact on the distribution of the void fraction. Studies have also been conducted on twophase flow past columns with different shapes (J. Lu et al., 2000), where the pressure distribution, oscillating lift and Strouhal number are evaluated under the variation of the void fraction and Reynolds number.
Scholars have also carried out in-depth experimental investigations to explore the flow past columns with different arrangements and spacing ratios. Some summarized the vortex shedding forms of single-phase flow with double-tandem cylinders and double-parallel cylinders with different spacing ratios (Zdravkovich, 1985). Scholars also explored the vortex shedding forms behind in-line or staggered tube bundles with various transverse and longitudinal spacing ratios (Ishigai, 1973). Vortex streets may appear behind a single column or be generated behind two columns, and a shedding boundary layer from the first column may be attached to the second column. The variation of drag coefficient and void fraction in a single column or tube bundles in two-phase flow has also been studied (Hazuku et al., 2020;Joo & Dhir, 1994). Owing to the complexity of gas-liquid interaction in two-phase flow, involving the vortex induced by bubbles coiling into the column (Breuer et al., 2000), the effect of bubbles on vortex shedding and the boundary layer (Sun et al., 2017), the numerical study of two-phase flow is difficult.
In recent years, the multiphase lattice Boltzmann method has been shown to have outstanding advantages in solving two-phase flow problems. The method selects appropriate model parameters and sets a certain density perturbation in the initial density field. By adopting the Shan-Chen force, an automatic separation between the two phases can be achieved (Y. Wu et al., 2018), then the flow field with an initialized uniform distribution can be obtained, thus enabling a simple and efficient simulation of the two-phase flow problem. Some researchers used the multiphase lattice Boltzmann method, and the results of validation cases were in good agreement with previous numerical results or theoretical solutions. The simulation on vertical and horizontal multiphase pipe flow was carried out (Cheng et al., 2018). On this basis, scholars continued to study the flow past square or cylindrical columns based on the multiphase lattice Boltzmann method. By analyzing the interaction between continuous and discrete phases, as well as the lift and drag forces acting on the solid, the intrinsic mechanism of the two-phase flow was revealed (Cheng et al., 2020).
Vortices play a key role in the turbulence mechanism, but there are still many issues that are not yet well resolved; in particular, there is no universally accepted definition of vortices. The first generation of vortex identification methods selected vorticity to identify and display vortices. However, subsequent studies have shown that there is a significant gap between vorticity and vortices in terms of size and direction (Wang et al., 2017). The correlation between vorticity and the vortex core is relatively low, because in some cases the vorticity at the location of the vortex core is small while the vorticity outside the core is large. In addition, the vector direction of the vorticity is often not aligned with the direction of the vortex core. Second-generation vortex identification methods based on eigenvalue-based criteria have been proposed, including the Q method (Hunt et al., 1973), method (Chong et al., 1990), λ 2 method (Jeong & Hussain, 1973), λ ci method (J. Zhou et al., 1999) and method (C. Liu et al., 2016). However, since these vortex identification variables are scalars, they cannot identify the direction of the rotation axis. Consequently, part of the information on the vortex vector is lost. Since the identified vortex structures depend on the selection of a threshold value, different threshold values will produce different vortex structures. which cannot accurately describe the rotational intensity of the fluid. So, the second-generation vortex identification methods cannot quantitatively study the vortex structure. Therefore, the above methods have certain defects in describing vortex structure and vortex direction.
Liu and co-workers  proposed and developed the third-generation vortex identification method, which provided a more accurate approach for vortex identification. The method is to decompose the vorticity into a rigid rotating part and a non-rotating part, so that it can identify the local vortex axis and is less sensitive to the threshold value. Based on the third-generation vortex identification method, scholars have extensively investigated the vortex structures, such as hairpin vortices (Gao et al., 2019), open channel turbulent vortices (Alam et al., 2022;Bai et al., 2019) and rotating jet vortices . Gui, Qi, et al. (2019) summarized the relationship between acceleration characteristics and the vortex core in swirl flow. The results showed that under different swirl number conditions, the instantaneous acceleration at different times and locations and the Liutex vectors are always perpendicular or parallel to each other, verifying the accuracy of the Liutex method for identifying the vortex core.
This paper investigates the two-phase flow through tube bundles based on the single-relaxation-time lattice Boltzmann model, and the multicomponent and multiphase pseudo-potential Boltzmann model. The two-phase flows past double-tandem circular cylinders, double-parallel circular cylinders, in-line tube bundles and staggered tube bundles are analyzed with various spacing ratios. The vortex variable distribution and vortex shedding patterns under different working conditions are further compared. The time-averaged drag, oscillating drag and Strouhal number are obtained and analyzed. The variation law of lift coefficient and lift power spectrum with different spacing ratios is revealed. The research results of this paper are likely to lay a solid research foundation for future two-phase flow issues in the nuclear engineering field.

Single-relaxation-time lattice Boltzmann model
Similarly to the traditional computational fluid dynamics methods, the lattice Boltzmann method discretizes in time and space. The special feature of this method is that it also discretizes the velocity space. Since the details of the microscopic particle motion prove to have an insignificant effect on the macroscopic motion of the fluid, the particle velocity can be reduced to a finite dimensional velocity space. The lattice Boltzmann method originates from the theory of gas dynamics, where Boltzmann proposed the governing equation describing the spatiotemporal evolution law of the nonequilibrium velocity distribution function, namely the Boltzmann equation: where f (x, v, t) denotes the velocity distribution function of the particles, representing the density of the particle number with a velocity along the v = (v x , v y , v z ) direction in the spatial position x = (x, y, z) at moment t.
(f ) represents the collision operator, which portrays the interaction between microscopic fluid particles. It plays a key role in the ability of the model to accurately capture the intrinsic physical laws of the fluid system.
Owing to the complex nonlinear integration involved in its computation, a linear approximation is generally employed through the single-relaxation-time model proposed by Bhatnagar et al. (1954). The essence of the method is to perform a linear approximation to the collision operator, considering that the effect of the collision converges the velocity distribution function to the Maxwell equilibrium distribution function f eq . Assuming that the change rate is proportional to the difference f eq − f , a scaling factor υ is introduced: where υ represents the total number of particle collisions with other particles per unit time, namely the collision frequency, and τ c characterizes the average time interval between two particle collisions, also known as relaxation time, which is related to the physical properties of the fluid. The Boltzmann-BGK equation can be obtained by combining the two formulations above. Since only one relaxation parameter is involved, the corresponding model is the single-relaxation-time lattice Boltzmann model. Owing to the simplicity of the model and ease of computation, the single-relaxation-time model has been widely employed. By discretizing in spatiotemporal and velocity space, the discretized single-relaxation-time lattice Boltzmann equation can be obtained: where c i represents the discrete particle velocity, f i (x, t) represents the particle density distribution function along the c i direction, and τ is the dimensionless relaxation time.
The discrete equilibrium distribution function can be calculated by the following equation: where c s is the lattice sound velocity, and ω i is the weighting factor corresponding to c i . In the implementation of the computational procedure, the solution of the lattice Boltzmann equation can be generally decomposed into two parts, the collision step and the steaming step, respectively: where f * i (x, t) denotes the distribution function after collision. The collision step may contain nonlinear calculation, but only involves the computation within the local lattice points, which has the characteristic of local nonlinearity. The streaming step distribution function propagates along the discrete velocity direction to the neighboring lattice points, but the calculation process is linear and has the characteristic of non-local linearity.
Through the Chapman-Enskog multiscale expansion method, the macroscopic equations corresponding to the single-relaxation-time lattice Boltzmann equation can be obtained: where kinematic viscosity is related to relaxation time, and pressure is related to density. Density, velocity and internal energy can be calculated by moments of the distribution function; that is, the integration for a continuous distribution function or summation for a discrete distribution function:

Multiphase pseudo-potential lattice Boltzmann model
In the multiphase flow simulation, the Shan-Chen pseudo-potential model (Shan & Chen, 1993) is widely employed in the research, as it can automatically track the phase interface in the calculation, which is more convenient and efficient for simulating a multiphase interface system. Consequently, the multicomponent and multiphase pseudo-potential model is chosen for calculation in this paper. The pseudo-potential model introduces the Shan-Chen force for the interaction between fluid particles, and the combined force on the particles located at x can be calculated by integrating the force on the particles at each surrounding location (Wu et al., 2021): where ψ(ρ) is the pseudo-potential function, and G(x, x ) is the Green's function that determines the strength and range of the interaction, the simplified form of which generally considers only the interaction between near-neighboring lattice points. For a multicomponent and multiphase system consisting of n components, the single-relaxation-time lattice Boltzmann evolution equation for each component is as follows: where , τ σ and F σ i denote the distribution function, the equilibrium state distribution function, the relaxation time and the external force applied to the σ component, respectively.
The interaction force among particles applied to σ components can be accumulated by the other components: To achieve the characteristic of different wetting properties of the walls, the interaction force between fluid and solid is introduced: (16) where s(x) marks the properties of the lattice point, the fluid lattice point is set as 0 and the solid lattice point is set as 1. G σ S characterizes the strength of the interaction between the fluid particles and the solid wall. To realize the hydrophilic or hydrophobic walls, the value of G σ S is set negative or positive. In this paper, G σ S is set to zero to achieve a non-wetting or neutral surface. The calculation model in this paper is a twocomponent multiphase system without considering phase change, so the value of n is equal to 2 and the value of G σ σ is positive, so as to achieve repulsion and separation among different components. The total density of the fluid, the pressure, the effective velocity of motion and the actual velocity of the fluid can be calculated as follows: The calculation in the lattice Boltzmann method uses lattice units, which are lattice length lu and lattice time lt. Corresponding to the actual physical process, the conversion of lattice units to actual units needs to be completed. The conversion principle is based on the similarity principle, so the dimensionless numbers need to be consistent. In multiphase and multicomponent systems, dimensionless parameters such as the Weber number, Reynolds number, viscosity ratio and density ratio produce major effects. Therefore, by keeping the relevant variables equal in lattice units and actual units, the correspondence between the simulation and the actual problems can be obtained. In this paper, lattice units are adopted for all numerical cases, including lattice length lu, lattice time lt and lattice velocity lu/lt.

Vortex identification method
The third-generation vortex identification method provides a more accurate means of vortex identification, which decomposes the vorticity ω into two parts, S and R. R represents the rigid rotating part Liutex R, and S represents the non-rotating part. This decomposition method provides the only tensor decomposition with Galilean invariance. The main component of ω in most cases is S, as the magnitude of R is smaller. The traditional method of characterizing fluid rotation by vorticity causes larger errors. The direction of the Liutex r represents the direction of the local rotation axis, which is determined by the real eigenvector in the velocity gradient tensor. The magnitude of Liutex represents the intensity of the local rotation.

Numerical method validation
In this section, the single-phase and the multiphase lattice Boltzmann model are thoroughly validated. The cases of single-phase flow past a single cylinder or two tandem cylinders are investigated.
The single-phase model is validated by the threedimensional lid-driven cavity flow under different Reynolds number values. The simulation of a cubic cavity was conducted with the Reynolds number at 100, 400 and 1000. The simulation results are shown in Figure 1(a)-(f), which are in good agreement with reference data (Ku et al., 1987).
The multiphase model is validated by the Laplace law and multicomponent Poiseuille flow. In the multicomponent Laplace test, a circular droplet of one phase with radius R is initially surrounded by another liquid phase.
As the pressure difference between the inside and outside of a stationary droplet is proportional to the surface tension, and inversely proportional to the droplet radius, the Laplace law takes the form of Δp = σ /R. The results of the Laplace test are illustrated in Figure 1(g) and (h), where the relationship between the measured pressure difference and radius is linear, with the coefficient of determination of the fitted lines above 0.999. In the multicomponent Poiseuille flow test (Banari et al., 2014), the center of the flow channel is occupied by liquid a with viscosity υ a and density ρ a , and the sides of flow channel are occupied by liquid b with viscosity υ b and density ρ b . The left and right sides of the flow channel are solid wall boundaries. As shown in Figure 1(i) and (j), the calculated velocity distribution fits well with the theoretical equation, and there is only a certain velocity error at the gas-liquid interface.
In addition, the correlation coefficients are verified for the single-phase flow past the two-dimensional circular cylinder. First, the laminar flow around a cylinder in the channel is simulated. As shown in Table 1, the accuracy of the numerical method is verified by comparing the calculation with previous simulation results (Schäfer et al., 1996), which are in good agreement with the reference values. Secondly, the laminar flow past two tandem cylinders with different space ratios is investigated. The Reynolds number is 200 and five sets of spacing ratio S/D are chosen, i.e. 1.5, 2.0, 3.0, 4.0 and 5.0. In Table 2, the simulation results are compared with reference data (S. Liu & Fu, 2000;K. Zhou et al., 2018), where C D1 and C D2 indicate the drag coefficient of the upstream and downstream cylinders, respectively. Concerning cases with a small spacing ratio, the drag coefficient of the downstream cylinder is negative, induced by the shielding effect from the upstream cylinder, in accordance with the previous experimental and simulation results. As the spacing ratio increases, the drag and lift coefficient of the downstream cylinder increases significantly, especially as the drag coefficient changes from negative to positive.

Numerical setting
Two-phase flow through tube bundles has a wide application in industrial equipment with heat exchange. The channel flow through nuclear fuel rod bundles in boilingwater reactors is a typical two-phase flow scenario. In this paper, the two-phase flow through double-tandem circular cylinders, double-parallel circular cylinders, in-line tube bundles and staggered tube bundles with different spacing ratios has been studied, which is valuable for engineering applications.
The solid wall surface is set to be non-wetting, and periodic boundaries are applied in the flow direction. The gravity in the vertical direction is applied to drive the fluid motion.
For the two-phase flow past double-tandem circular cylinders, six sets of spacing ratio S/D are chosen: 1.2, 1.6, 2.0, 3.0, 4.0 and 5.0. The computational domain (L x , L y ) is (256, 2048)lu, and the diameter of the cylinder D is 100lu. Applying the same vertical gravity value, the Reynolds number of the computational condition is 120.   For the two-phase flow past double-parallel circular cylinders, five sets of spacing ratio S/D are chosen: 1.2, 1.6, 2.0, 2.5 and 3.0. The computational domain (L x , L y ) is (640, 1536)lu, and the diameter of the cylinder D is 100lu.
Applying the same vertical gravity value, the Reynolds number of the computational condition is approximately 120.
For the case of two-phase flow through in-line tube bundles, three sets of longitudinal spacing ratios are chosen: 2.0, 4.0 and 8.0. For the case of two-phase flow through staggered tube bundles, two sets of longitudinal spacing ratio are chosen: 4.0 and 12.0. The transverse spacing ratio is 1.5 for the above cases, the computational domain (L x , L y ) is (640, 3072)lu, and the diameter of the cylinder D is 128lu. The Reynolds number for the cases is approximately 100, and grid independence tests have been performed for the above numerical calculation.  Figures 2 and 3 show the evolution of the vortex variables R and ω in the two-phase flow past double-tandem circular cylinders. The flow field velocity reaches equilibrium under a Reynolds number of 120. At a small spacing ratio (S/D = 1.2), the fluid flow between the two cylinders is slow, corresponding to a slender column, resulting a single vortex street in the wake of the downstream cylinder. As the spacing ratio increases (S/D = 1.6), the boundary layer separated from the upstream cylinder adheres to the rear column in a quasi-static form. When the spacing ratio continues to increase (S/D = 2 and 3), vortices are already formed behind the upstream cylinder, and detached vortices occasionally adhere to the rear cylinder or move downstream. When the spacing ratio is relatively large (S/D = 4 and 5), the shear layer separated by the upstream cylinder is not attached to the downstream column, and there are vortex formation and alternate detachment after both upstream and downstream columns, characterized by the form of double vortex streets. The R and ω have positive and negative discrete distribution areas on both sides of the column, which spread to the rear sides with the motion of bubbles. It can also be found that with the shedding of the vortex, its interference with the flow field will also affect the gas-liquid two-phase distribution pattern of the downstream columns. The disturbance of the flow field further affects the flow field distribution of the downstream column. With the increase in the spacing ratio, the gas-liquid two-phase separation is more obvious owing to the full development of the vortex. Figure 4 shows the evolution of the drag and lift forces, and the analysis of the lift force spectrum of the doubletandem cylinders. The drag and lift forces of the upstream cylinder are represented by FU D and FU L, which are shown as solid blue lines in the figure. FD D and FD L denote the drag and lift forces of the downstream column, which are indicated as red dashed lines in the figure.

Two-phase flow past double-tandem circular cylinders
The oscillating lift force in the two-phase flow process consists of three parts. The first one is the oscillating lift force triggered by the periodic alternating shedding vortex in the wake to form a vortex street, which has a large amplitude and concentrates the main energy of vibration. Secondly, as the two-phase flow has a large turbulence intensity, its random pulsation triggers the oscillating lift. The Reynolds number of the working conditions in this paper is low, so the resulting oscillating lift force is smaller. The third part is the oscillating lift triggered by the bubbles with different density, which continuously hit the cylinders. The first part of the oscillating lift triggered by vortex shedding has obvious periodicity, while the last two parts have a certain degree of randomness.
As can be seen in Figure 4, when the spacing ratio ranges from 1.2 to 2.0, the drag and lift forces change steadily. With the bubble impacting the cylinder to produce a certain pulsation, the drag of the upstream cylinder is obviously larger than that of the downstream cylinder. When the spacing ratio continues to increase, the drag forces of the upstream and downstream cylinders are comparable in magnitude. As for the lift force, it basically remains comparable. When the spacing ratio is in the range of 3.0-5.0, the magnitude of the lift force on the downstream cylinder is larger. Fast Fourier transform spectrum analysis of the lift force indicates that the oscillating frequencies of the upstream and downstream cylinders are basically the same at different spacing ratios. Consequently, the corresponding vortex street frequency is generally equal. There exists a difference between the power spectrum peak values of the upstream and downstream cylinders. When the spacing ratio ranges from 1.2 to 2.0, the main frequency peak value of the upstream cylinder is higher than that of the downstream cylinder, which is higher by 58.80%, 24.87% and 10.55%, respectively. When the spacing ratio is in the range of 3.0-5.0, the main frequency peak value of the downstream cylinder is higher than that of the upstream cylinder, which is higher by 17.13%, 27.12% and 29.95%, respectively. Table 3 presents the statistical characteristic values of the drag and lift for each cylinder at different spacing ratios, where U1.2 and D1.2 denote the upstream and downstream cylinders with S/D of 1.2, respectively, the other groups of marks have similar meanings, and N1 denotes the case of one single cylinder. The time-averaged drag coefficient C D , oscillating drag coefficient C D, time-averaged lift coefficient C L , and oscillating lift coefficient C L are defined as follows: Figures 5 and 6 present the comparison of the drag and lift for each cylinder at different spacing ratios. In Figure 5(a), the value of C D for a single cylinder is 2.26, while the C D of any column in the double-tandem circular cylinders is smaller than that owing to the shielding effect. When the spacing ratio is between 1.2 and 2.0, the downstream cylinder is located in the wake of the upstream cylinder. Consequently, the drag coefficient is significantly lower, and the C D of the upstream cylinder is higher than that of the downstream cylinder. With the increase in the spacing ratio, the C D of the upstream cylinder changes less, with a maximum increase of 4.98%. As the C D of the downstream cylinder gradually increases, the gap between them gradually decreases. Since the boundary layer separated from the upstream cylinder is no longer attached to the downstream cylinder, the C D of the upstream and downstream cylinders is nearly equal for S/D = 4. When S/D = 5, the C D of the downstream cylinder is higher than that of the upstream cylinder by 7.28%. In Figure 5(b), when the spacing ratio is between 1.2 and 2.0, the oscillating drag coefficient C D of the upstream cylinder is significantly higher than that of the downstream cylinder. As the spacing ratio increases, the difference between them gradually decreases. For S/D = 4 or 5, the C D of the downstream cylinder is higher than that of the upstream cylinder.
In Figure 6(a), the time-averaged lift coefficients C L of the upstream and downstream cylinders show that the      averaged lift forces are in the same or opposite directions, and the magnitude of C L is larger when the spacing ratio is between 3.0 and 5.0. In Figure 6(b), the variation pattern of the oscillating lift coefficient C L of the upstream and downstream cylinders is similar to C D in Figure 5(b). As the spacing ratio increases, the difference between them gradually decreases. With a larger spacing ratio, C L is higher for the downstream cylinder than the upstream cylinder. Figure 7 compares the lift power spectra under different spacing ratios, where N1 denotes the single cylinder case. Table 4 presents the lift frequency and Strouhal number for the upstream cylinder, where S/D = 0 indicates the single cylinder case. When the spacing ratio is small (S/D = 1.2), the lift power spectrum shows a bimodal structure with smaller peak values. Compared with the peak in the single cylinder case, there exists a higher peak on the left side and a lower peak on the right side, corresponding to different Strouhal numbers of 0.2148 and 0.2476. As the spacing ratio increases to 2.0, the lift power spectrum returns to a single-peak structure, and the peak value is higher than the bimodal peaks of the S/D = 1.2 case. But it is lower than that of the peak in the single cylinder case, and the frequency is reduced. When the spacing ratio ranges from 3.0 to 5.0, the single peak shifts to the left side, and the peak frequency and peak value both decrease.   Figures 8 and 9 show the evolution of R and ω for a two-phase flow past double-parallel circular cylinders. The flow field velocity reaches equilibrium at a Reynolds number of 120. Since the working conditions are similar for parallel cylinders, the interaction between the columns depends mainly on the disturbance of the flow in the wake region with vortex shedding. At a small spacing ratio (S/D = 1.2), the wake flow behind columns forms a single vortex street, followed by periodic vortex shedding generated behind the columns. The vortex shedding on the left side rotates in the clockwise direction, and the vortex shedding on the right side rotates in the counterclockwise direction, forming two rows of counter-rotating staggered vortices on both sides. As the spacing ratio increases (S/D = 1.6), the vortices become coupled and entangled with each other, and the gap flow from the narrow space is enhanced, then intermittently biased to the left or right column. The column closer to the biased gap flow has a narrower wake flow, causing a higher frequency of vortex shedding, and greater drag and lift. The column away from the biased gap flow has a wider wake flow, so its vortex shedding frequency, drag and lift are smaller. The unsteady characteristics of the flow field and vortex field are obvious. With the increase in the spacing ratio (S/D = 2, 2.5 and 3), the perturbation between the columns decreases and vortex street shedding is formed after both double cylinders. The phases of the two vortex streets represent the characteristics of synchronous in-phase or synchronous anti-phase, where the vortices with synchronous in-phase tend to form a vortex street, while vortices with synchronous anti-phase keep the original state, propagating to the rear field. The gas phase is mainly gathered in the vortex region behind the cylinder, and the gathering area is affected by the joint influence of the fluid domain, vortex position and vortex development. Since the pressure behind the cylinder is more uniform and the vortex shedding area is large enough, the gas-liquid two-phase distribution is more balanced. With the increase in the spacing ratio, the  gas-liquid two-phase separation is more obvious owing to the full development of the vortex. Figure 10 shows the drag and lift evolution process and the lift spectrum analysis of the two-phase flow past double-parallel circular cylinders. FL D and FL L represent the drag and lift of the left column, which are shown as blue solid lines in the figure. FR D and FR L denote the drag and lift of the right column, which are depicted as red dashed lines in the figure. As can be seen in the graph, at different spacing ratios, the drag of both left and right cylinders produces certain oscillations with comparable magnitude and pattern. As for the lift, when the spacing ratio is between 1.2 and 2.0, the lift is slightly larger than the drag. When the spacing ratio is in the range of 2.0-3.0, the oscillation amplitudes of the lift and drag are comparable. Fast Fourier transform spectrum analysis of the lift force for both side cylinders shows that the oscillation frequencies are basically similar at different spacing ratios. In particular, there is no obvious main frequency when the spacing ratio is 3.0. There exist certain differences in the spectrum peak values for the left and right columns, where the main frequency peak values of the right column are slightly higher than those of the left columns when the spacing ratio is between 1.2 and 2.0, being 2.25%, 10.77% and 2.86% higher, respectively. In contrast, when the spacing ratio is between 2.5 and 3.0, the main frequency peak values of the left column are higher than those of the right columns, by 17.96% and 14.93%, respectively. Table 5 indicates the drag and lift coefficients of each cylinder at different spacing ratios, where L1.2 and R1.2 denote the left and right columns, respectively, at a spacing ratio of 1.2. Other groups have similar meanings, and N1 denotes the single cylinder case. Figures 11 and 12 compare the statistical characteristic values of the drag and lift for each cylinder at different spacing ratios. In Figure 11(a), the time-averaged drag coefficients C D of both left and right cylinders are lower than that of a single cylinder. The difference in the values of the left and right cylinders is small with different spacing ratios. In Figure 11(b), the oscillation drag coefficients C D of both left and right cylinders are higher than that of a single cylinder, and C D increases significantly when the spacing is small. With the increase in the spacing ratio, the values of the left and right cylinders gradually move close to those of a single cylinder.  In Figure 12(a), the time-averaged lift coefficient C L of the double-parallel cylinders is approximately antisymmetrical. With the decrease in the spacing ratio, the magnitude of the lift coefficient increases. As for the flow past a single cylinder, the oscillation lift is periodic and the mean value approaches zero. For the parallel cylinders with smaller spacing ratio, the time-averaged pressure distribution on the surface of column no longer remains symmetrical, and the degree of mutual repulsion between the columns is stronger, so the value of C L is not zero. As the spacing ratio increases, both left and right cylinders gradually approach the single cylinder case, and the magnitude of C L gradually decreases. In Figure 12(b), the oscillation lift coefficient C L of parallel cylinders is     lower than that of the single column, and the magnitude first decreases and then increases with the increase in the spacing ratio.  Figure 13 compares the lift power spectra of the left and right cylinders, where N1 denotes the single cylinder case. Table 6 presents the lift frequency and Strouhal number of the cylinders, where S/D = 0 indicates the single cylinder case. When the spacing ratio is small (S/D = 1.2), a bimodal structure is exhibited, which corresponds to different Strouhal numbers. As the spacing ratio increases (S/D = 1.6), the position of the double peaks moves to the right side. Consequently, the corresponding oscillation frequency and Strouhal number rise a little. When the spacing ratio reaches 2.0, the power spectrum changes to a single-peak structure, and the peak frequency and value increase slightly. As the spacing ratio continues to increase, the main frequency of the peak rises slightly as the peak value decreases.

Two-phase flow through in-line tube bundles
Figures 14-16 present the R, S and ω distributions of the two-phase flow through in-line tube bundles for different longitudinal spacing ratios. As can be seen in Figure 14, when the longitudinal spacing ratio reaches 2.0, the separated boundary layer of the front row cylinder is attached to the back row cylinders. However, owing to the presence and interaction of multiple columns, the full development of the wake is restricted. Consequently, no obvious vortex shedding is formed. As can be seen in Figure 15, when the longitudinal spacing ratio reaches 4.0, an oscillating jet with vortex shedding is formed. When the spacing ratio is further increased, as shown in Figure 16, a deflecting jet forms at the rear of the cylinders. There are discrete distribution regions with positive and negative vortex variables, which spread to the rear regions with the motion of the bubbles. Owing to the enhanced velocity distribution inhomogeneity in the narrow flow channel among the tube bundles, the velocity in the channel center increases and the pressure decreases, so that the gas phase accumulates in the low-pressure region. Affected by the effect of shedding vortex, the volume fraction of gas phase near this region increases, while the gas phase near the walls on both sides decreases. There are also regions where the gas phase gathers near the downstream bundles, but the location and size of the gathering regions are affected by the way in which the tube bundles are arranged. With the increase in the spacing ratio, the gas-liquid phase separation is more obvious owing to the full development of the vortex.
For the convenience of comparison, the individual cylinders are numbered. The cylinder at the upper and left corner is named N1, the columns are labeled in order from left to right and from top to bottom, and the rules for column labeling are the same for each group of numerical cases. Table 7 shows the averaged values of the cylinders in the same row. For instance, P-1 indicates the average value of three cylinders in the first row. The drag and lift coefficients of each column of in-line tube bundles at different spacing ratios are shown in Table 8. From the data in the table, it is clear that the drag of the front row is greater than that of the back row, and the C D of P-1 is 2.33 times that of P-9. As the spacing ratio increases, the overall drag and lift coefficient gradually increases owing to the reduction in wake interference among the bundles.
The variation of the drag and lift is analyzed for the in-line tube bundles with a spacing ratio of 4.0. The statistical characteristic values of lift and drag coefficients for each cylinder are shown in Figures 17 and 18. The timeaveraged value and oscillation value of drag located at the center of horizontal direction (N-2, N-5, N-8, N-11 and N-14) are lower than those of the cylinders on both sides. The value of C D of the center-positioned cylinder is only 93.73%, 89.40%, 92.59%, 87.05% and 86.02% of those on both sides, respectively. This difference is more obvious in the middle and rear positions because the drag of the front columns is higher than that of the rear columns.
For the time-averaged lift, it can be seen that the directions are opposite for the left and right side columns. As for C L, the oscillation values of lift for the cylinders located at the center (N-5, N-8, N-11 and N-14) are lower than those of the cylinders on the left and right sides. The value of C L at the center is only 65.23%, 56.86%, 56.31% and 72.22% of the C L on both sides, respectively.  Cylinder number with different longitudinal spacing ratios. The longitudinal spacing ratio is 4.0 in Figure 19, and there exists an oscillating jet with vortex shedding. There are positive and negative discrete distribution regions on each side of the columns, which spread to the rear region with the motion of bubbles. When the spacing ratio is expanded to 12.0, as shown in Figure 20, the form of the wake flow is a deflected jet. The vortices behind the columns are coupled and entangled with each other, and the gap flow between the narrow space is enhanced. The gap flow is  intermittently deflected to the left or right side of the columns. Because of the higher flow velocity of the narrow channel among the tube bundles, the gas phase is accelerated and converges to the middle region. As the vortex velocity between the columns increases, the gas gathers toward the vortex center, so the gas phase on both sides decreases accordingly. As the spacing ratio increases, the gas-liquid two-phase separation becomes more pronounced owing to the full development of the vortex.

Two-phase flow through staggered tube bundles
The drag and lift coefficients for each column in the staggered tube bundles at different spacing ratios are shown in Table 9. Table 10 presents the averaged values for the same row. For instance, P-1 indicates the average value for the three columns in the first row. From the data analysis in the table, it can be seen that the odd rows with three columns have higher average values of time-averaged drag and oscillating drag, compared with the even rows with two columns. As the spacing ratio increases, the overall drag and lift coefficient gradually   increases owing to the reduction in wake interference among the tube bundles.
The change in drag and lift is analyzed taking the spacing ratio of 12.0 as an example. Figure 21 shows the variation of lift and drag for the first row (N-1, N-2, N-3) and the second row (N-4, N-5). From the graph, it can be seen that the drag of the first row is slightly higher than that of the second row, in general. With regard to the lift, the difference between the two rows is relatively small.
The statistical characteristic values of the drag and lift for each column (from N-1 to N-10) are shown in Figures 22 and 23. The columns located in the second and fourth rows (N-4, N-5, N-9 and N-10) are characterized by lower time-averaged drag and oscillating drag. The variation of time-averaged lift and oscillating lift is similar. This is due to the fact that the cylinders in the even rows are located in the wake of columns in the odd rows, which has reduced the drag force. For the cylinders in the first and third rows, the time-averaged lift on the left side is opposite to that on the right side. The time-averaged lift of the middle column is smaller, while the oscillating lift is higher than that of the left and right sides. For the second and fourth rows, the direction of the timeaveraged lift is opposite for the left and right column, but the magnitude of the oscillating lift is close.

Conclusions
In this paper, for the purpose of investigating the twophase flow problem in a thermal hydraulic application, thorough numerical calculations are performed, adopting the multiphase lattice Boltzmann method. By analyzing the distribution of vortex variables, the vortex shedding law, time-averaged and oscillating drag, lift  coefficient and lift power spectrum with various spacing ratios, the following conclusions are obtained.
(1) With the increase in spacing ratio, the mutual interference between the double-tandem cylinders is weakened, and the pattern of vortex shedding changes from a single vortex street or deflected gap flow to double vortex streets. The drag and lift force of the downstream cylinder increases, and the difference of the time-averaged drag from the upstream and downstream cylinder decreases. The lift power spectrum changes from a double-peak distribution to a single-peak structure, with both the peak frequency and value decreasing.
(2) As the spacing ratio increases, the amplitude of the time-averaged drag applied to the double-parallel cylinder decreases, and the difference in timeaveraged lift between adjacent cylinders decreases. The lift power spectrum changes from a double-peak distribution to a single-peak structure, with increasing peak frequency and Strouhal number, as well as decreasing peak value.
(3) For the flow through tube bundles, the vortex shedding law, the force of each cylinder and the average value of each row have been compared and analyzed. The results reveal that the time-averaged drag, oscillating drag and lift of the front row are higher than those of the back row, while the drag and lift of the center cylinder are lower than those of the cylinders on both sides. To ensure the stable operation of the reactor and other engineering equipment, attention should be paid to the structural strength of the front columns. (4) The time-averaged drag, oscillating drag and lift of the odd row (three cylinders) in the staggered tube bundles are higher than those of the even row (two cylinders), due to the fact that the latter are located in the wake of the former. For the stable operation in   the reactor and other engineering equipment, attention should be paid to the thermal hydraulic limits of odd-row cylinders.
The research in this paper is significant and meaningful for revealing the thermal hydraulics characteristics of two-phase flow, and ensuring the stable operation of the reactor or other engineering equipment. It will also contribute toward laying the foundation for subsequent related research.