In this work, the Rothman- Keller (RK) type LB model (Leclaire, S., et al., 2017; Xie et al., 2018b) is selected as the basic multiphase flow solver, and the D2Q9 model (Qian, Y.H., et al., 1992) is chosen to discrete the velocity space, which is defined as: c i = [ c i x , c i y ] = { [ 0,0 ] , i = 0 { cos [ π 2 ( i − 1 ) ] , sin [ π 2 ( i − 1 ) ] } c , i = 1,2,3,4 { cos [ π 4 + π 2 ( i − 5 ) ] , sin [ π 4 + π 2 ( i − 5 ) ] } 2 c , i = 5,6,7,8 (1) with the weight parameters W i being: W i = { 4 / 9 , i = 0 1 / 9 , i = 1,2,3,4 1 / 36 , i = 5,6,7,8 (2)

In the model, each lattice Boltzmann equation for fluid k is defined as: f i k ( x + c i ∆ t , t + ∆ t ) − f i k ( x , t ) = Ω i k [ f i k ( x , t ) ] + Γ i ( ρ k ρ F ( x , t ) ) (3) where f i k is the distribution function, the subscript i denotes the direction of the discrete velocity space; x represents the lattice position, and Δ t is the time step; Ω i k and Γ i are the collision term and the modified forcing term (Guo, Z., et al., 2002), respectively.

The fluid density and the flow velocity are calculated as follows: ρ k = ∑ i f i k = ∑ i f i k , e q (4) ρ = ∑ k ρ k , u = ( F / 2 + ∑ i ∑ k f i k c i ) / ρ (5) where ρ k represents the macroscopic phase density, ρ is the total density, and u denotes the flow velocity.

The equilibrium distribution function f i k , e q can be given by f i k , e q ( ρ k ,   u , α k ) = ρ k ( ϕ i k + W i ( 3 c i ⋅ u c 2 + 9 ( c i ⋅ u ) 2 2 c 4 − 3 u 2 2 c 2 ) ) + Φ i k (6) where ϕ i k = { α k , i = 0 ( 1 − α k ) / 5 , i = 1,2,3,4 ( 1 − α k ) / 20 , i = 5,6,7,8 (7) where α k ( 0 < α k < 1 ) is a parameter related to the density ratio between fluids 1 − α l 1 − α k = ρ k , 0 ρ l , 0 (8) where superscript 0 represents the initial value of density.

To obtain the right momentum for different density ratios, the Φ i k in Eq. 6 is set as (Leclaire, S. et al., 2013a): Φ i k = { − 3 ν e f f ( u ⋅ ∇ ρ k ) / c 2 , i = 9 4 ν e f f ( D k : c i ⊗ c i ) / c 4 , i = 1,2,3,4 ν e f f ( D k : c i ⊗ c i ) / c 4 , i = 5,6,7,8 , (9) with D k = [ ( u ⊗ ∇ ρ k ) + ( u ⊗ ∇ ρ k ) T ] / 8 .

The collision term Ω i k in Eq. 3 is combined by three sub operators as: Ω i k = ( Ω i k ) M [ ( Ω i k ) S + ( Ω i k ) I ] , (10) where ( Ω i k ) S is the single-phase operator, ( Ω i k ) I is the perturbation operator and ( Ω i k ) M is the recoloring operator.

The perturbation operator related to the interfacial effect of multiphase flow is written as: ( Ω i k ) I = f i k + ∑ l ( l ≠ k ) A k l C k l 2 | G k l | ( W i ( c i ⋅ G k l ) 2 | c i | 2 | G k l | 2 − B i ) ∆ x , (11) where B i = { − 4 / 27 , i = 0 2 / 27 , i = 1,2,3,4 5 / 108 , i = 5,6,7,8 (12) and G k l is the “color” gradient perpendicular to the interface between the phases k and l : G k l ( x , t ) = ρ l ρ ∇ ( ρ k ρ ) − ρ k ρ ∇ ( ρ l ρ ) (13)

A concentration factor that limits the range of interfacial action (Leclaire, S., et al., 2013b) is defined as: C k l = min { 10 6 ρ k ρ l ρ k , 0 ρ l , 0 , 1 } (14) and A k l is related to interfacial tension γ k l : A k l = A l k = 9 2 γ k l τ e f f c 2 Δ x (15)

The recoloring operator, which ensures that each phase satisfies the conservation of mass perfectly, is written as ( Ω i k ) M = ρ k ρ ∑ k f i k + ∑ l ( l ≠ k ) [ β k l ρ k ρ l ρ 2 cos ( φ i k l ) × ∑ k f i k , e q ( ρ k , 0 , α k ) ] (16) where the parameter β k l controls the interface thickness of phases k and l , and the φ i k l is the angle between c i and G k l .

To obtain the apparent viscosity of decane-water emulsion, we also simulate the single-phase flow of water as a reference. Based on the comparison between the single- and two-phase flow, the apparent viscosity can be calculated as (Xie and Balhoff, 2021): μ a p p = μ s F t / Q t F s / Q s (17) where the subscripts t and s denote the two-phase and single-phase flow states, respectively, μ s is the dynamic viscosity of the single-phase reference fluid, F and Q represent the body force and flow rate, respectively.

In our simulations, we apply a body force of the same magnitude for both single- and two-phase flows. Therefore, the above equation to obtain the apparent viscosity is simplified as: μ a p p = μ s Q s Q t (18)

The computation domain is a homogeneous porous media with a total length of 14 mm and a width of 9.5 mm, containing evenly distributed round solid grains, as shown in Figure 4. The diameter of the solid grain is 2.8 mm. The periodic boundary condition is applied in the horizontal direction, and the fluids are driven by a horizontal body force from left to right. For the two-phase flow simulations, the fluids are water and decane, while for the single-phase flow simulation, the fluid is water. Their physical properties and other key parameters in the LBM simulation are listed in Table 1.

42 simulation cases in total are prepared to discuss the effect of initial distribution, saturation, interfacial tension (IFT), and contact angle (CA) on the apparent viscosity of the decane-water emulsion in porous media. We compare three types of initial fluid distributions, including the upstream-decane distribution, the downstream-decane distribution, and the random distribution. The saturation, interfacial tension, and contact angle are varied from 10% to 90%, 10 mN/m to 50 mN/m, and 0–180°, respectively. In the default case, the decane saturation is 20%, initially randomly distributed in the porous media with a contact angle of 45°, and the interfacial tension between the decane and water is 20 mN/m. A summary of parameters discussed in all simulation cases is given in Table 2.

For each case, the simulation is terminated once the flow stabilizes and the fluid distribution does not change. The inlet and outlet flow rates of the single-phase case and the default two-phase flow case in the last 40000 steps (the gap is 10000 steps) are shown in Table 3. In both cases, the inlet flow rates are approximately equal to the outlet flow rates, and the gaps between the inlet and outlet rates stabilize in the last 40000 steps, which illustrates that the total time step we set is adequate for the flow stabilization.

We discuss three kinds of initial fluid distributions, which include the random distribution, upstream distribution, and downstream distribution as shown in Figures 5A,C,E. The simulation results at the last steps are also shown in Figures 5B,D,F. The results show that the different initial states lead to different distributions at the last time step. However, all of the last-time distributions are qualitatively random, no matter how the fluids are distributed at the initial time step. This result indicates that the impact of initial phase distribution on the apparent viscosity is weak.

The effect of decane saturation on the apparent viscosity of the emulsion system is shown in Figure 6. When the decane saturation is 0, the apparent viscosity is the water dynamic viscosity of 1 mPa s. While the decane is added to form the decane-water emulsion, the apparent viscosity increases. This is because the original single-phase flow state becomes the two-phase flow state, and the interfacial interaction between water and decane occurs and the capillary trapping force arises, leading to the increase in the apparent viscosity, even though the dynamic viscosity of decane is lower than that of water. The apparent viscosity of the decane-water reaches its maximum value of 4.27 mPa s when the decane saturation is 20%, after which, the apparent viscosity gradually decreases with the increase in the decane saturation. The emulsion’s apparent viscosity decreases to the water viscosity of 1 mPa s when the decane saturation reaches around 95%. It is worth mentioning that this value may vary in a small range if different fitting functions are applied. If the decane saturation is higher than 95%, the emulsion’s apparent viscosity is lower than 1 mPa s, and finally reaches its minimum value, that is the dynamic viscosity of decane (0.85 mPa s).

Figure 7 shows the effect of interfacial tension between the water and decane on the emulsion’s apparent viscosity. Here we also change the decane saturation from 20% to 80%. We find that the emulsion’s apparent viscosity increases with the increase in interfacial tension for all satuations. This is because the capillary dragging force increases with interfacial tension. Therefore, the increase in interfacial tension decreases the effective flow rates, which leads to higher apparent viscosities of the emulsion. The maximum value of apparent viscosity is 22.3 mPa s, which is about 20 times greater than the dynamic viscosity of water. The minimum of that is 1.37 mPa s, also greater than the dynamic viscosity of water.

Figure 8 shows the change in the emulsion’s apparent viscosity with contact angles. Here we also consider the change in the decane saturation from 20% to 80%. As is seen, for all satuations, higher apparent viscosities are found in decane-wet and water-wet systems, and the apparent viscosity reaches the lowest value in intermediate-wet porous media. This is because the decane-wet and water-wet porous geometry have stronger adsorption capacities for decane and water, respectively, than that of intermediate-wet media, and the stronger adsorption capacity leads to the higher apparent viscosity.