Liutex-based analysis of drag force and vortex in two-phase flow past 2-D square obstacle using LBM on GPU ＊

Peng-xin Cheng, Nan Gui, Xing-tuan Yang, Ji-yuan Tu, , Sheng-yao Jiang, Hai-jun Jia

1. Institute of Nuclear and New Energy Technology, Collaborative Innovation Center of Advanced Nuclear Energy Technology, Key Laboratory of Advanced Reactor Engineering and Safety of Ministry of Education, Tsinghua University, Beijing 100084, China

2. School of Engineering, RMIT University, Melbourne, Australia

Abstract: The two-phase flow past a square is a ubiquitous phenomenon widely encountered in industries and engineering, where the interaction of disparate phases coupled with the influence of the solid is rather complicated. In this context, the flow characteristics and the vortex field are investigated to reveal the mechanisms of the two-phase drag and the vortex variation. The lattice Boltzmann method (LBM) is utilized to study the multi-component two-phase flow. The computation implemented on the GPU is remarkably accelerated thanks to the natural parallelism of the LBM. The process of the two-phase flow past a square is thoroughly examined. The drag and lift forces, including the total force and the components caused by the dispersed phase and the continuous phase, respectively, are obtained, and their variation mechanisms are explained. Meanwhile, the vortex-identification approaches based on the Liutex as well as the traditional methods are compared. The relationship between the bubble breakup and coalescence processes and the extremums of different vortex identification variables is analyzed.

Key words: Two-phase flow, flow past a square, lattice Boltzmann method (LBM), vortex identification, Liutex, two-phase drag

Introduction

The flow past obstacles has a wide application in industries ranging from chemical process, food processing, power engineering, refrigeration to metallurgy. For decades, numerous experiments and numerical simulations were carried out focusing on the drag and lift forces, the vortex separation, the flow-induced vibration, the fluidelastic instability and so on. The two-phase flow past obstacles plays a significant role in practical applications, such as the packed bed reactor in petroleum chemical process and the heat exchangers in power industry, including evaporators, condensers, reboilers and steam-generators.

In the field of nuclear engineering, the flow between bundles in boiling-water reactors is a typical scenario of two-phase flow past cylinders. The vortex shedding and the flow-induced vibration may lead to severe damages in the forms of fatigue and fretting of materials. In addition, the water-cooled pebble bed reactor is a promising new concept reactor under research, which combines the advantage of the desirable economy features of the water reactor and the reliable safety performance of the spherical fuel element. In this case the coolant experiences a transition from the single-phase flow to the two-phase flow past the pebble bed. It is essential to have a thorough understanding of this issue to determine the thermal-hydraulic characteristics.

Nevertheless, the two-phase flow past obstacles involves three phases, namely, solid, liquid and gas,which are coupled and in an intricate state, due to the impact of the bubbles on the solid, the influence of the bubble movement on the vortex shedding, the involvement of the bubbles in the vortex, and some other factors. Under such circumstance, using the traditional computational fluid dynamics (CFD)methods like the finite volume method (FVM) we might encounter difficulties and fail to acquire the details of the flow field in the macroscopic scale.

Endowed with the edge of the flexible geometry characteristics, the inherent parallelism, the high precision and the simplicity of implementation, the lattice Boltzmann method (LBM) has witnessed a rapid development in the past decades, particularly in the field of multiphase flows and flows with a complex geometry. Based on the molecular kinetic theory, the LBM provides a novel and potent approach with a solid background of physics and efficient algorithms in the multi-scale analysis in the mesoscopic perspective. In contrast to the traditional CFD methods that involves a nonlinear convection term in the Navier-Stokes equation and where the Poisson equation is solved to obtain the pressure, the convection operator is linear in the Lattice Boltzmann equations and the pressure could be obtained conveniently by solving the simple equation of state.Furthermore, with the computation largely restricted to local nodes, the LBM is ideal for the implementation of parallel algorithms and on the related hardware.

Graphics processing units (GPU) were originally designed to rendering the graphics and recently has been employed in a broad scope of domains. Modern GPUs are massively parallel computational devices with an apparent advantage over the CPUs in the field of high performance computing (HPC). The implementation of the LBM on the GPU may boost the computation speed and achieve a reliable acceleration performance.

In perspective of the basic research, the vortex is an issue of vital significance in understanding turbulence. After decades of intensive researches,there are still plenty of ambiguous issues left to be solved, including a widely-acknowledged problem of the definition of the vortex. The vortex identification methods have evolved from the first generation intuitive vorticity-based methods, to the second generation methods that rely on eigenvalue-based parameters includingQ, Δ,2λ, andciλ, and finally to the third generation Liutex method. The first generation method fails to discriminate the rotational motion from a shear boundary layer in viscous turbulent flows. The maximum vorticity magnitude does not necessarily imply a vortex core and the local vorticity vector is often not aligned with the vortex core direction. The second generation methods are contaminated by shearing, thus they can not be used to capture the vortical strength to some extent. Since the parameters are all scalars, involving no directions to identify the rotation axis. In addition, they require the use of iso-surfaces to display the vortex structure and depend on some case-related uncertain thresholds.

Recently Liu and Gao et al.[1-2]put forward the third generation Liutex method, which provides a perspective of defining the vortex in a more precise way. According to the method, the vorticity could be decomposed into the non-rotational part, mainly the shearSand the rotational part called LiutexR,namelyω=S+R. The direction of the Liutex vector, representing the local rotation axis, is defined via the real eigenvector of the velocity gradient tensor.The magnitude of the Liutex, used to quantify the local rotational strength, could be simply determined by an explicit formula proposed by Wang et al.[3],which follows the idea of getting rid of the non-rotational part from the vorticity. Liu et al.[4]extensively compared the Liutex and the Omega-Liutex with previous vortex identification methods, highlighted the essential advantage of the third generation methods, and it was shown that the core issues of the vortex identification could thus be mathematically defined. The Liutex was applied in various cases, including the hairpin vortices[5]and the vortices in swirling jets[6-7], the open channel flow[8]and the turbine rotor passages[9].

This paper will present a comprehensive investigation of the flow and vortex field and the drag and lift forces of multiphase flows past obstacles. The variation of the multiphase drag force in the movement of the fluid will be discussed and interpreted. The detailed evolution of the primary variables during the bubble deformation will be analyzed for its relationship with the vortex generation and evolution. A wide range of identification parameters including the Liutex will be evaluated. The relationship between the extremums of different vortex identification variables and the bubble deformation process, mainly the breakup as well as the coalescence, will be analyzed.

1. Numerical method models and methods

1.1 Basics of lattice Boltzmann method

The LBM was originated from the lattice cellular automaton and the lattice gas model. One could find its root in the Boltzmann equation based on the kinetic molecular dynamics.

wheref(x,ξ,t) denotes the particle distribution function also called thepopulation, representing the density of the particles with the velocityξat the positionxand the timet. The collision operatorΩfcharacterizes the microscopic interaction between molecules and involves a nonlinear integral, which could be simplified by the Bhatnagar-Gross-Krook(BGK) single-relaxation-time approximation in the form of a linearized expression[10].

whereτis the relaxation time related to the transport properties such as the viscosity. In addition to the BGK model, other models, including the multiple relaxation time (MRT) model, were proposed to enhance the numerical stability.

Similar to other traditional CFD methods, the LBM inherits the method of space and time discretizations, yet the LBM distinguishes itself in discretizing the velocity space, with the zero velocity and a finite scope of directions linking a node and its closer neighbors. The widely adopted discretization model is usually named as DdQm, wheredandmdenote the numbers of dimensions and velocities.In the present simulation, the prevailing D2Q9 model is adopted and illustrated in Fig. 1.

Fig. 1 (Color online) The velocity set of D2Q9 scheme

Through discretizing the physical space, the velocity space and the time, the discrete lattice Boltzmann equation (LBE) could be expressed as

wherecirepresents the discrete velocity direction,andfi(x,t) is the particle distribution function along the direction ofci.

The equilibrium distribution function can be described as follows

wherewiis a set of weights corresponding tociandcsis the sound speed

By means of the multi-scale expansion techniques such as the Chapman-Enskog analysis, the LBGK model could recover the continuity equation and the Navier-Stokes equation in the macroscopic scale in the limit of a very small Mach number, i.e.

The macroscopic parameters such as the density,the velocity and the internal energy can be calculated statistically from the moments of the discrete distribution function.

1.2 Multiphase lattice Boltzmann methods

For the multiphase flow simulation, the Shan-Chen pseudo-potential method proposed by Shan and Chen[11]is the most widely applied. This model incorporates a repulsive or attractive Shan-Chen force without the need of any explicit interface tracking, and it is efficient and is adopted in the present simulation.

For the single-component multiphase flow, the interaction force density acting on the fluid atxcan be expressed as an integral over all possible interaction sitesx′.

whereψ(ρ) also called the pseudo-potential,denotes the effective density, and its prevailing form is proportional to the effective mass as follows

where the reference density0ρis usually set to unity andψ()ρis close toρfor0ρ?ρ. As another common form, it is equal to the fluid density itself,namely,ψ(ρ)=ρ.G(x,x′) is a Green function to determine the strength and the range of the interaction,and in the simplest form for considering the nearest lattice neighbors, it is expressed as

The discretized Shan-Chen force can be summed over the neighbors.

In a multiphase fluid system constituted byncomponents, the LBE of the Shan-Chen model can be expressed as

The total interaction force acting on theσcomponent can be computed by summing the forces imposed by all components.

In the present work we focus on the two-component multiphase system without phase change, namely,n=2 andGσσ= 0.Gσσ′is set as positive for the repulsive force between different components. The overall density, the weighted average velocity and the pressure of the fluid are calculated as follows:

For the force scheme, the exact difference method (EDM) is adopted based on the kinetic theory,where the force is introduced via shifting thefiin the velocity space.

1.3 Momentum exchange algorithm of the bounce-back method

The widely adopted method to implement the no-slip wall boundary condition in the LBM is the bounce-back method. It originates from the simple concept that the fluid particle streaming towards the wall would be reflected back to the node where it comes from, and the method is simple to implement and ensures no flux across the boundary.

In general, the force acting on an obstacle could be calculated via an integral of the stress tensor t on the surface.

wherewσrepresents the stress tensor andnis the unit normal vector that points from the surface to the fluid. Unfortunately, in the Bounce-Back scheme, this method suffers from the loss of accuracy. Hence the momentum exchange algorithm (MEA) is preferable in the LBM and is adopted in this paper.

The basic idea of the MEA is to identify the populations across the boundary and calculate the net momentum transfer consisting of two parts, namely,the momentum carried by the fluid to the walland the momentum transported from the wall to the fluid. A schematic diagram is demonstrated in Fig. 2. The procedure includes the following steps[12].

(1) Identify the momentum links between the solid and boundary nodes.

(2) Evaluate the incoming populationsand the bouncing populationsof each momentum link.

(3) Calculate the transferred momentum during the time steps and the acting force on the boundary.

Fig. 2 (Color online) The schematic diagram of bounce-back method

The drag and lift forces are calculated using the MEA in the following section. The drag forceFDquantifies the resistance in the opposite direction of the flow. The lift forceFLis in the direction perpendicular to the flow direction. Thus they can be evaluated by the components of the total force inXandYdirections. Since the force is calculated via the populationf, for the multi-component case in the present paper, the force of theσcomponent could be obtained based on the population of each componentfG.

1.4 Parallel implementation of LBM on GPU

For practical convenience, the main process of the calculation of the LBE can be decomposed into two parts, namely, the collision and the streaming.

In the collision step, the collision is implemented in the local lattice site, which brings about non-linearity, called the local non-linearity, because only the local calculation is involved . In the streaming step, the movement of the particles towards neighboring lattice nodes is considered, which requires the next neighbor interaction, which is apparently linear, called the linear non-locality. These characteristics render the LBM naturally suitable for the computation on parallel architectures including the GPU.

The GPU is initially intended for graphics applications and is currently widely applied in scientific computations. The CUDA is specifically developed to implement the parallel computation on the GPU by the Nvidia Corporation with a rapid progress in recent decades. In the CUDA architecture,a heterogeneous programming model is adopted. The CPU is regarded as the host optimized for the single instruction, single data (SISD). The GPU is viewed as adevice, capable of the single instruction, multiple thread (SIMT), i.e., executing a multitude of threads in parallel simultaneously, to take advantage of thousands of cores in the GPU and to remarkably boost the computation.

The implementation of the LBM on the GPU has attracted a wide research interest and is proved successful in achieving a reliable acceleration and is ideal for big data applications. In this paper, the open-source software Sailfish is adopted in the simulation[13]and the computation is implemented on the NVIDIA Tesla K20 GPU.

1.5 Vortex identification methods and Liutex

Several well-acknowledged identification parameters are selected including the vorticityω[14],Q[15], LiutexR[5], OmegaΩ[16-17], Omega-LiutexRΩ[4], which are defined as:

2. Validation

In this section, the LBM model is firstly validated by a 2-D single-phase flow past a square cylinder for the accuracy of the force evaluation. Then the Laplace test is conducted to examine the multiphase model. Breuer et al.[18]simulated the flow past a square cylinder by the LBM as well as by the FVM and their results are adopted as the reference.The comparison of the present simulation results with the validation reference is displayed in Fig. 3 in logarithmic coordinates and the data are listed in Table 1. On the whole, the simulation results are in desirable conformity with the reference. In the 2-D case, with the characteristic length corresponding to the obstacle’s diameter, the drag coefficient is calculated as follows

Fig. 3 (Color online) The validation of drag coefficient of flow past a 2-D square cylinder

Table 1 The validation of drag coefficient of flow past a 2-D square cylinder

In the multi-component Laplace test, a circular droplet of one phase with radiusRis initially surrounded by another liquid phase. The interface is determined to be where the density reaches the average of the bulk densities inside and outside the droplet. The result of the Laplace test is shown in Fig.4, where the relationship between the measured pressure difference Δpand 1/Ris linear, with the coefficient of determination of the fitted lines above 0.999, in an excellent agreement with the Laplace Law, namely Δp=σ/R.

Fig. 4 (Color online) The validation of Laplace test for multicomponent case

The two-phase model validation was completed in our earlier work[19], by validating the bubble departure diameter and the release period.Additionally, the GPU performance of the LBM is also tested and validated via a cavity driven flow case in another earlier work[20]. As a whole, the LBM mode,the two-phase scheme and the GPU performance are well validated.

3. Results and discussions

3.1 Numerical setup

The flow domain is discretized by 256×1 024 Cartesian mesh grids. The square obstacle consisting of 60×60 nodes is placed right in the center of the channel and is referred to as the block or the square in the following sections. Initially a circular bubble of the dispersed phase is set up at the upside of the channel and it is of 100 nodes in diameter, which will be called the bubble later. The rest of the field is filled with the continuous phase. Note that the density contour denotes the density of the continuous phase.Based on the previous multi-component multiphase research[21], the viscosity of the two phases are set as the same and the parameterGis determined viaGρ= 1.8. The positive directions of the initial velocity of the bubble and the force is defined as in the downward direction

3.2 The evolution of two-phase flow

The evolution of the density contour is shown in Fig. 5, where the typical cases are selected. Under the influence of the initial velocity and the imposed force,the bubble gradually moves towards the square. As the nearly round bubble approaches the obstacle, the reaction from the obstacle flattens the bubble at the bottom. Thereafter it is separated by the solid (atT=7 300) and stretches a great extent (atT=10 200), till the filament between the bubble becomes unstable and breaks up into several satellite bubbles (atT=10500). Two small bubbles adjacent to the side move along the lateral periphery of the square (atT=11900) and then get involved in the vortex behind the obstacle (atT=14 700). With the distance between the two bubbles reduced, they finally merge into one bubble (atT=15 400). The bubble on the top of the square moves along the right side of the square and detaches from it afterwards.

3.3 The evolution of drag force and lift force

The evolutions of the drag force and the lift force as well as the detailed distribution of the forces of the disparate phases are analyzed in this section.

Fig. 5 (Color online) The contour of density at different time steps

Firstly the evolution of the drag force is displayed in Fig. 6(a). As the bubble approaches the cylinder, the drag force of the continuous phaseincreases sharply while that of the dispersed phasedecreases correspondingly, resulting in a net increase of the total drag force. Note that the decrease of the drag force componentcaused by the dispersed phase in the process of the two-phase flow past the square is as compared to the earlier stage process, regarded as a locally continuous-phase flow past over the square since the bubble is far away in the upstream of the square. To say specifically, when the bubble is far away, the fluid around the square is occupied by a flow of continuous phase, the density of which is larger than that of a dispersed-phase flow past the square when the bubble reaches the square.

Fig. 6(a) (Color online) The evolution of drag force

On the other hand, compared to the locally single-phase flow past the square when the bubble is far away from the square, the movement of the bubble in the locally two-phase flow past the square compresses the fluid between the bubble and the solid.Also the velocity of the continuous-phase around the square in the two-phase flow is larger than that in the single-phase flow past the square. The reasons include the following: (1) when the bubble reaches the square,the volume of the dispersed phase occupies a large part of the region around the square, (2) because of the mass flow conservation, the total mass flow rate of the locally two-phase flow (equivalent to the two-phase density) past the square should be equal to that of the continuous phase alone, which is a locally single-phase flow when the bubble is far away in the upstream, (3) due to the effects of (1) and (2), the velocity of the continuous-phase in the two-phase flow past the square must be larger than that of the continuous-phase when the bubble is far away from the square, which is regarded as a locally single-phase flow past the square. As the drag force is proportional to the square of the velocityFD～u2, a larger flow velocity must imply a larger momenta of the continuous phase, as well as a larger continuous-phase component of the drag force. That explains the mechanism why the drag force componentincreases sharply in the two-phase flow past the square, as well as the total effective drag forcesince the drag force of the continuous phaseincreases dominantly, which boosts.

After that, the curves flatten out (atT=7 301 and atT=7 449), corresponding to the nearly constant amount of the fluid surrounded by the bubble.Subsequentlygradually increases and finally reaches the local maximum (atT=10 547). In the meantime,decreases till it reaches the former minimum (atT=10 222) and the latter minimum (atT=10561). The curve oftakes the minimum atT=10 220. The extremums coincide with those in the break-up progress of the bubble. Then the curves have a similar plateau with the remaining small bubble settled on the top of the block. Afterwardsgradually rises as the bubble starts to move.andreach a local maximum (atT=14 741) when the bubble detaches from the solid.

In the evolution of the lift force shown in Fig.6(b), the forces of both phases are fixed at the value of zero initially related to the symmetric flow. As the bubble gradually comes into contact with the solid, the lift force of the continuous phaseturns into positive while that of the dispersed phasebecomes negative. Then at the moment of the bubble breakup,reaches its maximum (atT=10163)andreaches its minimum (atT=10174), with the total lift forcereaching its maximum (atT=10154). Since the parameters of the breakup process are not perfectly symmetric, such as the location of the bubble fracture, the portion of the volume of each bubble, the position and the motion of each small bubble, the asymmetry will make the lift force not in balance in the horizontal direction.Subsequently as the two small bubbles on the lateral side of the square move,grows into positive and reaches its maximum (atT=11194) when the bubbles are about to detach from the solid. Similarly because of the asymmetry in the motion of the bubbles on two sides as well as the unequal distances between the solid and each bubble, the lift force is not balanced and a net value will be generated.

Fig. 6(b) (Color online) The evolution of lift force

Finally the asymmetry comes mainly from the last remaining bubble resting on the center top of the square, which is in an unstable equilibrium. The fluctuations in the flow would contribute to the destabilization and the movement of the bubble towards either side, which might further aggravate the asymmetry since only one bubble is on the lateral side of the square. The movement of the remaining bubble makes the lift force fluctuating. As the bubble detaches from the upper periphery of the block,reaches its maximum (atT=14698). Afterwards,rises sharply with only one bubble on the right side of the square whileanddrop abruptly. Until the bubble detaches from the influence area of the square, the degree of asymmetry diminishes and the lift force falls back. The moment when the last bubble departs from the lateral area of the block,reaches its peak atT=15349 and then drops, whileandreach the minimum atT=15384 andT=15418, respectively.

3.4 The evolution of R, S, ω and comparison of vortex identification methods

Based on the relationR+S=ω, the evolution ofR,Sand the vorticityωat different time steps are shown in Fig. 7. In the evolution process, as shown in Fig. 7, before the bubble comes into contact with the solid, zones with negative Liutex are observed on the top of the block, as well as zones with positive Liutex on the lateral top of the block and behind the block. Some regions with negative Liutex are also observed inside the bubble. The compression of the bubble on the adjacent fluid leads to a concentration of the Liutex. AtT=7 300, the deformation also affects the distribution of the Liutex,with a negative zone in the center and the positive zones on the sides. AtT=10 200, when the bubble breaks up, zones with relatively high magnitude of the Liutex are observed near the position of the bubble break-up. In addition, the stretched bubbles change the distribution of the Liutex remarkably. Zones with high positive Liutex are observed between the bubbles and ahead of the bubbles, as well as zones with negative Liutex laterally adjacent to the bubbles. AtT=10500, after the bubble break-up process, the surface tension deforms the bubbles into round shapes.Zones with high magnitude of the Liutex are observed in the surrounding fluid near the ends of the bubbles.AtT=15 400, when two small bubbles merge into one, as a result of the distortion, zones with high level of magnitude of the Liutex are observed inside the bubbles, as well as zones with alternatively positive and negative Liutex outside the bubbles.

We will study the changes of the Liutex across the interface between the two phases. It is worth mentioning that the Shan-Chen multiphase LBM is a diffusive interface method in contrast to the sharp interface method. In the Shan-Chen model the order parameter profiles such as that of the density vary smoothly across the interface between two bulk values.In multiphase systems, the density varies from the liquid to gas phases across the interface. The interface width represents the length scale characterizing the variations in the density profile across the interface.The interface position is usually defined by the average ofgρa(bǔ)ndlρ. The Shan-Chen model will result in a relatively wide interface and in our simulation the interface width is nearly 4-6 nodes.

Fig. 7 (Color online) The evolution of R, S and ω at T=7 300, 10 200, 10 500 and 15 400, respectively

As shown in Fig. 8, the detailed distribution of the density and the Liutex atT=7 300 (Fig. 8(a) ) is chosen as an example to see the variations of the parameters across the interface. The distributions of the density, the velocity components and the Liutex on two linesX=65 andY=545 are considered. In Figs. 8(b) and 8(e), the sharp variations of the density clearly indicate the location of interfaces, as is roughly consistent with the distribution shown in Fig. 7(a),while it is cut off by the interface zone. Near the interface, some disconnected regions are observed with positive or negative Liutex values. In Fig. 8(c),from the curves of the velocity, one sees fluctuations near the interface, probably because of the averaging process of Eq. (13) in calculating the velocity where the influence of the density variation is involved. In Figs. 8(c) and 8(f), the scattered positive or negative Liutex values are witnessed near the interfaces,possibly due to the drastic variation of the velocity across the bubble surface.

As a whole,SandΩare in a good consistency with each other, both in distributions and magnitudes. The bubble’s motion and shape dominate the variations of the shearSand the vorticityΩ.Their distributions and magnitudes are so close to each other that the variations ofRplays only a secondary role. In other words, because of the rela-tively low velocity, the shear deformation is dominant in the vortex evolution characteristics of the two-phase flow past the square. Nevertheless, as is indicated byR, the pure rotational motion (the real role of the vortex) also has clear features dominated by the formation of bubbles, their deformation,breakup and coalescence.

Fig. 8 (Color online) The distribution of density, velocity and R parameters at T =7300

Fig. 9 (Color online) The contour of R, Q, ω, S, Ω and RΩvariables at =10150 T

Moreover, taking the very moment of the bubble breakup for example, Fig. 9 shows a comparison of the varied vortex identification approaches defined by Eq. (18). Compared to the density fieldρ, allR,Q,ω,S,ΩandRΩshow some basic characteristics of the vortex structure, with, however,differences.

It is clear that, at the moment of the bubble breakup, the rotational vortices are formed around the three small satellite bubbles above the square.Therefore, there are mainly three rotational vortices in the left, middle and right sides of the square top.Meanwhile, in the downside of the square, some small pure rotational vortices are also formed. The rotational vortices are mainly formed around the locations with sufficiently large or locally maximum curvatures of the interface, i.e., around the ends of the stretched bubbles. It is because the vortex can only be generated when and where the shear of the fluids is so large that the fluid cannot resist the shear deformation.Therefore, during the bubble breakup into an absolutely different geometrical topology or near the region with sufficiently large curvatures of the interface, vortices may be generated.

In addition, although with different vortex identification methods, similar structures of the vortex are obtained,RandQ, andΩandRΩgive comparatively similar distributions, respectively, in two different pairs. The reason is thatRandQevaluate the vortex by the absolutely rotational strength, whereasΩandRΩdetermine the vortex by relative levels of strengths of the rotation against the shear.

3.5 The evolution of R max and S max

The statistics of the primary variables (denoted by the superscripts max,min and avefor maximum, minimum and average) of the dispersed and continuous phases (distinguished by the subscriptsdandc, respectively) will be discussed in this and the following sections. The local maximums and minimums are represented by means of red and blue marks in the figures. In order to ease the burden of computation, the data are collected every ten time steps, resulting in the exhibited tenfold extremes.

Firstly, the maximum value ofRis shown in Fig. 10. In Fig. 10(a), it is noted that the data ofare only selected from the interior of the bubble(dispersed phase), whereas in Fig. 10(b) the data ofare selected from the exterior of the bubble(continuous phase). In Fig. 10(b), the suddenly sharp increase ofat aboutT=10150 is caused by the fact that the bubble begins to break up into small satellite bubbles at the very moment, whereas the sharp increase ofatT=15300 corresponds to the moment when two small bubbles coalesce into a large one (recalling Fig. 7 atT=10 200 and

T=15 400). In Fig. 10(a) for the dispersed phase, the sharp increase ofcorrespondingly coincides with that of the continuous phase in Fig. 10(b). The small differences of the peak times ofandin Figs. 10(a) and 10(b) are caused by the tiny differences of the rotational motion of the vortices of the two phases.

Fig. 10(a) (Color online) The evolution of

Fig. 10(b) (Color online) The evolution of

On the other hand, with regard to the maximum of the shear, the evolution ofSmaxis shown in Fig.11. In Fig. 11(a) for the dispersed phase, the two remarkable peaks are atT=10380 and atT=15330, with the locations approximately indicated in Fig. 11(b). The two relatively lower peaks are atT=8 440 and atT=14 670. In Fig. 11(b) for the continuous phase, the curve has four peaks atT=8 420,T=10170,T=10530 andT=14 670.Comparing Figs. 11(a) and 11(b) to Figs. 10(a) and 10(b), it is also clear that the peaks ofRmaxandSmaxalmost take place at the same moments. This confirms the former conclusion that the peak ofmaxRor the occurrence of the rotational vortex is caused by the condition that the deformation is too large for the fluids to resist even more shear. In the traditional applications, the vortex is always identified by the vorticity, and in the current situation, the flow is mostly dominated by the shear rather than the pure rotation. It is because the flow is very slow and the Reynolds number is very small. When the Reynolds number is large, the vortex shedding will appear and the Karman vortex street will turn unstable, under which circumstance the Liutex will become dominant.However, as shown in Fig. 10, the maximum of the real vortex indicated bymaxRcan also be identified by the maximum shearmaxS, and the vorticity identification methods are somewhat usable if a very accurate description of the vortex is not required.Nevertheless, usingRto quantify the pure rotational vortex is mostly encouraged in two-phase flows, either for the vortices in the dispersed phase or in the continuous phase.

Fig. 11(a) (Color online) The evolution of

Fig. 11(b) (Color online) The evolution of

3.6 The evolution of maxQ

The evolution ofQmaxis shown in Fig. 12.With respect tofor the dispersed phase in Fig.12(a), three remarkable peaks are observed atT=10 250,T=14 730 andT=15300. Referring to the previous figures, these peaks are near the points of the bubble break-up and coalescence. In Fig. 12(b)for the continuous phase, the three major peaks are atT=6850,T=10140,T=14 460 andT=15310,which cover the peaks of.

Moreover, a detailed inspection also indicates that the peaks ofcome later than those of.It may possibly mean that the rotational vortex or the vortical motion is transferred from the exterior to the interior of the bubble. In other words, the vortex inside the bubble is induced by the external shear or rotational motions.

Fig. 12(a) (Color online) The evolution of

Fig. 12(b) (Color online) The evolution of

To sum up, it is found thatQmaxcould also be used to indicate the two-phase vortical motion. In other words, it is also an applicable indicator for the shape deformation of the bubbles and the rotational motion of the fluids.

3.7 The evolution of Ω ave and

Finally, the evolution ofΩaveandare shown in Figs. 13 and 14. In this section, the averaged values ofΩandΩRover the dispersed or continuous phases are considered.in Fig. 13(a)first hikes to a maximum and subsequently falls back,later shows drastic fluctuations. This kind of variations ofclearly indicates the variations of the vortical motion inside the bubbles, clearly more smooth and slower than the variations ofin Fig.13(b). As a whole, neithernorcan reflect the instantaneous variations of the vortex motion of the two-phases as good as the maximum values in the vortex identification methods mentioned in the previous sections.

Fig. 13(a) (Color online) The evolution of

Fig. 13(b) (Color online) The evolution of

Figures 14(a) and 14(b) show the variations ofof the dispersed and continuous phases. Under this condition, the variation ofshows much clearer features than that ofΩave. To say specifically,always fluctuates more evidently and intermittently thandoes. Moreover, the peaks ofare also around the points of suddenly sharp changes obtained by some other vortex identification methods, e.g.,RmaxorQmax. For example, the peak atT=10 480 corresponds to the bubble breakup,and the peak atT=15 230 is related to the bubble coalescence. In other words, we can say thatis a much better choice to use thanΩave.

The reasons are as follows: in general, the averaged values of the vortex identification may be not a good index to characterize the vortex structures,since the locally larger rotation of the fluid vortex may become weaker and obscured by the averaging process. That is whyaveΩin both Figs. 13(a) and 13(b) fails to show the abrupt change of the vortex related to the bubble breakup or coalescence. In contrast, although averaged by respective phase fields in the same way,of both the dispersed and continuous phases can still indicate the abrupt change of the vortex via evident peaks. Referring to Eq. (18),it is suggested that this mechanism lies in the definition of. In general,αandβare better thanAandBin defining the pure rotational motion of the fluids.αandβare used to defineRor the pure rotational vortical motion. On the contrary,AandB, for the shear deformation, do not indicate the pure rotational motion features.Therefore, it is better to evaluate even the total level or strength of the rotational motion of the fluids by usingαandβ.

Fig. 14(a) (Color online) The evolution of

Fig. 14(b) (Color online) The evolution of

4. Conclusions

This study shows the characteristics of the vortex and the drag force of two-phase flows past a square obstacle using the Liutex-based analysis through the validated multiphase LBM model by the GPU computation. The main findings can be briefly summarized as follows:

(1) With respect to the evolution feature, the bubble may be deformed and stretched greatly until its final breakup during the process of the flow past the square. After that, the small satellite bubbles may coalesce to form larger ones in the downstream.

(2) Dominated by the bubble behavior, the drag force components caused by the dispersed bubble phase and the continuous fluid phase may be differentiated, with sharp increase offor the continuous phase and abrupt drop offor the dispersed phase. The mechanism of such kind of variations is explained in details by the specific variations of the velocity and the density of the two-phases. As regards the lift force, the bubbles on the lateral side of the square generate an imbalanced force and the extremums are attained at the moment of the bubble breakup or the detachment of the bubble from the square.

(3) Also dominated by the two-phase interaction behavior, the vortex identification methodRsuccessfully describes the region of the pure rotation of the fluids, although its magnitudes are small compared to the shear of the fluids. The coincidence of the distribution and the levels between the shear magnitudeSand the vorticityωshows the dominating role of the shear/deformation over the rotation in the two-phase flow.

(4) Through averaging over the respective phase fields,RΩshows clear advantages compared toΩ,sinceRΩwill not be obscured or weakened by the averaging process even in the complex two-phase flows past the obstacle, whereas theaveΩfails to show such kind of features. SinceΩandQare based on the Cauchy-Stokes decomposition and might not work as well asRΩ, whileRΩis based on the Liutex which is not contaminated by the shear.