A New Solution of Drag for Newtonian Fluid Droplets in a Power-Law Fluid

: Understanding flow behaviors of multiple droplets in complex non-Newtonian fluids is crucial in many science and engineering applications. In this study, a new and improved analytical solution is developed based on the free surface cell model for the flow drag of swamp of Newtonian fluid drops through a power-law fluid. The developed solution is accurate and compares well to the numerical solutions. The improvement involves a new quantification of shear stress boundary condition at the interface and a more consistent approximation in linearizing the power-law fluid flow governing equation. The Newtonian fluid solutions can be reasonably used to linearize the flow governing equation. The approximation of the boundary conditions at the interface, however, has a major impact on the model prediction. The main improvement in the new solution is observed under the condition of comparable viscosities of the Newtonian drops and the outside power-law fluid when the results are sensitive to the interface boundary condition. Under the two extreme conditions of high viscosity ratio (approaching particles) and low ratio (approaching bubbles), the present and existing solutions converge.


Introduction
The interaction and motion of dispersed liquid droplets and gas bubbles in another continuous phase liquid are basic processes encountered in many applications related to the manufacture of foodstuffs, foams, polymers, cosmetics, toiletries, and many other products.In other applications, the falling or rising of liquid droplets and gas bubbles in another fluid is of critical importance.An adequate understanding of the hydrodynamic aspects of the motion and interaction of multiple droplets and bubbles is needed for the modeling of transport processes in such multiphase systems [1].Since the motion of droplets and bubbles is characterized by the drag force and drag coefficient, efforts have been made to determine the effects of physical and geometrical parameters on the drag.
In recent decades, it has been realized that most fluids do not follow the simple Newtonian fluid behavior in terms of the relation between stress and deformation rate.These fluids are known as non-Newtonian fluids.One common type of non-Newtonian fluids is the pseudoplastic fluid whose apparent viscosity decreases with the increasing shearing rate, which is also called shear-thinning behavior [2].Some non-Newtonian fluids, known as viscoelastic fluids, exhibit both viscous and elastic characteristics during deformation since they have the ability to store energy and thus show partial recovery upon the removal of stress, while other non-Newtonian fluids may have time-dependent characteristics, in which stress varies with the duration of flow [2].Therefore, it is highly important to obtain knowledge of droplet behaviors in not only a Newtonian fluid but also in a non-Newtonian one, because the investigation of the droplets in non-Newtonian fluids provides useful and essential information for optimal process designs and operations.In recent decades, there have been many studies to examine the behaviors of multiple flow problems involving non-Newtonian fluids.The approaches to the analytical modeling of the motion of droplets, particles, and bubbles in non-Newtonian liquids were summarized with emphasis on the non-Newtonian models and their effects [3].More recently, the studies on the behaviors of drop motions through non-Newtonian fluids have been reviewed [4].
Slow flows of pseudoplastic Carreau fluid over Newtonian spherical drops and bubbles were solved numerically to assess the effects of pseudoplastic behavior, holdup, and viscosity ratio on the drag and to evaluate the effects of the Newtonian plateau seen in the viscosity functions of some pseudoplastic fluids [5,6].Kishore et al. [7] elucidated the role of power law rheology on the velocity of an ensemble of Newtonian droplets translating in a power-law fluid by numerically solving the governing equations at moderate Reynolds numbers.Flow and drag characteristics of a single bubble and a swarm of spherical bubbles in contaminated power-law fluids were numerically investigated using the spherical stagnant cap model [8,9].The gas holdup of bubble swarms in shear-thinning fluids was experimentally studied at various superficial gas velocities [10].
There have been many theoretical studies to solve flow problems of a pseudoplastic fluid over bubbles, drops, and particles, in recent decades.Among them, many studies have adopted the free surface model originally proposed by Happel [11] as a conceptualized representation of flow problems over multiple particles, droplets, or bubbles to simplify the interacting effects among them.The free surface cell model was used in combination with variational principles to obtain the upper and lower bounds on the drag coefficient of a swarm of Newtonian fluid drops in a power-law fluid [12] and in a Carreau fluid [13].Similarly, the combination of free surface cell model and variational principles was adopted to predict the rising velocity of spherical bubbles in a Carreau fluid [14] and in a powerlaw fluid [15].The free surface cell model and variational principles were also used to determine the drag bounds of flows of a power-law fluid [16] and a Carreau fluid [17] through an assemblage of solid particles.However, the variational principle results have the main drawback that the predictions are sensitive to the choice of trial functions used in the analysis and it is impossible to estimate a priori the potential errors inherent in the predicted bounds.
The free surface cell model has been utilized to derive approximate analytical solutions of non-Newtonian pseudoplastic fluid through the assemblage of bubbles, drops, and particles in recent decades.For complex non-Newtonian fluids, exact solutions are generally not possible.An approximate approach was proposed to linearize the governing equation of a non-Newtonian fluid; this approach uses Newtonian fluid solutions in part of the non-Newtonian fluid equation so that analytical solutions become possible [18].Based on the linearization, approximate solutions of power-law fluid flow through an assemblage of particles and bubbles were derived [19].An analytical solution of power-law fluid flows over a swarm of bubbles was developed by applying a similar but simpler and consistent approximation in linearizing the flow governing equation [20].The linearization approach was also applied to obtain an analytical solution for a power-law fluid flow over a swarm of drops [21].The flow problem of a power-law fluid through a swarm of Newtonian fluid drops was solved by using the finite-difference method [22].
Previous analytical solutions deviated significantly from numerical solutions for fluids with strong pseudoplastic effects.In this study, a new and more accurate analytical solution is developed for creeping power-law fluid flow over a swarm of Newtonian fluid drops.In particular, the solution involves a new treatment of the shear stress boundary condition at the interface and a more consistent approximation in linearizing the power-law fluid governing equation.The developed solution is compared with the finite-difference numerical solutions from the literature and the existing solutions, and the differences among the models are then discussed.

Methods
The drag of swamp of Newtonian fluid drops in a power-law fluid is solved using the free surface model [11] of two concentric spherical cells.The inside cell is a spherical Newtonian fluid droplet and the outside cell is a power-law fluid (Figure 1).The radius of the outside cell is then determined by the requirement that the volume fraction of drop in the cell of two concentric spheres should be equal to the holdup of the system being considered.For the free surface model, the shear stress on the outside fluid cell surface is zero (i.e., the fluid surface is frictionless).The free surface cell model has been widely used for estimating the rise velocity of swarms of spherical bubbles and drops [5-7, [23][24][25].For the rise velocity of a swarm of spherical bubbles in power law liquids, the upper and lower bound results using the cell model generally agreed with the scant experimental results in the literature [2].
Newtonian fluid droplet and the outside cell is a power-law fluid (Figure 1).The radius of the outside cell is then determined by the requirement that the volume fraction of drop in the cell of two concentric spheres should be equal to the holdup of the system being considered.For the free surface model, the shear stress on the outside fluid cell surface is zero (i.e., the fluid surface is frictionless).The free surface cell model has been widely used for estimating the rise velocity of swarms of spherical bubbles and drops [5-7, [23][24][25].For the rise velocity of a swarm of spherical bubbles in power law liquids, the upper and lower bound results using the cell model generally agreed with the scant experimental results in the literature [2].
In the free surface cell model adopted in this study, the size of the cell is determined by the overall droplet holdup requirement where R2 is the radius of the inside droplet, R1 is the radius of the outside fluid cell, and e is the holdup, which is the volume fraction of droplets, in the system.For a power-law fluid in the outer cell, the dimensionless governing equation is [22] For the inside Newtonian fluid droplet, the corresponding governing equation is In the free surface cell model adopted in this study, the size of the cell is determined by the overall droplet holdup requirement where R 2 is the radius of the inside droplet, R 1 is the radius of the outside fluid cell, and e is the holdup, which is the volume fraction of droplets, in the system.For a power-law fluid in the outer cell, the dimensionless governing equation is [22] For the inside Newtonian fluid droplet, the corresponding governing equation is where (and throughout the text) the subscripts 1 and 2 correspond to the outer cell of a power-law fluid and the inside cell of Newtonian fluid, ψ is the stream function, ξ is the dimensionless radial distance, θ is the polar angle, and n is the power-law fluid index.In the above equations, the operator is For the free surface cell model, the boundary conditions are specified as follows where v ξ , v θ are the velocity components, D mij represents (i, j) the component of the rate of deformation tensor for the outside cell (m = 1) and inside cell (m = 2), respectively, and a = µ 2 / K(V/R 2 ) n with V being the superficial velocity, and K the consistency index.
With the governing equation given in Equation (2) and the boundary conditions given in Equations ( 5)-( 9), a new approximate solution for the flow and drag is developed in this study.
Under the weakly non-Newtonian behavior condition, all the terms in Equation (2) that are multiplied by (n − 1) or (1 − n) are approximated using the solutions when the outside cell is also a Newtonian fluid, expressed in Equations ( 10)-( 14) below, to linearize and simplify the governing equation given in Equation (2).In the study of Jarzebski and Malinowski [21], the two terms involving the third-order derivative of ψ 1 in the second term on the left-hand side were retained in the governing equation.All the terms on the right-hand side were, however, evaluated using the Newtonian fluid solutions as an approximation.Since these two terms involving the third-order derivative of ψ 1 on the lefthand side are also multiplied by |n − 1|, they are treated in the same way as the Π 1 terms by using the Newtonian fluid solutions in this study.Therefore, the governing equation in the present study is simpler than that in the study of Jarzebski and Malinowski [21].
For Newtonian fluid drops in another Newtonian fluid with a viscosity ratio b = µ 2 /µ 1 , the analytical solutions could be derived and expressed as [21] where the coefficients N 1 , N 2 , N 3 , and N 4 are Applying the Newtonian fluid results for the terms involving (n − 1) or (1 − n) in Equation (2) and only retaining the highest order term of ξ or q as they are always greater than 1, one can obtain the following dimensionless flow governing equation of a power-law fluid in the outside cell Fluids 2024, 9, 99 5 of 11 The general solution to Equation ( 15) can then be obtained as The Π 1 term in Equation ( 7) is evaluated using the Newtonian fluid results since it is raised to the power (n − 1)/2, which can be expressed as One of the main differences in this study from that of Jarzebski and Malinowski [21] lies in the evaluation of Π 1 in the boundary condition in Equation ( 7).Note that all the terms lower than N 3 ξ −2 were neglected and it was also approximated that cos 2 θ = 1 in the study of Jarzebski and Malinowski [21] while applying the boundary in Equation ( 7).While it is true that N 3 ξ −2 has highest order in most regions in the outside cell, it can be shown that 2N 1 − N 3 − 3N 4 = 0 when the inside cell is a solid particle and N 1 + N 4 = 0 when the inside cell is a bubble at the interface where ξ = 1.Therefore, at ξ = 1, simply retaining only one term involving N 3 in the evaluation of Π 1 for the boundary condition at the interface could be a source for significant errors.Therefore, in this study, all the terms in Equation ( 17) are included for the evaluation of Π 1 at ξ = 1.On the other hand, approximating cos 2 θ as 1 could also introduce large errors when the average of cos 2 θ is only 0.5, with θ varying from 0 to π.In this study, the average value of 0.5 for both cos 2 θ and sin 2 θ is adopted, and the following equation is used for the stress boundary condition Equation (7): In Equation ( 18), the factor c has the following expression: Jarzebski and Malinowski [21] only retained the highest order term of q and assumed cos 2 θ and sin 2 θ were approximated to be equal to 1 and 0, respectively.As a result, the factor c in Equation (18) had a simpler form, as follows: In this study, both approaches of approximating the stress boundary condition at the interface discussed above are examined and quantitatively compared.The solutions based on Equations ( 19) and (20) are denoted as "model 1" and "model 2", respectively, in the subsequent analysis and comparison.
From the required boundary conditions, the coefficients in Equation ( 16) can then be determined as where N = n (n − 1)N 3 .
The drag, F D , for the power-law fluid flow through a swarm of Newtonian fluid droplets can then be determined: From the solutions, the pressure and stress distributions in the flow field can be derived.The flow drag, F D , can then be determined analytically as The correction factor for the drag coefficient, Y D , is defined as where C D = F D / 1 2 ρV 2 πR 2 2 and R e = ρV 2−n (2R 2 ) n /K.The Y D can then be derived as In summary, one new aspect in the developed solutions lies in the way of dealing with the stress boundary condition and the quantification of the factor c in Equation (18).To analyze the differences, both expressions of the factor c in Equations ( 19) and (20) (i.e., model 1 and model 2, respectively) are evaluated and discussed.These two models and the model of Jarzebski and Malinowski [21] are compared with the numerical solutions by Zhu and Deng [22].The finite-difference technique was adopted to obtain numerical solutions of Equations ( 2) and (3) and the associated boundary conditions Equations ( 5)-( 9) by Zhu and Deng [22].The central space differences were used and the resulting finite-difference equations were solved using the successive over-relaxation method.The governing equations of both the inside cell Newtonian fluid flow and the outside cell power-law fluid flow were iteratively solved, since two flow systems are coupled through the boundary conditions in Equations ( 5)-(7).Since the numerical study solved the same governing equations and boundary conditions without any terms being neglected, the direct comparison with the numerical solution can quantify the extent of errors in the analytical solutions due to the omission of high order terms and the linearization of the governing equation using the Newtonian fluid results, and offer insights on the dominant factors that affect the flow drag on the droplets.The performance and difference of analytical model predictions are then quantified and discussed in the next section.

Results and Discussion
Figure 2 shows the comparison of the two models (model 1 and model 2) developed in the present study and the model of Jarzebski and Malinowski [21] with the numerical solutions of Zhu and Deng [22] when the drop holdup is e = 0.6.Three viscosity ratio (i.e., parameter a) values of 0.1, 1, and 10 are used in the comparison, which span a spectrum from close to bubbles (a = 0.1) to close to solid particles (a = 10).The difference between model 1 and model 2 lies in the treatment of the shear stress boundary condition in Equation (7).In model 1, all the terms are retained in the evaluation of the second invariant of the rate of deformation tensor Π 1 in the boundary condition at the interface of the droplet and the outside non-Newtonian fluid (i.e., at ξ = 1).In model 2, however, only the highest ξ order term is used in the evaluation of Π 1 .In the outside region, where ξ is greater than 1, neglecting the lower order term of ξ might be reasonable.At the interface where ξ is equal to 1, however, the lower order terms of ξ could be as important as the highest order term of ξ, and the simple treatment of the shear stress boundary condition in Equation ( 7) by retaining only the highest order could lead to significant errors.It can be seen from the comparison that the solutions of model 1 compare more favorably with the finite-difference numerical solutions.If the power-law fluid viscosity and the Newtonian fluid droplet viscosity have the same order of magnitude when the value of a is in the order of 1, the improvement using model 1 over the other two models is especially significant, which indicates that the treatment of the stress boundary condition at the interface is important.For a = 1, the solution of Jarzebski and Malinowski [21] could significantly over-estimate the drag at n = 0.6, while model 1 is closest at over-predicting by less than 4%.The prediction from model 2 from the present study is only slightly better that that from Jarzebski and Malinowski [21], which illustrates that the treatment of the boundary condition at the interface is the main contributor of drag prediction errors, while the contribution of approximating the governing equation is secondary.When a = 1 and the viscosities of the Newtonian fluid droplets and outside powerlaw fluid are comparable, the difference of model 1 in the present study and Jarzebski and Malinowski [21] is largest.In this case, the contribution from both normal stress and shear stress to the drag is comparable.In the study of Jarzebski and Malinowski [21], only the highest order term on the normal stress portion was retained.When a is large (approaching particles) and a is small (approaching bubbles), the effect of approximating the c factor diminishes and, therefore, the results from this study and Jarzebski and Malinowski [21] converge.The comparison of the correction factor of the flow drag of a power-law fluid through a swarm of Newtonian fluid droplets among the models, with the finite-difference numerical solutions at three different levels of drop holdup when a = 1, is shown in Figure 3.The closest agreement with the finite-difference numerical solutions is observed in model 1.
The main difference between model 2 and the model of Jarzebski and Malinowski [21] is the different linearization approximation in the main governing equation expressed in Equation (2).The model results of the correction factor for the drag coefficient from Jarzebski and Malinowski [21] and those from model 2 are close to each other, with model 2 being slightly better.Compared to the study of Jarzebski and Malinowski [21], the flow governing equation expressed in Equation (15) in the present study is the result of further simplification by using the Newtonian fluid solutions in evaluating all the terms that are multiplied by the factor (1 − n) or (n − 1).In the study of Jarzebski and Malinowski [21], however, the terms involving the second-order derivative of the unknown stream function were not approximated by using the Newtonian fluid results, although these two terms are also similarly multiplied by the factor (n − 1).Therefore, the resulting final approximate governing equation in Jarzebski and Malinowski [21] was more complex than that in the present study.With the same stress boundary condition approximation in dealing with Equation ( 7), the model of Jarzebski and Malinowski [21] produces similar results as model 2, as both significantly over-estimate the drag, which illustrates that the treatment of the stress boundary condition at the interface has the dominant effect on the model prediction.The results of the correction factor for the drag coefficient of a power-law fluid flow over a single droplet are obtained by using a small value of the drop holdup of e = 10 −10 , as shown in Figure 4.In the case of a single droplet, the radius ratio of the outside cell over the droplet radius (i.e., the q parameter defined in Equation ( 1)) approaches infinity, and retaining only the highest order produces similar results as including all terms in dealing with the stress boundary condition (Equation ( 7)) at the interface between the droplet and the outside power-law fluid.The difference among the models mainly comes Malinowski [21] is largest.In this case, the contribution from both normal stress and shear stress to the drag is comparable.In the study of Jarzebski and Malinowski [21], only the highest order term on the normal stress portion was retained.When a is large (approaching particles) and a is small (approaching bubbles), the effect of approximating the c factor diminishes and, therefore, the results from this study and Jarzebski and Malinowski [21] converge.
The results of the correction factor for the drag coefficient of a power-law fluid flow over a single droplet are obtained by using a small value of the drop holdup of e = 10 −10 , as shown in Figure 4.In the case of a single droplet, the radius ratio of the outside cell over the droplet radius (i.e., the q parameter defined in Equation ( 1)) approaches infinity, and retaining only the highest order produces similar results as including all terms in dealing with the stress boundary condition (Equation ( 7)) at the interface between the droplet and the outside power-law fluid.The difference among the models mainly comes from the different ways of linearizing the governing equation in Equation (2).Unlike the high droplet holdup cases, the results from all three models are closer to those from the finite-difference numerical results for the single drop case for the scenarios considered.Therefore, the approach of linearizing the governing equation works better under the low droplet holdup conditions.
Fluids 2024, 9, x FOR PEER REVIEW 10 of 12 droplet, retaining only the highest order term in the shear stress boundary condition treatment at the interface is more reasonable compared to the higher droplet holdup cases.Therefore, the difference among the models with varying approaches of linearizing the governing equation and approximating the boundary conditions is relatively small.Overall, the new approaches of approximating the boundary condition result in improvements in predicting the drag experienced by a swarm of droplets in a non-Newtonian powerlaw fluid.As expected, as the power-law fluid index deviates further from 1 (i.e., further from Newtonian fluid behavior), the results from all three models also show increasing difference from the finite-difference numerical solutions because the approximation errors are all proportional to |1 − n|.From the results and discussion above, the solutions developed in the present study perform well in comparison with the finite-difference numerical solutions.Since the drag experienced by the droplets is intimately related to the stress and pressure at the interface between the droplets and the surrounding non-Newtonian fluid, the accurate evaluation of the second invariant of the rate of deformation tensor at the interface is essential.Therefore, the new treatment of the stress boundary condition at the interface by including all For a single droplet in a non-Newtonian power-law fluid, the results of the correction factor of the drag coefficient from model 2 are closer to the finite-difference numerical solutions than those from the study of Jarzebski and Malinowski [21] when a = 10.On the other hand, however, the results of Jarzebski and Malinowski [21] are closer to the finite-difference numerical solution than model 2 when a = 0.1.For a single droplet, there is an improvement by using model 1, but the improvement is not significant, especially when the viscosity ratio is small, which means that the droplet is close to a bubble.For a single droplet, retaining only the highest order term in the shear stress boundary condition treatment at the interface is more reasonable compared to the higher droplet holdup cases.Therefore, the difference the models with varying approaches of linearizing the governing equation and approximating the boundary conditions is relatively small.Overall, the new approaches of approximating the boundary condition result in improvements in predicting the drag experienced by a swarm of droplets in a non-Newtonian power-law fluid.As expected, as the power-law fluid index deviates further from 1 (i.e., further from Newtonian fluid behavior), the results from all three models also show increasing difference from the finite-difference numerical solutions because the approximation errors are all proportional to |1 − n|.
From the results and discussion above, the solutions developed in the present study perform well in comparison with the finite-difference numerical solutions.Since the drag experienced by the droplets is intimately related to the stress and pressure at the interface between the droplets and the surrounding non-Newtonian fluid, the accurate evaluation of the second invariant of the rate of deformation tensor at the interface is essential.Therefore, the new treatment of the stress boundary condition at the interface by including all terms in the second invariant of the rate of deformation tensor and a simpler linearization approach for the governing equation used in the present study can produce an accurate prediction of drag compared to the finite-difference numerical solutions.It should be noted, however, that there are limitations inherent in the underlying assumptions.Since the Newtonian fluid results are used for the terms that are multiplied by the factor (n − 1), the solutions are only applicable when |n − 1| is small, which means that the pseudoplastic effect of the power-law fluid is not strong.

Concluding Remarks
In this study, a new and improved solution is developed for the drag of a non-Newtonian power-law fluid flow through a swarm of Newtonian droplets.Both the new solution and existing model in the literature are compared with the finite-difference numerical solutions.
The most significant improvement in the present study is observed under the condition of comparable viscosities of Newtonian droplets and the outside power-law fluid.
Using Newtonian fluid results in part of the terms to linearize and simplify the governing equation is reasonable.However, the approximation of the shear stress boundary condition has a major impact on the model prediction of drag when the droplets and non-Newtonian fluids have a similar magnitude of viscosity because the drag results are sensitive to the treatment of the interface shear stress boundary condition.
Under the two extreme conditions of high viscosity ratio (approaching particles) and low viscosity ratio (approaching bubbles), the drag results are less sensitive to the treatment of the interface shear stress boundary condition, and the present and existing solutions produce similar results.

Figure 1 .
Figure 1.Schematic diagram of a free surface cell model.

Figure 1 .
Figure 1.Schematic diagram of a free surface cell model.

Figure 2 .
Figure 2. Comparison of theoretical models with numerical solutions when the holdup is e = 0.6 at three different levels of viscosity ratio.Factor c is determined by Equation (19) for model 1 and is given in Equation (20) for model 2. "J&M" is Jarzebski and Malinowski [21] and "numerical" means the finite-difference solution of Zhu and Deng [22].

Figure 2 .
Figure 2. Comparison of theoretical models with numerical solutions when the holdup is e = 0.6 at three different levels of viscosity ratio.Factor c is determined by Equation (19) for model 1 and is given in Equation (20) for model 2. "J&M" is Jarzebski and Malinowski [21] and "numerical" means the finite-difference solution of Zhu and Deng [22].

Fluids 2024, 9 , 12 Figure 3 .
Figure 3.Comparison of theoretical models with numerical solutions when the viscosity ratio a = 1 at three different levels of holdup.Factor c is determined by Equation (19) for model 1 and is given in Equation (20) for model 2. "J&M" is Jarzebski and Malinowski [21] and "numerical" means the finite-difference solution of Zhu and Deng [22].

Figure 3 .
Figure 3.Comparison of theoretical models with numerical solutions when the viscosity ratio a = 1 at three different levels of holdup.Factor c is determined by Equation (19) for model 1 and is given in Equation (20) for model 2. "J&M" is Jarzebski and Malinowski [21] and "numerical" means the finite-difference solution of Zhu and Deng [22].When a = 1 and the viscosities of the Newtonian fluid droplets and outside power-law fluid are comparable, the difference of model 1 in the present study and Jarzebski and

Figure 4 .
Figure 4. Comparison of theoretical models with numerical solutions for a single drop in powerlaw fluid.Factor c is determined by Equation (19) for model 1 and is given in Equation (20) for model 2. "J&M" is Jarzebski and Malinowski [21] and "numerical" means the finite-difference solution of Zhu and Deng [22].

Figure 4 .
Figure 4. Comparison of theoretical models with numerical solutions for a single drop in power-law fluid.Factor c is determined by Equation (19) for model 1 and is given in Equation (20) for model 2. "J&M" is Jarzebski and Malinowski [21] and "numerical" means the finite-difference solution of Zhu and Deng [22].