1. Introduction
Phase separation and multiphase flow are critical phenomena across a wide range of fields, including aerospace, energy engineering, and the natural sciences [1]. In aerospace, multiphase flow is vital for investigating the fuel transfer in space vehicles, energy systems, enhanced heat transfer devices, and space stations [2]. In energy engineering, it is integral to the functioning of boiling tubes, gas-liquid separators, heat exchangers, and other equipment in nuclear stations. Additionally, multiphase flow offers valuable insights into natural phenomena such as fast-moving landslides and mudflows. Therefore, a comprehensive study of the complex phenomena associated with multiphase flow can facilitate more rational and economical design and operation of various types of industrial equipment. Numerous researchers have explored multiphase flow from various perspectives using a range of research approaches, including theoretical, experimental and numerical methods. While both experimental and theoretical methods offer distinct advantages, experimental approaches tend to be expensive and less flexible, whereas theoretical models are usually difficult to solve or under quite simple conditions. To the contrary, numerical simulation has become an essential tool for solving the complex multiphase problem, reducing experimental costs, and overcoming some application limitations of theoretical methods.
As a versatile numerical method, lattice Boltzmann method (LBM) has made significant progress in simulations of multiphase flow [3-7]. Based on the kinetic theory, LBM characterizes the flow of macroscopic fluids by calculating the distribution function of particles. It maintains continuity at the macroscopic level while achieving discretization at the microscopic level [8, 9]. The mesoscopic characteristics of LBM enable a deep analysis of the interactions between particles at the microscopic level and the flow processes of fluids at the macroscopic level. LBM for multiphase flow has developed into various models, which can be mainly classified into four categories: color-gradient [10], pseudopotential [11], free-energy [12], and mean-field models [13]. Specifically, the most notable feature of the pseudopotential model lies in the introduction of an interparticle potential. This characteristic allows the model to spontaneously realize the mixed or separated state between two phases. In 2009, Wang et al utilized the traditional pseudopotential LBM to simulate phase separation under shear flow [14]. In 2014, Lycett-Brown et al introduced a kind of multiphase central-moments-based lattice Boltzmann method (CLBM) to reduce spurious velocities [15, 16]. Then Lycett-Brown and Luo constructed a forcing scheme for the pseudopotential method to further enhance the CLBM model [17, 18]. In 2020, Fei et al implemented an improved three-dimensional CLBM to simulate the liquid-vapor phase change process [19]. In 2021, Saito et al modified the mechanical stability condition in CLBM to align with the Maxwell construction when simulating two-dimensional forced-convection boiling on a cylinder [20]. Recently, Zhang et al investigated the flow boiling in complex microchannel with U-bend using a pseudopotential model [21]. In the current study, the phase separation under shear is investigated using the pseudopotential model combined with a CLBM [22].
The rest of the paper is organized as follows. Section 2 provides a comprehensive description of the CLBM pseudopotential model, incorporating the consistent forcing scheme and the Peng-Robinson equation of state (EOS). The LBM is validated in section 3 using four benchmarks, including Poiseuille flow, freefall, shear flow, and gas-liquid coexistence. In section 4 , we simulate the phase separation under shear, followed by detailed analyses and discussions. Finally, section 5 concludes the paper.
2. Pseudopotential CLBM
In this paper, we employ the CLBM pseudopotential model with the consistent forcing scheme [20]. This model is partitioned into two main steps, which are collision:
$\begin{eqnarray}\begin{array}{rcl}|{f}_{i}^{\ast }({\boldsymbol{x}},t)\rangle & = & |{f}_{i}({\boldsymbol{x}},t)\rangle +{\boldsymbol{\Lambda }}[|{f}_{i}^{{\rm{eq}}}({\boldsymbol{x}},t)\rangle -{f}_{i}({\boldsymbol{x}},t)]\\ & & +\,{\delta }_{t}({\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{\Lambda }}}{2})|{F}_{i}({\boldsymbol{x}},t)\rangle \,\rm{,}\,\end{array}\end{eqnarray}$
and streaming: $\begin{eqnarray}{f}_{i}({\boldsymbol{x}}+{\delta }_{t}{{\boldsymbol{c}}}_{i},t+{\delta }_{t})={f}_{i}^{* }({\boldsymbol{x}},t)\,\rm{,}\,\end{eqnarray}$
where the symbol $\left.| \cdot \right\rangle $ denotes a column vector, * post-collision quantities, $\left|{F}_{i}\right\rangle ={\left[{F}_{0},{F}_{1},\ldots ,{F}_{8}\right]}^{{\rm{T}}}$ the force term, $\left|{f}_{i}\right\rangle ={\left[{f}_{0},{f}_{1},\ldots ,{f}_{8}\right]}^{{\rm{T}}}$ the distribution function, $\left|{f}_{i}^{{\rm{e}}{\rm{q}}}\right\rangle \,={\left[{f}_{0}^{{\rm{e}}{\rm{q}}},{f}_{1}^{{\rm{e}}{\rm{q}}},\ldots ,{f}_{8}^{{\rm{e}}{\rm{q}}}\right]}^{{\rm{T}}}$ the equilibrium distribution function, and I the identity matrix. The collision matrix Ʌ is defined as Ʌ = T-1ST, where S is a 9 × 9 diagonal matrix and T is the transformation matrix that will be defined later. As shown in figure 1, the discrete velocity for D2Q9 takes the form: $\begin{eqnarray}\left\{\begin{array}{l}{{\boldsymbol{c}}}_{{}_{i}}=\left(0,{\rm{}}0\right)c\,\rm{,}\,\,i=0\,\rm{,}\,\\ {{\boldsymbol{c}}}_{{}_{i}}=\left(\cos \frac{(i-1)\pi }{2}\,\rm{,}\,\sin \frac{(i-1)\pi }{2}\right)c,{\rm{}}i=1,{\rm{}}2,{\rm{}}3,{\rm{}}4\,\rm{,}\,\\ {{\boldsymbol{c}}}_{{}_{i}}=\sqrt{2}\left(\cos \frac{(2i-1)\pi }{4},{\rm{}}\sin \frac{(2i-1)\pi }{4}\right)c\,\rm{,}\,i=5,{\rm{}}6,{\rm{}}7,{\rm{}}8\,\rm{,}\,\end{array}\right.\end{eqnarray}$
where the subscript i = 0, …, 8 indexes the discrete velocity and c = δx/δt is the lattice speed, with the space step δx = 1 and time step δt = 1.
Figure 1. Sketch of D2Q9 discrete velocities. |
In addition, the central moments of the distribution function, equilibrium distribution function and force term satisfy the following relationships [22, 23]: 5 ) can be given as:
$\begin{eqnarray}\left|{k}_{i}\right\rangle ={\boldsymbol{T}}\left|{f}_{i}\right\rangle \,\rm{,}\,\end{eqnarray}$
$\begin{eqnarray}\left|{k}_{i}^{{\rm{e}}{\rm{q}}}\right\rangle ={\boldsymbol{T}}\left|{f}_{i}^{{\rm{e}}{\rm{q}}}\right\rangle \,\rm{,}\,\end{eqnarray}$
$\begin{eqnarray}\left|{R}_{i}\right\rangle ={\boldsymbol{T}}\left|{F}_{i}\right\rangle .\end{eqnarray}$
Here $\left|{k}_{i}\right\rangle ={\left[{k}_{0},\ldots ,{k}_{i},\ldots ,{k}_{8}\right]}^{{\rm{T}}}$ and $\left|{k}_{i}^{{\rm{e}}{\rm{q}}}\right\rangle ={\left[{k}_{0}^{{\rm{e}}{\rm{q}}},\ldots ,{k}_{i}^{{\rm{e}}{\rm{q}}},\ldots ,{k}_{8}^{{\rm{e}}{\rm{q}}}\right]}^{{\rm{T}}}$ represent the central moments of the distribution function $\left|{f}_{i}\right\rangle $ and the equilibrium distribution function $\left|{f}_{i}^{\mathrm{eq}}\right\rangle $, respectively. $\left|{R}_{i}\right\rangle ={\left[{R}_{0},\ldots ,{R}_{i},\ldots ,{R}_{8}\right]}^{{\rm{T}}}$ denotes the central moment of the discrete force term. The transformation matrix is expressed by T = [T0, T1, …, T8], where the elements take the form Ti = [ti,00, ti,10, ti,01, ti,20 + ti,02, ti,20 - ti,02, ti,11, ti,21, ti,12, ti,22]T with ${t}_{i,mn}={({c}_{ix}-{u}_{x})}^{m}{({c}_{iy}-{u}_{y})}^{n}$, with ux and uy being the horizontal and vertical components of the velocity, respectively. The elements of $\left|{k}_{i}^{\mathrm{eq}}\right\rangle $ in equation ( $\begin{eqnarray}\left\{\begin{array}{c}{k}_{0}^{{\rm{e}}{\rm{q}}}=\rho \,\rm{,}\,\\ {k}_{3}^{{\rm{e}}{\rm{q}}}=2\rho {c}_{s}^{2}\,\rm{,}\,\\ {k}_{8}^{{\rm{e}}{\rm{q}}}=\rho {c}_{s}^{4}\,\rm{,}\,\\ {k}_{i}^{{\rm{e}}{\rm{q}}}=0,{\rm{}}i={\rm{1}},{\rm{}}2,{\rm{}}4,{\rm{}}5,{\rm{}}6,{\rm{}}7\,\rm{,}\,\end{array}\right.\end{eqnarray}$
through calculation based on the Hermite polynomials, where ${c}_{s}=\sqrt{1/3}$ is the lattice sound speed.In the pseudopotential model, the interaction force that mimics molecular interactions plays a crucial role when studying phase separation. The specific expression for the intermolecular force is defined as [24]: 9 ) remains positive. ψ(x) is a pseudopotential function, the value of which is determined only by the density. $w({\left|{{\boldsymbol{c}}}_{i}\right|}^{2})$ are the weights, which are given by w(1) = 1/3 and w(2) = 1/12 in the D2Q9 model. We disregard w(0) as particles do not exert interaction forces on themselves.
$\begin{eqnarray}{{\boldsymbol{F}}}_{{\rm{int}}}({\boldsymbol{x}})=-G\psi ({\boldsymbol{x}})\displaystyle \sum _{i}w\left({\left|{{\boldsymbol{c}}}_{i}\right|}^{2}\right){\rm{}}\psi ({\boldsymbol{x}}+{\delta }_{t}{{\boldsymbol{c}}}_{i}){{\boldsymbol{c}}}_{i}\,\rm{,}\,\end{eqnarray}$
where G is a parameter that controls the intensity of the interaction force and ensures the term under the square root in equation (The choice of the pseudopotential function is pivotal in constructing the pseudopotential model. In this model, we adopt the pseudopotential function in the square root form [25]: 6 ) is written as follows [20]:
$\begin{eqnarray}\psi ({\boldsymbol{x}})=\sqrt{\displaystyle \frac{2({p}_{\mathrm{EOS}}-\rho {c}_{s}^{2})}{G{c}_{}^{2}}}\,\rm{,}\,\end{eqnarray}$
where pEOS denotes the pressure. Moreover, equation ( $\begin{eqnarray}\left\{\begin{array}{c}{R}_{1}={F}_{x}\,\rm{,}\,\\ {R}_{2}={F}_{y}\,\rm{,}\,\\ {R}_{3}=4\alpha {\left|{{\boldsymbol{F}}}_{{\rm{int}}}\right|}^{2}\,\rm{,}\,\\ {R}_{6}={F}_{x}{c}_{s}^{2}\,\rm{,}\,\\ {R}_{7}={F}_{y}{c}_{s}^{2}\,\rm{,}\,\\ {R}_{8}=\frac{4}{3}\alpha {\left|{{\boldsymbol{F}}}_{{\rm{int}}}\right|}^{2}\,\rm{,}\,\\ {R}_{i}=0,{\rm{}}i=0,{\rm{}}4,{\rm{}}5{\rm{}}.\end{array}\right.\end{eqnarray}$
Here ${\left|{{\boldsymbol{F}}}_{{\rm{int}}}\right|}^{2}={F}_{{\rm{int}},x}^{2}+{F}_{{\rm{int}},y}^{2}$. Fx and Fy represent the body force (the external force proportional to the volume) in the x and y direction, respectively. The parameter α = σ/δtψ2 is used to approximately restore the thermodynamic consistency [26].Combining equations (1 ) and (4 )-(6 ), the collision process in the central moment space can be computed as follows: 11 ) can be simplified as [20]: 2 ). Then, the density and flow velocity are updated by [20]:
$\begin{eqnarray}\begin{array}{l}\left|{k}_{i}^{* }\right\rangle =({\boldsymbol{I}}-{\boldsymbol{S}}){\boldsymbol{T}}\left|{f}_{i}\right\rangle +{\boldsymbol{ST}}\left|{f}_{i}^{\mathrm{eq}}\right\rangle +\left({\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{S}}}{2}\right){\boldsymbol{T}}\left|{F}_{i}\right\rangle \\ =({\boldsymbol{I}}-{\boldsymbol{S}})\left|{k}_{i}\right\rangle +{\boldsymbol{S}}\left|{k}_{i}^{\mathrm{eq}}\right\rangle +\left({\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{S}}}{2}\right)\left|{R}_{i}\right\rangle \,\rm{,}\,\end{array}\end{eqnarray}$
where S is a 9 × 9 diagonal matrix, its specific expression is as follows: $\begin{eqnarray}{\boldsymbol{S}}={\rm{diag}}\left[1,1,1,1,\omega ,\omega ,1,1,1\right]\,\rm{,}\,\end{eqnarray}$
where $\omega ={\left(\upsilon /{c}_{s}^{2}{\delta }_{t}+0.5\right)}^{-1}$ with ν the fluid kinematic viscosity. Using appropriate non-orthogonal basis vectors, the expression of equation ( $\begin{eqnarray}\left\{\begin{array}{c}{k}_{0}^{* }=\rho \,\rm{,}\,\\ {k}_{1}^{* }=\frac{1}{2}{F}_{x}\,\rm{,}\,\\ {k}_{2}^{* }=\frac{1}{2}{F}_{y}\,\rm{,}\,\\ {k}_{3}^{* }=\frac{2}{3}\rho +2\alpha {\left|{{\boldsymbol{F}}}_{{\rm{int}}}\right|}^{2}\,\rm{,}\,\\ {k}_{4}^{* }=(1-\omega ){k}_{4}\,\rm{,}\,\\ {k}_{5}^{* }=(1-\omega ){k}_{5}\,\rm{,}\,\\ {k}_{6}^{* }=\frac{1}{6}{F}_{y}\,\rm{,}\,\\ {k}_{7}^{* }=\frac{1}{6}{F}_{x}\,\rm{,}\,\\ {k}_{8}^{* }=\frac{1}{9}\rho +\frac{2}{3}\alpha {\left|{{\boldsymbol{F}}}_{{\rm{int}}}\right|}^{2}.\end{array}\right.\end{eqnarray}$
Through the inverse operation of the matrix, we obtain the distribution functions after the collision: $\begin{eqnarray}\left|{f}_{i}^{* }\right\rangle ={{\boldsymbol{T}}}^{-1}\left|{k}_{i}^{* }\right\rangle \,\rm{,}\,\end{eqnarray}$
which are streamed through equation ( $\begin{eqnarray}\begin{array}{rcl}\rho & = & \displaystyle \sum _{i}{f}_{i}\,\rm{,}\,\\ \rho {\boldsymbol{u}} & = & \displaystyle \sum _{i}{f}_{i}{{\boldsymbol{c}}}_{i}+\frac{{\boldsymbol{F}}}{2}{\delta }_{t}\,\rm{,}\,\end{array}\end{eqnarray}$
where u = [ux, uy] represents the flow velocity. Thus, the process of the LBE with the forcing term for a single time step is completed.In this model, we employ the Peng-Robinson EOS [27],
$\begin{eqnarray}{p}_{\mathrm{EOS}}=\displaystyle \frac{\rho RT}{1-b\rho }-\displaystyle \frac{a\varphi (T){\rho }^{2}}{1+2b\rho -{b}^{2}{\rho }^{2}}\,\rm{,}\,\end{eqnarray}$
where T represents the temperature and φ(T) = [1 + (0.37464 + 1.54226ω - 0.26992${\omega }^{2})\left(1-\sqrt{T/{T}_{c}}\right){]}^{2}$. The parameter ω is known as the acentric factor, which value is related to the type of fluid. For water, the acentric factor is chosen as ω = 0.344. The critical values in the Peng-Robinson equation can be calculated as $a=0.4572{R}^{2}{T}_{c}^{2}/{p}_{c}$ and b = 0.0778RTc/pc. In this paper, the parameters are set as a = 3/49, b = 2/21, and R = 1. Correspondingly, we can determine the critical density ρc = 2.657, temperature Tc = 0.1094, and pressure pc = 0.08936.3. Numerical verification
In this section, we perform numerical simulations to evaluate the validity of the above-outlined CLBM for both single-phase flow and multiphase flow.
3.1. Poiseuille flow
Poiseuille flow refers to the laminar flow in an infinitely long, straight circular pipe. To simulate Poiseuille flow, the initial density of the fluid is set to ρ(x, t) = ρc within a computational domain Lx × Ly = 256 × 256. The Reynolds number is set as ${\rm{Re}}={u}_{\max }{L}_{y}/\upsilon \,=\,100$, where ν refers to the fluid kinematic viscosity. The periodic condition is used in the x direction and the no-slip boundary condition is applied in the y direction. Moreover, an external force Fx is exerted along the horizontal axis. This external force drives the fluid to achieve the anticipated flow condition. The expression of the force can be written as:
$\begin{eqnarray}{F}_{x}=\frac{8\upsilon {u}_{\max }}{{L}_{y}^{2}}\,\rm{,}\,\end{eqnarray}$
where ${u}_{\max }$ signifies peak velocity when reaching a steady state. Theoretically, in the steady state, the velocity distribution of the fluid across the pipe's cross-section forms a paraboloid as follows [20]: $\begin{eqnarray}{u}_{x}=\frac{4y}{{L}_{y}}\left(1-\frac{y}{{L}_{y}}\right){u}_{\max }.\end{eqnarray}$
Figure 2 illustrates the profile of horizontal velocity along the y axis. The solid line represents the exact solution and the circles indicate the simulation results. As can be seen, the two sets of results overlap, indicating a high degree of accuracy in the simulation. Consequently, it is validated that the LBE with the forcing term is correct.
Figure 2. Velocity profile of poiseuille flow. |
3.2. Freefall
To further verify the effect of the force term, in this subsection, the freefall of a fluid in a gravitational field is simulated, and its velocity variation over time is analyzed. The fluid is initially at rest with a density of ρ = 1. Periodic boundary conditions are employed in both the vertical and horizontal directions. The fluid is subjected to a gravity with an acceleration of a.
Figure 3 depicts the variations of fluid velocity over time under different accelerations. The circles, triangles and squares denote the numerical results for a = 0.005, 0.01 and 0.02, respectively. The lines represent the corresponding theoretical solutions uy = at. It can be observed that the results obtained using the present scheme closely match the analytical predictions.
Figure 3. Evolution of the vertical velocity of a fluid in the freefall process. |
3.3. Shear flow
Shear flow refers to fluid motion induced by the relative movement of parallel boundaries without a pressure gradient. In simple shear flow situations, the streamlines are parallel. The velocity along each streamline remains constant, while a distinct velocity gradient exists perpendicular to the streamlines.
The shear effect in single-phase flow is studied in a domain with size Ly × Ly = 50 × 50. The initial density is set to ρc. When the system reaches a steady state, the theoretical solution for the amplitude of flow velocity is as follows:
$\begin{eqnarray}u(y)={u}_{y,\max }\frac{y}{h}\,\rm{,}\,\end{eqnarray}$
where ${u}_{y,\max }$ represents the maximum velocity in the vertical direction. Symbol h denotes the distance from the wall to the center of the fluid, and y represents the vertical distance to the center.Figure 4 displays the horizontal flow velocity along the y direction. The solid line represents the theoretical solution and circles indicate the numerical results. Clearly, the simulation results agree well with the theoretical data, indicating that the model can effectively simulate the shear flow.
Figure 4. Numerical and theoretical solutions for shear flow velocity. |
3.4. The gas-liquid coexistence curve
To verify the thermodynamic consistency of the current CLBM scheme, the evolution of a one-dimensional gas-liquid coexistence system is simulated. In the computational domain Lx = 500, Ly = 1, the kinematic viscosity is set to ν = 0.1. Initially, the fluid is at rest dividing into two parts: a liquid phase within the region 2Lx/5 < x < 4Lx/5 and a gaseous phase in the other part. Furthermore, the periodic boundary condition is applied to each boundary.
After the system reaches a steady state, a comparative analysis was conducted between the numerical simulation results and theoretical solutions derived from Maxwell's equal area rule. Figure 5 illustrates different outcomes when selecting different σ in equation (10 ). Circles, triangles and crosses represent LBM results of σ = 0, σ = 0.25, and σ = 0.45 respectively. The solid line denotes the analytical solution. As shown in figure 5, the introduction of the parameter σ effectively adjusts the thermodynamic consistency of the system. Furthermore, within the adjustable range, increasing σ enhances the model's capacity to tune thermodynamic consistency. For subsequent simulations, we adopt σ = 0.45.
Figure 5. Temperature versus density in the steady gas-liquid system. Symbols denote simulation results for different σ and the solid line represent the analytical solution. |
4. Numerical experiment
In this section, we employ the LBM to delve into the phase separation of gas and liquid under shear. The effects of shear velocity, temperature, viscosity coefficient, and physical domain size on the phase separation are investigated in detail.
4.1. Initial configuration
The two-dimensional phase separation process under shear boundary conditions is studied. As shown in figure 6, the upper plate moves rightward at a speed of u0, and the lower plate moves leftward at the same speed. Between the two plates is the fluid with initial density ρ = ρc(1 ± RAND), where RAND represents a random function within the range (-0.1, 0.1), introducing minor fluctuations that promote phase separation.
Figure 6. Initial density distribution and boundary conditions of shear flow. |
In this manuscript, the random function for the fluctuation in density is kept identical for all cases. The periodic boundary condition is adopted in the x direction, and the shear boundary condition is imposed on the plates [14]. This configuration is relevant to certain real physical scenarios, such as gas-liquid separation processes in nuclear power plant circuits. The CLBM can also be adapted with other boundary conditions and computational domains to represent corresponding physical processes.
4.2. Case I: droplet
To begin with, let us consider a physical domain of Lx × Ly = 256 × 256, with the shear velocity u0 = 0.1, temperature T = 0.82Tc, and viscosity ν = 0.1. Figure 7 illustrates the snapshots of the density field during the evolution of phase separation under shear. In the early stage, the density difference between the liquid (high density) and the gas (low density) is minimal at t = 20, and it gradually increases as the phase transition progresses, reaching a more pronounced difference at t = 80. As time advances, adjacent droplets begin to merge and grow larger, with the gaps between them also expanding, as shown from t = 500 to 5000. Under the influence of surface tension and viscous shear, the liquid transitions from a complex shape to an ellipsoidal form, as seen at t = 10,000 and 30,000.
Figure 7. Density contours in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.82Tc, and ν = 0.1. |
Additionally, figure 8 depicts the rate of change in density at the same time instants as in figure 7. Blue regions indicate negative values, signifying a decrease in density over time, while red regions represent positive values, indicating an increase. A comparison between figures 7 and 8 shows that density changes occur most rapidly at the interface between the gas and liquid, while remaining nearly constant within the gas or liquid phases. In the early stage at t = 20, the absolute values of both the maximum and minimum change rates are relatively small, indicating a slow phase transition. The transition speed then accelerates, reaching its peak around t = 80. Afterwards, the absolute values of the change rate decrease, signifying a deceleration in the phase transition. As time progresses, both the maximum and minimum change rates approach zero, reflecting the system's approach to a steady state.
Figure 8. Change rate of density in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.82Tc, and ν = 0.1. |
Furthermore, figure 9 depicts the evolution of average density gradients corresponding to figure 7. The solid, dashed, and dotted lines represent the gradients $\overline{{\rm{\nabla }}\rho }$, $\overline{{\rm{\nabla }}{\rho }_{x}}$, and $\overline{{\rm{\nabla }}{\rho }_{y}}$, respectively. From figure 9, we can see two distinct stages, i.e. spinodal decomposition and domain growth. In the spinodal decomposition stage, from t = 0 to t = 80, gas-liquid separation lead to rapid growth of the interfaces, which result in increase of the density gradient simultaneously. It can be noticed that densities reach their peaks at about t = 80, marking the transition point between the two stages. During the domain growth stage beginning at t = 80, small droplets and bubbles coalesce to minimize free energy, forming larger droplets and bubbles. The summation of all gas-liquid interfaces gradually decreases, leading to a corresponding reduction in the overall density gradient.
Figure 9. Profiles of average density gradients in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.82Tc, and ν = 0.1. |
4.3. Case II: band affected by temperature
Next, let us investigate the influence of the temperature on the process of phase separation under shear. Numerical experiments showed that temperature significantly affects the final pattern of the liquid. The simulation is under the conditions of Lx × Ly = 256 × 256, u0 = 0.1, and ν = 0.1. As shown in table 1, the tests demonstrate that a band structure emerges when the temperature satisfies T ≥ 0.824Tc, whereas a droplet formation is observed when T ≤ 0.823Tc. This indicates that the temperature T = 0.8235Tc serves as a critical threshold in this process.
Table 1. Final patterns for various temperatures. |
| Temperature | 0.81Tc | 0.82Tc | 0.823Tc | 0.824Tc | 0.825Tc | 0.85Tc | 0.9Tc |
|---|---|---|---|---|---|---|---|
| Final pattern | Droplet | Droplet | Droplet | Band | Band | Band | Band |
Figure 10 displays the density distribution for the case of T = 0.9Tc in table 1. Compared to the temperature T = 0.82Tc in figure 7, the temperature is increased to T = 0.9Tc in figure 10. Both similarities and differences can be observed between the two figures. The overall pattern of density variation remains generally consistent. Initially, numerous small liquid-gas interfaces appear, as illustrated at t = 140. Over time, nearby droplets start to coalesce and increase, while the spaces between them widen progressively, as shown between t = 500 and t = 5000. However, in the final equilibrium state, the liquid distinctly forms a band-like structure. This outcome arises from two competing mechanisms: the shear effect, which stretches the fluid and induces internal deformation leading to the formation of bands, and the surface tension at the gas-liquid interface, which favors the development of droplets. As the temperature increases, the repulsive force between molecules is enhanced and the attractive force is weakened, causing a reduction in surface tension and phase separation. Consequently, the liquid forms a band-like structure.
Figure 10. Density contours in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.9Tc, and ν = 0.1. |
Moreover, figure 11 illustrates the temporal change rate of density under the same conditions as in figure 10. In comparison with figure 8 where the temperature is lower, the changing process of the absolute value follows a similar trajectory, beginning close to zero, followed by a rapid increase and then a gradual decrease. The distinction lies in the timing, as the change rate of density in figure 8 is relatively high at t = 80, whereas the one in figure 11 achieves this around t = 140. In fact, this finding is further confirmed by the density gradient in figure 12.
Figure 11. Changing rate of density in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.9Tc, and ν = 0.1. |
Figure 12. Profiles of average density gradients in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.9Tc, and ν = 0.1. |
Figure 12 displays the profile of average density gradients at T = 0.9Tc. Similar to figure 9, figure 12 shows two distinct stages. During the first stage, the density gradient increases rapidly. In the second stage, it gradually decreases. However, the density gradient reaches its peak at the time instant t = 144 in figure 12, compared to t = 80 in figure 9, indicating an extended spinodal decomposition phase. With an increasing temperature, the phase separation slows, hence the cumulative effect of shear increases. Thus the liquid is more likely to form bands as the temperature rises, which aligns with theoretical predictions. Additionally, a comparison between figures 9 and 12 indicates that the maximum density gradient becomes lower for T = 0.9Tc than for T = 0.82Tc. This difference arises due to the increase in temperature, which reduces the density difference between the two phases, thereby lowering the peak density gradient.
In addition, it should be emphasized that this critical temperature is subject to change with variations in other factors such as flow velocity, viscosity, and domain size, which will be explored in subsequent parts.
In fact, from the point of view of two competing mechanisms, e.g., with greater velocity or viscosity, the shear force will be more dominant, and it will be easier to form a band. Numerical tests show that the critical temperatures are 0.8235Tc and 0.8215Tc for u0 = 0.1 and u0 = 0.104, respectively. And the critical temperatures are 0.8245Tc and 0.8235Tc for ν = 0.09 and ν = 0.1, respectively. Therefore, the critical temperature decreases with larger velocity or viscosity.
4.4. Case III: band affected by velocity
Next, we examine the influence of shear velocity on the phase separation process. Numerical experiments showed that velocity has a significant effect on the final pattern of the liquid. The test results are presented in table 2. The simulation is under the conditions of Lx × Ly = 256 × 256, T = 0.82Tc, and ν = 0.1. Here the shear velocity varies while other parameters remain consistent with those in figure 7. We can see that a band structure forms for u0 ≥ 0.105, while droplet formation occurs when u0 ≤ 0.104. This suggests that u0 = 0.1045 acts as a critical threshold in this process. In fact, the critical threshold of flow velocity varies depending on other physical factors, including the temperature, viscosity, and domain area.
Table 2. Final patterns for various velocities. |
| Shear velocity | 0.09 | 0.1 | 0.104 | 0.105 | 0.106 | 0.13 | 0.15 |
|---|---|---|---|---|---|---|---|
| Final pattern | Droplet | Droplet | Droplet | Band | Band | Band | Band |
Figure 13 illustrates the snapshots of the density field under the conditions of u0 = 0.15 in table 2. Compared to figure 7, where the flow velocity is lower, figure 13 depicts the scenario with an elevated flow velocity of u0 = 0.15. A phase transition is observed in figure 13, akin to the process depicted in figure 7, albeit with distinct differences in the pattern of the physical field. Initially, several small liquid-gas interfaces emerge, as shown at t = 80. As the system evolves, the smaller droplets gradually vanish, while the larger droplets continue to grow, as shown between t = 500 and t = 3000. Subsequently, the liquid begins to form a band-like structure between t = 3000 and t = 30,000, which is different from the droplet shape in figure 7. This phenomenon is attributed to the increased shear velocity, which intensifies the shear effect, stretching the fluid and inducing internal deformation, ultimately leading to the formation of bands.
Figure 13. Density contours in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.15, T = 0.82Tc, and ν = 0.1. |
Moreover, figure 14 illustrates the temporal change rate of density at u0 = 0.15. In comparison with figure 8, where the shear velocity is lower, the absolute value of the temporal change rate of density follows a similar trajectory, starting near zero, rising rapidly, and then gradually decreasing. In both figures, the average value is higher at t = 80. The difference is that the absolute value of the temporal change rate is relatively higher in figure 14 over time. This is due to the increased shear velocity, which accelerates the density changes, resulting in higher temporal change rate.
Figure 14. Change rate of density in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.15, T = 0.82Tc, and ν = 0.1. |
Figure 15 presents the corresponding profile of average density gradients at u0 = 0.15. Similar to figure 9, figure 15 exhibits two different stages. The density gradient rises rapidly in the early stage and declines gradually in the later stage. The density gradient reaches its peak at t = 80 in both figures, indicating that the duration of the spinodal decomposition stage is identical. The maximum gradient values are nearly the same in both cases. However, the density gradient declines more rapidly in figure 15, which is attributed to the increased shear velocity that accelerates the phase separation process.
Figure 15. Profiles of average density gradients in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.15, T = 0.82Tc, and ν = 0.1. |
4.5. Case IV: band affected by viscosity
Now, we explore the impact of the viscosity upon the phase separation under shear. Numerical experiments indicates that viscosity greatly influences the final pattern of the liquid. Table 3 displays the test results. The simulation is under the conditions of T = 0.82Tc, Lx × Ly = 256 × 256, and u0 = 0.1. Here the viscosity coefficient is adjusted while maintaining all other parameters the same as in figure 7. Clearly, a band structure forms at viscosity of ν ≥ 0.12, whereas droplet formation is observed when ν ≤ 0.11. This implies that ν = 0.115 functions as a crucial threshold. Similarly, the crucial threshold of the viscosity coefficient depends upon the temperature, shear velocity and domain size as well.
Table 3. Final patterns for various viscosity coefficients. |
| Viscosity coefficient | 0.05 | 0.1 | 0.11 | 0.12 | 0.14 | 0.2 | 0.3 |
|---|---|---|---|---|---|---|---|
| Final pattern | Droplet | Droplet | Droplet | Band | Band | Band | Band |
Figure 16 displays the density distribution under the conditions of ν = 0.3 in table 3. Compared to figure 7, the viscosity coefficient is increased to ν = 0.3 in figure 16. In the early stage (from t = 0 to 100), multiple interfaces between the two phases emerge, similar to the behavior observed at ν = 0.1. As time goes on, droplets exhibit a tendency to absorb surrounding molecules, resulting in fewer but larger droplets (from t = 500 to t = 5000). The difference is that the liquid tends to form a droplet in figure 7 and a band-like structure in figure 16. This demonstrates that the increasing the viscosity coefficient promotes band formation by shearing liquid.
Figure 16. Density contours in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.82Tc, and ν = 0.3. |
Additionally, figure 17 illustrates the rate of change in density corresponding to figure 16. A comparison between figures 17 and 8 reveals both similarities and differences. The absolute value of the density's temporal change rate follows a similar rise-and-fall pattern, but with a timing difference: the average changing rate peaks at t = 80 in figure 8, whereas it occurs at t = 100 in figure 17. This observation is further supported by the density gradient shown in figure 18.
Figure 17. Change rate of density in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.82Tc, and ν = 0.3. |
Figure 18. Profiles of average density gradients in the evolution of phase separation under shear: Lx × Ly = 256 × 256, u0 = 0.1, T = 0.82Tc, and ν = 0.3. |
Figure 18 displays the density gradient under the same conditions. Compared to figure 9, figure 18 displays two stages: an initial rapid increase in the density gradient followed by a gradual decline. However, the peak is delayed, occurring at t = 98 in figure 18 compared to t = 80 in figure 9, indicating an extended spinodal decomposition stage. Moreover, the inserts suggest that the peak value decreases with a higher viscosity coefficient. As viscosity increases, the evolution process slows, lengthening the time of phase separation and resulting in a smaller density gradient.
4.6. Case V : band affected by domain size
Finally, let us probe the influence of the physical domain size on the dynamic process of phase separation under shear. Numerical experiments demonstrated that domain size plays a substantial role in determining the final pattern of the liquid. The test results are shown in table 4. The simulation is under the conditions of Lx × Ly = 128 × 128. It is important to highlight that under different domain sizes, with all other parameters held constant as in figure 7, a band structure appears when the domain length Lx = Ly ≤ 230, while droplet formation occurs for the domain length when the domain length Lx = Ly ≥ 240. This suggests that Lx × Ly = 235 serves as a key threshold. Moreover, it is important to note that the critical threshold of domain size is also influenced by temperature, shear velocity, and viscosity.
Table 4. Final patterns for various domain sizes. |
| Domain length | 128 | 200 | 230 | 240 | 250 | 256 | 280 |
|---|---|---|---|---|---|---|---|
| Final pattern | Band | Band | Band | Droplet | Droplet | Droplet | Droplet |
It should be mentioned that the initial configuration of density field is set up using the same random function. The density variation pattern varies when the physical domain changes. However, from the perspective of statistic physics, the physical laws and characteristics are identical, although the initial density variation pattern is different. For the sake of verification, we have conduct another simulation with a density variation pattern different from figure 19, and find that the final pattern is also a band, which is similar to figure 19. Consequently, the conclusion remains the same.
Figure 19. Density contours in the evolution of phase separation under shear: Lx × Ly = 128 × 128, u0 = 0.1, T = 0.9Tc, and ν = 0.1. |
Figure 19 depicts the density contours under the conditions of ν = 0.1, T = 0.82Tc, Lx × Ly = 128 × 128, and u0 = 0.1. Although the physical domain is smaller than in figure 7, the overall density variation pattern remains consistent. For example, mini liquid-gas interfaces initially form at about t = 80. As the system progresses, smaller droplets gradually disappear, while larger ones continue to expand, as observed between t = 500 and t = 5000. However, significant differences emerge in the steady state, where the liquid distinctly forms a band-like structure. In smaller domains, with the shear velocity held constant, the confined liquid molecules experience intensified shear effects, resulting in a larger velocity gradient. This gradient promotes the formation of a band-like structure within the liquid.
Moreover, figure 20 illustrates the rate of density change in the evolution for Lx × Ly = 128 × 128. Compared to figure 8, which features a larger computational size, the absolute value follows a similar trajectory, with a rapid rise followed by a gradual decline. However, the maximum absolute value of the density change rate is smaller in figure 20, due to the less gas-liquid interface with the reduced computational domain.
Figure 20. Change rate of density in the evolution of phase separation under shear: Lx × Ly = 128 × 128, u0 = 0.1, T = 0.82Tc, and ν = 0.1. |
Figure 21 further depicts the density gradient corresponding to figure 20. The gradient rises rapidly initially and then follows a gradual decrease, which is comparable to figure 9. The density gradient reaches its maximum at t = 81, marking the transition between the two stages. The timing of the peak is nearly identical to that in figure 9, suggesting a comparable spinodal decomposition phase. It is also interesting to find that the maxima of the average density gradient are similar in figures 21 and 9. This indicates a good spatial isotropy in the phase separation under consideration.
Figure 21. Profiles of average density gradients in the evolution of phase separation under shear: Lx × Ly = 128 × 128, u0 = 0.1, T = 0.82Tc, and ν = 0.1. |
5. Conclusion
In this paper, we utilize a CLBM incorporating the pseudopotential model and the Peng-Robinson EOS to simulate liquid-vapor phase separation under viscous shear. The model is first validated through numerical simulations of Poiseuille flow, free fall, shear flow, and assessments of thermodynamic consistency. Subsequently, we conduct an in-depth investigation into the influence of temperature, shear velocity, viscosity, and domain size on phase separation under shear. The system's evolution reveals both commonalities and distinctions across various conditions, which are explored in detail to elucidate the underlying physical mechanisms.
On the one hand, the general pattern of phase separation remains consistent: initially, numerous small liquid-gas interfaces form. As the system progresses, smaller droplets gradually disappear while larger droplets continue to grow. Additionally, the temporal change rate of density follows a similar trajectory, beginning near zero, rising rapidly, and then gradually decreasing. Moreover, the average density gradient increases rapidly during the spinodal decomposition stage and decreases gradually during the domain growth stage. The maximum of the average density gradient is used as a physical criterion to distinguish the two stages.
On the other hand, altering parameters leads to distinct evolution of the multiphase shear flow, especially in the final state. The different outcomes result from two competing mechanisms: surface tension and shear force. Specifically, the liquid tends to form a droplet when the surface tension dominates; the liquid tends to form a band when shear force dominates. In fact, surface tension is weakened by increasing temperature, while viscous shear is amplified by increasing shear velocity and viscosity or by reducing the physical domain size. Consequently, during the evolution of phase separation under shear, the liquid tends to form droplets under conditions of low temperature, shear velocity, and viscosity, and in larger domain sizes. Otherwise, it tends to form bands.