{'fraction_non_alphanumeric': 0.03497496450190569, 'fraction_numerical': 0.011845153575965921, 'mean_word_length': 4.668926074984114, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Drag coefficient of a rising bubble in a shear-thinning fluid using the power-law scheme coupled with a Cahn-Hilliard equation with a variable mobility: A lattice Boltzmann study and comparison with experiment

This study aims to investigate the behavior of multicomponent fluid flows consisting of Newtonian and non-Newtonian components, especially terminal velocity of a rising bubble in a power-law fluid. A recent lattice Boltzmann (LB) model is extended using power-law scheme to be able to simulate both Newtonian and non-Newtonian fluid flows at high density and viscosity ratios. Also, a variable mobility is introduced in this study to minimize the unphysical error around small bubbles in the domain. A three-component fluid flow system is examined using a constant and variable mobility. It is shown that each component has more stability using variable mobility while constant mobility causes interface dissipation, leading to mass loss gradually. In addition, two test cases including power-law fluid flows driven between two parallel plates are conducted to show the accuracy and capability of the model. To find a grid-independent computational domain, a grid independency test is carried out to show that a 200*400 domain size is suitable for our computations. Then, terminal velocity of a rising bubble is compared to an existing correlation in the literature, indicating that the results are in good agreement with existing study so that average relative error in six different cases is 5.66 %. Also, the simulated examples show good conformity to experimental results over a range of the Reynolds and Eotvos numbers.


Introduction
Multiphase flows, especially rising bubble in Newtonian and non-Newtonian liquids, are widely applicable in the industry. Multiphase flow systems exist in many contexts such as water treatment processes, biomedical engineering, pharmaceutical and food industries, and so on [1,2]. Also, multiphase flow systems are found in the droplet-based microfluidic devices which is applicable in drug delivery and cell analysis [3,4,5]. Therefore, many research efforts are targeted on predicting interaction between liquid and gas not only in Newtonian liquids but also in non-Newtonian ones. Since gas phase is in the form of bubble in many systems, it is of importance to find the behaviors and characteristics of the bubble in Newtonian and non-Newtonian liquids. As a matter of fact, understanding two-phase flow behavior can help improve and maintain such systems. Terminal velocity of a bubble is an important parameter in two-phase flow systems which plays an important role in optimization and designing of gas-liquid systems [6], showing the necessity of investigating bubble dynamics in Newtonian and non-Newtonian liquids. In addition, terminal velocity can exhibit equilibrium in a system, as it occurs when drag and buoyant force are balanced.
Numerical analysis has become very popular because of its prediction and comprehension of phenomena in nature. Many numerical approaches have been applied to analyze multiphase flows in macroscopic scale like front-tracking [7], volume of fluid [8], level set method [9], etc. Battistella et al. [6] used a fronttracking model to simulate rising bubble in a power-law fluid. They used a drag closure proposed by dijkhuizen et al. [10] to obtain terminal velocity of bubble, and compare it with numerical simulation results. In most cases, the accuracy of the model was 20%. Li et al. [11] experimentally investigated rising bubble in a non-Newtonian liquid. They proposed a drag coefficient correlation by considering the Reynolds number, aspect ratio, and rheological properties of the liquid. Also, they plotted a map of shape regimes for bubbles in a non-Newtonian fluid using Reynolds, Eötvös, and Morton numbers. Islam et al. [12] used volumeof-fluid (VOF) approach to investigate rising bubble behavior in a Newtonian and non-Newtonian liquid. Behavior of a non-Newtonian liquid was described by power-law equation. They studied terminal velocity and deformation shape of bubble using numerical simulations. They showed that their results are in good agreement with empirical correlations. Ohta et al. [13] used coupled level-set and volume-of-fluid (CLSVOF) to simulate bubble rising in a non-Newtonian liquid.
The viscosity of liquid was determined by Carreau model. They investigated bubble rise velocity and bubble shape in non-Newtonian fluids by changing Carreau model parameters. They concluded that in case of implementing effective parameters, Bhaga-Weber map can predict bubble shape in the non-Newtonian fluid.
Premlata et al. [14] carried out VOF method using the Carreau-Yasuda model to simulate rising bubble in a non-Newtonian fluid. They validated their model by comparing with results in the literature. They showed that rising velocity increases by increasing shear-thinning tendency and deformation of bubble increases at high Gallilei number and low Eötvös number. Pang and Lu [15] used VOF method to simulate bubble rising in Carreau fluid. They investigated the influence of surface tension, rheological index, and different gravity field on bubble dynamics. They found that bubble deformation increases by increasing Gallilei and Eötvös numbers and decreases with the increase of rheological index. Also, they showed that terminal velocity increases by increasing Gallilei number and decreasing rheological index.
In the recent decades, lattice Boltzmann method (LBM) has been used as an alternative way to simulate fluids dynamics in mesoscopic scale. The LBM has several advantages over traditional computational fluid dynamics (CFD) methods such as ease of implementation, capability of scaling on parallel computers, and handling complex geometries [16,17,18]. Several LB models have been proposed for two-phase fluids such as color gradient [19], pseudopotential [20], free energy [21], and mean-field [22]. All of these early models confront with some drawbacks like limited density ratio and spurious currents (non-physical velocity fields). Inamuro et al. [23] proposed a model based on free energy that can capture large density ratio, but the Poisson equation in the model reduced simplicity and efficiency of the LBM. Lee [24] proposed a two-phase model using  [27]. Beside above models, many efforts have been done to capture non-Newtonian behavior of fluids in LBM. Gabbanelli et al. [28] extended LBM to non-Newtonian formulation using a truncated power-law model to capture non-Newtonian characteristics of singlecomponent fluids. Boyd et al. [29] proposed a second order LB model to simulate non-Newtonian fluids. Their model is efficient computationally and more accurate compared to previous models. Moreover, further works were done for multiphase non-Newtonian fluids [30,31,32].
To the authors' best knowledge, the model introduced in [25] is strong due to many numerical tests conducted by the authors. Therefore, in this study, this model is used for the computation of multiphase flows. In this paper, we enhance the model using variable mobility to minimize little numerical dissipation around small bubbles reported in [27]. Also, Variable and constant mobilities are compared to show how mass conservation is satisfied in the present study. Variable mobility can effectively reduce dissolutions of bubble and droplets to surrounding fluid [26,33,34]. Also, we couple the mentioned model to the power-law scheme, using the truncated power-law model for viscosity distribution, to simulate non-Newtonian fluids at high density and viscosity ratios.
The paper is organized as follows. Section 2 represents lattice Boltzmann equations (LBEs) together with CH equations and explanations about the truncated power-law model. Benefits of variable mobility and validation of method, which includes two non-Newtonian test cases studied by the present method and compared with analytical solution to show accuracy and robustness of the extension, are carried out in section 3. Then in section 4, grid convergence test is carried out and rising behavior of a bubble in a shear-thinning fluid is examined. Also, terminal velocity and deformation shape of the bubble are compared with empirical results. Finally, this work is concluded in section 5.

Interface capturing model
Consider a domain Ω of three incompressible, immiscible fluids. Each fluid is distinguished by its phase-field variable (C i i = 1 − 3) along with the following constraint: Free energy for three-component flows is as follows [26]: where ξ is the interface thickness which is usually of 3 to 5 lattices, Σ i are capillary coefficients, and E 0 is the bulk free energy. For consistency with diphasic systems, the capillary coefficients are defined as Σ i = σ i j + σ ik − σ jk (i, j, k = 1 − 3) where σ i j is the surface tension coefficient between fluid i and fluid j. The bulk free energy is defined as where Λ is a positive parameter for total spreading scenarios.
The three-component CH equations can be written as following: where t is the time, u is the velocity, and µ i is the chemical potential of the fluid i.
the mobility of the fluid i, D i , and Σ T are defined as where D 0 = ξ/300 is an auxiliary parameter [25].

Lattice Boltzmann equations
The LBE with the single relaxation time can be written as [24] f α (x + e α δt, where f α is the particle distribution function in the α-direction, λ is the relaxation time, ρ is the density of mixture, and F is the intermolecular force modeling nonideal gas effects. f eq α is the equilibrium distribution function which are defined as follows: In this study, the D2Q9 lattice structure is used with the following weight coefficients, w α , and microscopic velocity, e α , set: In the above equations, c s = c √ 3 is the speed of sound in the D2Q9 lattice in which c = δx/δt is the lattice speed with δx and δt being the length scale and time scale, respectively. The force term F is expressed as where F b = −ρg − → j is the body force with g being the gravitational acceleration.
To obtain a LBE for recovering the hydrodynamics properties of the mixture, the following transformation is used: where Γ α (u) = f eq α /ρ. Taking the total derivative of g α and using the trapezoidal rule results in the following LBE: where the dimensionless relaxation time τ = λ/δt is related to kinematic viscosity ν by τ = ν/c 2 s . In the above equations, whereḡ α andḡ eq α are modified particle and equilibrium distribution functions, respectively, which are defined as: In order to track the interfaces among three fluids, we use the following distribution functions Note that C 3 is obtained via the constraint in Eq. (1).
Like previously, by taking the total derivative of h i α and using the trapezoidal rule, we have the following LBE for interfaces capturing: where τ h = 0.5 is the dimensionless relaxation time.h i α andh i,eq α are the modified particle and equilibrium distribution functions, respectively, defined as where h i,eq As mentioned earlier, in this study we employ a variable mobility as it shows better mass conservation properties. Equation (6) defines both constant and variable mobilities. As can be seen, the variable mobility is not constant and dependent of C i . As such, when taking the total derivative of Eq. (17) and considering the variable mobilities, the diffusion term on the right-hand side of Eq. (4) simplifies to which the last term is related to the variable mobilities and also appears in the interface capturing LBE, i.e., Eq. (18).
In the above equations, second-order central and mixed finite differences and directional derivatives of a scaler variable are implemented to enhance the numerical stability [27].
The phase-field variables, momentum and pressure can be calculated by taking the moments of the modified distribution functions

Strain rate tensor and relaxation time calculation
In the continuum fluid flows, the Navier-stokes (NS) equations are recovered from LBE through the Chapman-Enskog (CE) analysis. The strain rate tensor is calculated by the second moment of non-equilibrium distribution function as follows [35,36]: where u t denotes the transposed matrix to u and Π (1) as first-order momentum flux for the D2Q9 lattice is given as where f (1) i is a first-order perturbation distribution function [35]. Using the transformation in Eq. (13), the relationship between non-equilibrium particle distribu-tion functions can be obtained as follows: By using Eqs. (15) and (16) we have As such, the non-equilibrium particle distribution function of f can be replaced as follows: As a result, knowing the magnitude of non-equilibrium distribution function of f , Constitutive equation for a generalized Newtonian fluid is defined as [37] T = 2η(γ)D (30) where η(γ) and T are the apparent viscosity function of shear rate and shear stress tensor, respectively. The non-Newtonian fluid in this study is modeled as a powerlaw fluid with consistency k and index n. Therefore, the apparent viscosity in Eq. (30) is defined as follows: The shear rate, as an invariant of strain rate tensor, is expressed asγ = and can be calculated locally by the non-equilibrium particle distribution function using Eq. (25).
In the non-Newtonian fluids, viscosity is not constant, so relaxation parameter is calculated at each lattice site by recalling the relationship between fluid viscosity and relaxation parameter: Therefore, the relaxation parameter obtained by Eq. (32) is used in Eq. (14).The above equation is an implicit equation which can be solved by iterative method numerically. However, in this study,γ(τ) is updated by the relaxation parameter at the previous time step [35,38]. Such approach uses much less computational cost.
Apparent viscosity corresponds to shear-thinning fluid, shear-thickening fluid, and Newtonian fluid when n is less than 1, more than 1, and equal to 1, respectively. Therefore, when shear rate increases, the viscosity of shear-thinning fluid decreases. As such, the shear rate in this study is computed by particle distribution function rather than using derivatives of velocities by finite difference scheme.
This approach, deriving parameters locally, is more consistent with the philosophy of the LBM [35].
Direct implementation of a power-law fluid flow in the LBM confronts an important problem. In a shear-thinning fluid (n < 1) and at zero shear rates (γ = 0), the relaxation parameter diverges. On the other hand, in a shear-thickening fluid (n > 1) and at zero shear rates (γ = 0), the relaxation parameter becomes zero, leading to instability in solution as well. Hence, in this study the truncated powerlaw viscosity model is implemented to limit the range of viscosity as follows [28]: Therefore, in the truncated power-law model, range of relaxation parameter is limited and prevents solution from instability.

Effect of variable mobility
In this section, we examine the effect of mobility in our simulations. We set our variable and constant mobility as in Eq. (6). As can be seen, the C i (1 − C i ) term, which was applied in [39], is multiplied by the constant mobility introduced in [25]. In fact, instead of a constant mobility applied to the whole domain, a variable mobility is implemented on the interfaces only. Figure 1 shows initial condition of three components so that red, green, and blue colors correspond to component 1, 2, and 3, respectively. Green bubble rises because of buoyant force to stick to the red droplet. This situation, bubble-droplet rising, is a critical situation in our simulation which can cause mass loss. Table 1

Single phase power-law fluid flows between two parallel plates
In this subsection and the following one, two benchmarks are conducted to validate the implementation of our code. First benchmark is a pressure-driven flow between two fixed parallel plates at a distance of 2L in the y-direction while fluid flow moves in the x-direction as shown in Fig. 4. The analytical solution for steady state velocity profile with power-law index is as follows: [6] u where u in is the inlet velocity, y is the distance from the center of two parallel  The relative L 2 -norm of the velocity has been calculated at steady state condition between analytical and numerical results as presented in the following equation: where u x is the x component of velocity, u analytical x is the x component of analytical velocity, and E rel is the relative L 2 -norm of the velocity.   ||u(x,t)|| 2 < 10 −8 . As discussed in [31], the analytical velocity in the x-direction is obtained as follows: It should be noted that the analytical maximum velocity is u * x (y = 0). The maximum velocity should not be greater than 10 −2 because of compressibility error [31]. Several simulations with different λ P and power-law indices are performed.
Analytical and numerical velocity profiles are compared as shown in Fig. 7. The relative L 2 -norm which was presented in Eq. (35) is obtained in each simulation to investigate the accuracy of this model. Table 3 shows the relative errors of the simulated and analytical velocity. As it can be seen in the table 3, λ P and λ are used for non-Newtonian and Newtonian systems, respectively. The relative errors are in order of 10 −2 in each simulation which shows that the numerical results are in good agreement with the analytical solution. It should be noted that as long as the maximum velocity decreases, the compressibility error will be also smaller, therefore, the relative errors are lower.  Table 3: Relative errors of the velocity for different power-law indices and different λ P and λ in the two-phase powerlaw layered flow driven by pressure gradient.

Drag coefficient and dimensionless numbers of a bubble rising in a liquid
In this section, the terminal velocity of a bubble rising in a power-law fluid is investigated. When the bubble rising (subscript g) in a liquid (subscript l) reaches its terminal velocity, buoyancy and drag forces exerting on the bubble become equal. Therefore, the drag coefficient can be written as [10,11]: The Reynolds number, Morton number, and Eötvös number in power-law fluids are defined as [11,42] Re = ρ l U 2−n d n e k (38) Where d e is the equivalent diameter which is defined as d e = 3 6V b π using bubble total volume, V b . In the Newtonian fluids, the drag coefficient is only a function of Reynolds number, but in the power-law fluids the index n has also effect on the drag coefficient (C D = f (Re, n)) [42].
Li et al. [11] proposed a drag correlation for a bubble rising in a shear-thinning liquid, considering deformation of bubble and rheological properties of fluids besides Reynolds number. They showed that the proposed correlation is capable of handling spherical and deformed bubble.
where E is the aspect ratio of the bubble, expressing the ratio of vertical and horizontal diameters of the bubble. They showed that this correlation have great agreement with experimental results under the finite range of 0.05 < Re < 300.
It should be noted that in our study, a bubble rising by buoyant force, the bubble 27 experiences E ≥ 1 because the only force exerted on the bubble is from bubble rear. But if a bubble rises in other non-Newtonian liquids applying forces from both left and right like viscoelastic fluids, bubble aspect ratio can be less than unity.

Grid independency
Grid independency test in four domains is performed to find the grid independent computational domain. A domain of L × 2L is considered which L is the length of it. The no-slip boundary condition is imposed on the top and bottom and periodic boundary condition is applied on the left and right boundaries. Time evolution of gravity center of the single bubble, which is under buoyant force, is plotted for each grid in Fig. 8. Table 4 shows rheological properties of liquid and physical properties of gas-liquid system, where τ g , d g , and g represent relaxation time of the disperesed component, the rising bubble, diameter of the bubble, and gravitational acceleration, respectively.  The parameters used in the Table 4 are selected so that dimensionless numbers, Reynolds, Eötvös, and Morton numbers, remain constant in each grid. The gas-liquid system is assumed to be quiescent initially and initial condition for concentration is as follows: In this study, the velocity of the bubble and gravity center of the bubble are recorded in average-form as follows: Also, we define dimensionless form for velocity and location for better comparison by dividing them to gd g and 2L, respectively. To choose a suitable domain, dimensionless velocity and center of gravity of the bubble are examined. Magnitude of dimensionless terminal velocity of the bubbles are 0.198, 0.243, 0.268, and 0.273 for grid 1, grid 2, grid 3, and grid 4, respectively. Since the error between grid 3 and grid 4 is less than 1%, grid 3 is suitable for the rest of our simulations.
Also, curves correspond to the grid 3 and 4 are so close as shown in Fig. 8

Bubble rising in a shear-thinning liquid
In this section, the terminal velocity of a bubble in a shear-thinning liquid, whose viscosity decreases by applying shear stress, is investigated. To investigate pre-   Velocity of the bubble increases gradually and reaches a constant value (steady state). To show this process, for instance, in case 4, lattice unit and dimensionless velocity of the bubble are shown in Fig. 9 (a) and (b), respectively.  To compare our results with the mentioned correlation, a relative error, (|numerical− correlation|/correlation), for each simulation is provided in Table 7 with the maximum relative error of 19%, minimum relative error of 0, and average relative error of 5.6%. Table 7 shows that terminal velocity of the bubble at different conditions are in good agreement with empirical correlation proposed by Li et al. [11]. One of the reasons behind small deviation between the results is that in the present model, the shear-thinning liquid is assumed to be pure, unlike experimental results which were done using carboxymethyl cellulose (CMC) solutions.
Case n E U t (Li et. al [11]) U t (present) Relative Error  In the case 1, the surface tension is very low; therefore, the bubble tends to deform more as shown in Fig. 11   In the case 2, the physical properties are the same as the case 1 but the surface tension. Therefore, the final shape of the bubble in the case 2, which has more surface tension, is less deformed with respect to the final shape of the bubble in the case 1 as shown in Fig. 12 (a). Also, Fig. 12 (b) demonstrates that by increasing surface tension, the kinematic viscosity around the bubble increases because a low deformation indicates a low shear rate. In the case 3, Reynolds and Eötvös numbers decrease by increasing consistency index and surface tension coefficient, respectively. As a result, the bubble is capable of maintaining spherical shape over time as shown in Fig. 13 (a). Also, Fig. 13 (b) shows that kinematic viscosity distribution is higher with respect to the cases 1 and 2. Also, the two regions on the left and right of the bubble, having large viscosity, is bigger with respect to the cases 1 and 2. In the cases 4 and 5, higher power-law indices, n = 0.6 and n = 0.63, are implemented while all other parameters are the same. Since in our simulation the shear rate magnitude is less than unity, shear thinning behavior increases by increasing power-law index [43]. Therefore, as Fig. 14 (b) and Fig. 15 (b) shows, n = 0.63 exhibits more shear-thinning characteristics than n = 0.6. In addition, due to the kinematic viscosity distribution in the domain, the Reynolds number in the case 5 is approximately two times larger than that in the case 4.  In some cases, for instance case 6, as shown in Fig. 16 (b), two regions at the top and bottom have constant, high kinematic viscosity because these regions experience minimum and equal shear rate at specific times. When the bubble rises, it exerts large shear stress around itself, resulting in a decrease in the viscosity of surrounding shear-thinning fluid. Therefore, the bottom region retrieve its high, constant kinematic viscosity when the bubble take away from it, and this constantviscosity region get larger as long as the bubble goes upward. As a result, when the bubble gets close to the top boundary, shear rate becomes variable and larger at the top region, transforming the top region to a small, high-viscosity region. The rising bubble pushes the surrounding fluid upward, forcing it to circulate.
As a result, two wakes appear, experiencing low velocity as shown in Fig. 11 (c)-16 (c), which make the bubble rear flattened or skirted, especially in high Eötvös numbers.

Conclusions
In this paper, a power-law scheme was implemented in a LB model to capture both shear-thinning and Newtonian behavior of a fluid. A variable mobility was introduced in this study to reduce interface dissipation between phases. A comparison between constant and variable mobility in a three-component system was presented, indicating that the variable mobility can preserve mass of components better.
Two cases were investigated to test the accuracy of the model. The dynamics of a rising bubble in the shear-thinning liquid was examined numerically using a LB model. Results show that the rising bubble deforms at high Eötvös number and reaches a steady state shape and velocity. Also, the bubble moves upward faster and has higher location at high Reynolds number. The terminal velocity obtained by an experimental drag coefficient was compared to the simulated terminal velocity, demonstrating that the results of the present study are in good agreement with the experiment. The average, maximum, and minimum relative errors for six cases were 5.6%, 19%, and 0, respectively. When the bubble rises in a shear-thinning liquid, the viscosity around the bubble decreases because of the increase of shear rate. Therefore, the viscosity nearby the bubble, which experiences the highest shear rate, is the lowest and it increases gradually when gets farther away from the bubble. Also, the regions which are too far from the bubble have constant and maximum viscosity. The velocity contour and streamlines show that the rising bubble generates two wakes around itself, containing low or negative velocity.
In future study, we will use multiple relaxation time (MRT) to capture higher viscosity ratios. Also, higher Reynolds numbers should be captured to investigate wider range of dimensionless numbers.

Conflict of interest
The authors declare no conflict of interest.

Data Availability
The data used to generate the results in the manuscript are available from the corresponding author upon reasonable request.