Stepwise mathematical derivation of the Herschel–Bulkley laminar fluid flow equations—in pipes

Stepwise derivation of flow equations of the Herschel–Bulkley (HB) model is not available in the literature. These equations are crucial for mechanical, chemical and petroleum engineering academia and industries where fundamental works on non-Newtonian fluids may be done to reach future models and estimation methods. Therefore, this work focuses on derivation of laminar flow equations and estimation methods of HB fluids through pipes. In this work, first, stepwise derivation of the HB fluid flow parameters consisting of fluid velocity, flow rate, average velocity and relative velocity equations is presented, followed by a straightforward mathematical model for use in numerical solution. Next, stepwise mathematical derivation of the laminar pressure drop equations by Merlo et al. (An innovative model for drilling fluid hydraulics. Paper presented at the SPE Asia Pacific oil and gas conference, Kuala Lumpur, Malaysia, 1995) and Gjerstad and Time (SPE J 20:1–18, 2014) is presented, and finally practical and user-friendly calculation procedures for different estimation methods are presented. The step-by-step derivation procedures presented in this work contribute to effective learning for engineering students and practitioners in addition to providing a clear example derivation guideline for future researchers to reach other more accurate non-Newtonian hydraulics models and estimation methods.

Effective diameter defined by Reed and Pilehvari (1993) for non-Newtonian flow, i.e., the diameter of a circular pipe that would have the identical pressure drop for flow of a Newtonian fluid with viscosity equal to the apparent viscosity and the same average velocity as the non-Newtonian flow [m] D hyd Hydraulic diameter used for transformation of pipe flow to annular flow equations [m] f YP A parameter defined by  and  almost equal to DP HÀP G Coefficient first used by Reed and Pilehvari (1993)    and  U HB Dimensionless velocity of the HB model fluid defined by  and  U PL Dimensionless velocity of the PL model fluid defined by  and Gjerstad  Greek letters DP HÀP The difference between the wall shear stress ratios in the HB and PL models (P HB and P PL , respectively) as defined by  and     and  for estimation of pressure drop in the HB fluid n A constant used by  and  for estimation of pressure drop in the HB fluid w A constant used by  and  for estimation of pressure drop in the HB fluid (for pipe flow only) c Shear rate [1/s]

Introduction
There are several non-Newtonian rheological models available in the literature. The most well-known ones are the Bingham plastic (BP) model; power-law (PL) model; Herschel-Bulkley (HB) model, also called the modified or yield power-law model; the Robertson-Stiff (RS) model; and the Sisko model. The BP model and the PL model are two-parameter rheological models and are still commonly used in mechanical, chemical and petroleum industries. However, some fluids especially drilling fluids, slurries or suspensions are not well approximated by these two models. Some modifications to these models may be applied to increase their prediction accuracy, and an example is a work done by Ashena et al. (2021) by applying a correction coefficient; in another example, roles of different affecting parameters on accuracy were quantified (Ashena et al. 2022). Alternatively, non-Newtonian fluids can be better represented by the HB model (Hanks 1979) and the RS model (Robertson and Stiff 1976). Both HB and RS models are 3-parameter rheological models which are both viscoplastic (with a yield stress before flow can be initiated) and shear thinning, i.e., their apparent viscosity reduces with increasing shear rate (Gjerstad et al. 2012). The RS model can be applied to describe behavior of cement slurries and drilling fluids, but it is slightly less accurate than the HB model (Mitchell and Miska 2011). The Sisko model was mentioned in a study by Weir and Bailey as a recommended model for drilling fluids (Sisko 1958, cited by Weir andBailey 1996). Although the HB model was introduced long time ago in 1926 (Herschel and Bulkley 1926), it yields mathematical expressions which are not readily solvable. Therefore, several works such as Chin (2012) focused on numerical estimations. Besides, its semi-analytical pressure drop equations were not found for long time. Instead, numerical methods were used to iteratively calculate pressure drops; Chilton and Stainsby (1998) presented a simpler integration version of that first presented by Govier and Aziz (1972). In addition, there was not a consensus on greater accuracy of HB advantage over the other models. Being considered a complex model, its use was not widespread particularly in the drilling industry. A turning point for the HB model was year 2006 when the American Petroleum Institute (API) recognized and announced the HB model as the most accurate rheological model and recommended its usage for the petroleum industry. This trend change by API is recognizable by comparing API (2003) with API (2009). Since then, applications of the HB model increased significantly in drilling hydraulics. Some serious works have been done on the HB model during decades. Metzner andReed (1955) andMetzner (1957) carried out some pioneering works in 1950s on non-Newtonian fluid flow through pipes. They developed a method called generalized shear rate which provided a way to generalize their results of PL fluids to all time-independent non-Newtonian fluids including HB model (cited by Whittaker and Staff 1985;Reed and Pilehvari 1993). Govier andAziz (1972, 2008) and Cheng (1975) presented an implicit equation for the HB fluid flow in pipes. Mitchell and Miska (2011) and  and  provided the velocity profile of HB model fluids in laminar flow, but clear stepwise proofs were not provided. Indeed, explicit pressure drop equations have not been found in their works yet. Instead, just estimations of the pressure drop in the laminar flow regime were found. Merlo et al. (1995) presented an approximation equation of the HB model's pressure drop equation, but without providing any proof. In their paper, they presented the calculation methodology in a practical example. Their work was later cited by Guo and Liu (2011). However, again no derivation of their equations was given to prove where they come from. In 2014, using two limiting cases and an analogy between the HB and PL models,  presented another approximation equation for the HB's fluid pressure drop equation in the laminar flow regime. Their methodology is based on taking limitations from special cases without providing proofs to the limits taken. In addition, their methodology is complex and userunfriendly to apply in a real example.
In summary, previous researchers working on laminar equations of HB fluids have not provided stepwise derivation of pipe flow equations and user-friendly estimation methods for pressure drop calculations. These matters are indicative of the gaps in the literature filling which is not only important to improve effective learning for engineering students and practitioners, but also to provide a clear example procedure guideline useful for future researchers to reach other more accurate non-Newtonian hydraulics models and estimation methods. Therefore, in this work, the laminar pipe flow equations of the HB model are derived. The equations consist of fluid flow velocity, flow rate, relative velocity and pressure drop equations. The simplest mathematical model for use in numerical solution is provided. Derivation of equations for estimated pressure drops by Merlo et al. (1995) and  is provided in a user-friendly manner for practical field or academia's applications. Finally, using an example, results of pressure drop estimations by the different approximation methods are compared with the numerical solution and the best matching one is identified. It is noted that the forms of equations presented in this work are in SI units.

Fluid velocity
First, a relation between shear stress and pressure drop is given for fluid flow through pipes. Mitchell and Miska (2011) have presented a derivation of the equation without showing all steps of the derivation. Here, we present the derivation in a more stepwise and simple manner. Assuming a force balance for steady flow (i.e., no acceleration), the following balance relation holds between the shear force stress (due to shear stress s) and the normal force (due to pressure drop dP): Therefore: Applying Eq.
(2) to the radius at which there is an unsheared portion of the fluid (r p ) and also to the radius at the wall (R): The equation of shear stress versus shear rate for the HB model is: Combining Eqs.
(1) and (4), we have: Rearrangement and using a dummy variable ''m'' instead of ''1/N'': Integrating from both sides: Journal of Petroleum Exploration and Production Technology (2023) 13:625-643 627 Thus: Using the boundary condition of V ¼ 0 at r ¼ R: Thus: Combining Eqs. (9) and (10), the annular velocity is: It is noted that Eq. (11) holds for ðr p r RÞ; however, there is a constant velocity plug flow or V p in the central part of the pipe ð0 r r p Þ as shown in Fig. 1. It is noted that there is a plug zone in HB fluids as is the case with any non-Newtonian fluids with yield stress. In this zone, the fluid moves as solid bodies within sheared flows (Chin 2012). The radius of this zone is: Replacing ''r'' in Eq. (11) with ''r p '' from Eq. (12) gives: Comparing Eqs. (10) and (13), it is inferred that Combining Eqs. (11) and (13), the pipe flow velocity can be generally expressed as follows: The above equations have been presented by Mitchell and Miska (2011) with a typo.

Flow rate and average velocity
In this section, we find the flow rate of flow in pipes. We know the following relation holds between incremental flow rate dq and velocity V:: To find the flow rate, we need to integrate the velocity with respect to the flow area as follows: As the velocity in the plug region (V 1 ) is constant in the plug flow region, q 1 is: Replacing r p from Eqs. (12), (13) becomes: Next, it is essential to evaluate the second component of the flow rate, q 2 : Taking the constant coefficients out: There are two integrals within q 2 to solve as follows: The solution of q 2À1 is simply: Replacing the original term of s y ðdP=dLÞ=2 for r p gives: Journal of Petroleum Exploration and Production Technology (2023) 13:625-643 629 To solve q 2À2 , a minor rearrangement leads to: Solving q 2À2 is not straightforward and requires using change of variables as follows: The term s y ðdP=dLÞ=2 is the same as the constant r p . We temporarily set it so in the integral to simplify the form of the integration. Thus: Performing the polynomial integration while using the equivalent integral limits from r ¼ r p to U ¼ 0, and r ¼ R to U ¼ R À r p , we have: Reverting U parameter and constant r p , respectively, back to their originals ðr À s y ðdP=dLÞ=2 Þ and s y ðdP=dLÞ=2 : Using Eqs. (19) and (22), q is found: Rearrangements give: Multiplying ðdP=dL=2Þ 3 by the numerator and denominator gives: gives: Inserting ð 1 1þm Þ into the bracket and multiplying by the right-hand side terms yields: Adding and subtracting some terms to the right-hand side: Factoring ''ð½dP=dLR=2Þ 2 À 2½dP=dLR=2 Â s y þ s y 2 '' into ð½dP=dLR=2 À s y Þ 2 gives: Factoring the coefficients of ½dP=dLR=2 À s y Þ 2 gives: Þ '' and some rearrangement: Rearrangement of the coefficient for 2s y ð½dP=dLR=2 À s y Þ gives: Replacing ''m'' with its equivalent ''1/N'' gives: Equation (23) is the same equation presented by Mitchell and Miska (2011). Therefore, this work presented the methodology of reaching this equation in a stepwise manner.
For ease of notation, we prefer to have a positive sign for q in the equations and thus we assumed that s y and ''dP=dL'' are all positive. However, to prevent any chance of confusion, we prefer to replace ''dP=dL'' notation with ''DP=DL'': As Eq. (24) shows, the flow rate q ¼ f ð DP DL ; R; s y ; NÞ depends on the following four parameters: pipe radius (R), flow behavior index (n) which is found using laboratory measurements, yield stress (s y ) which is measured at the laboratory, and pressure drop per length ( DP DL ). Normally, the flow rate will be available together with the geometry characteristic (radius) and fluid properties (s y and N), with pressure drop being the output interesting for engineering applications. This is possible either numerically or by approximation methods.
The average velocity is found by dividing q by the crosssectional area of the pipe (R 2 p): Flow rate and average velocity expressed in terms of shear stress Knowing ''½DP=DLR=2'' is equal to the shear stress at the wall s w , we can write Eq. (24) as follows: Next, we convert ''ðdP=dL=2Þ 3 '' to s w by multiplying the numerator and the denominator by ''R 3 '': Therefore, the average velocity V is: Note: For purpose of solving the HB's model numerically, it is beneficial to assume s w s y ¼ X where X is considered as the unknown parameter to be found numerically by a computer. Therefore, we have: Mitchell and Miska (2011) made further modifications in order to write the average velocity in terms of the Bingham number '' s y s w .'' We also use their term in a more stepwise manner. To do that, we first factor out a coefficient of s w 1þ1=n . Next, multiplying by the cross-sectional area, the flow rate q is found as: Allowing application of s w 2 to the parentheses gives: Next, we make a simplification of the form of Eq. (29) by using a dummy variable called generalized Bingham number / (Peixinho et al. 2005). Knowing / ¼ s y s w , we have: Using the equivalent of ''''½DP=DLR=2'' for s w gives: Factoring ''ðN þ 1) (2N þ 1) (3N þ 1)'' gives: Combining the polynomial terms in the bracket {} gives: To find the relative velocity profile, first the average velocity of the fluid in the pipe is found by dividing q in Eq. (31) by the cross-sectional area of the pipe ''R 2 p'': Relative velocity profile First, the velocity equations are essentially expressed in terms of shear stress. As the first step, we multiply and divide Eq. (13) by R 1þm : Knowing s y ¼ ðdP=dLÞR=2 and replacing }m} by its equivalent ''1=N'': Again, using the dummy variable / for s y s w (i.e., Bingham number): Next, the same process above is repeated for V 2 . Thus, we first multiply and divide Eq. (11) by R 1þm : Replacing }dP=dLR=2} with }s w } and }m} with }1=N} and s y s w with /: Having simplified and found V 1 and V 2 ðrÞ in terms of /, we can summarize the relative velocities as follows: Replacing dP=dL R 2 with s w : The relative V 2 V is found by: Summarizing Eqs. (37) and (38) for relative velocities gives: The above equations are presented by Mitchell and Miska's equation 5.123-a and b (2011) with a typo.

Pressure drop equations
There are four important approximation/estimation methods in the literature for HB fluids' pressure drops which are discussed as follows (with stepwise derivation of the first two ones): Merlo et al. (1995) and Guo and Liu (2011) presented laminar pressure drop equation of HB fluids for flow in the pipe, but they did not present any derivation. This task is handled in this work: Replacing s w with s HB and K with K HB Eq. (27) for the average velocity in the HB model:

Merlo et al. (1995)
The pressure drop equation of a PL fluid DP DL for pipe flow is: The shear stress of the PL fluid s PL is: Therefore, Case (a) s HB ! ¥ (infinite shear stress and shear rate) Therefore, when s HB ! 1, the limit is: When s HB ! 1, the limit is: In this imaginary case, N tends to zero. In this case, the fluid behaves like a solid with constant shear stress. In this case, the limit is: With N ¼ 0, we have: s HB ¼ s y þ K and K ¼ s PL . In other words: s HB ¼ s y þ s PL . Next, Merlo et al. (1995) have apparently assumed that as the HB's flow behavior index tends to zero (N ! 0), the PL model's flow behavior index would also tend to zero (n ! 0). This assumption may cause some error in low shear rates, but that has been ignored. Therefore, we have: Cancelling out ''ð 1 N Þ} by ''n'' and replacement of n and N equal to zero: Considering the above limits for the investigated cases (a, b and c), an equation structure can be made for s HB Àsy K HB 1=N s PL K PL 1=n as follows: where We already showed in Eq. (41) that for flow through pipes, s PL K PL 1=n is: Therefore: Cancelling out terms: Taking both sides to the power of N: Therefore, s HB is: we can find ðDP=DLÞ HB as follows: Some further rearrangements give: Merlo et al. (1995) and Guo and Liu (2011) Þs PL þNs y Nþ1 ð Þðs PL þs y Þ .'' C c is the circular or pipe correction coefficient. Reed and Pilehvari (1993) called a similar term Þs PL þNs y Nþ1 ð Þðs PL þs y Þ R'' as the effective diameter.

DP DL
This is the same equation given by Merlo et al. (1995) and Guo and Liu (2011), except that they replaced average velocity ''V'' with its equivalent '' q pR 2 '' and ''R'' with ''D=2.'' The correction coefficient C c is simplified as follows: Replacing the equivalent of s PL using Eq. (41) gives the same equation as in Merlo et al. (1995) or Guo and Liu (2011): It is noted that the above equations are in SI units.

Gjerstad et al. (2014)
Similar to Merlo et al. (1995),  used comparison of HB model prediction with that of the PL model. However, they first proceeded to non-dimensionalizing the velocity and shear stress parameters in the HB model. Next, they related the HB model parameters to the PL one. Next, having found some limiting cases, a general structure was reached as the pressure drop equation of the HB model. The provided equations are not step-by-step and straightforward. Therefore, in this section, a stepwise derivation is provided: First, the dimensionless forms of two important parameters of the pressure gradient (DP=DL) and average velocity (V) are defined as follows: According to Eq. (27), the average velocity of a HB fluid in terms of shear stress is rewritten as: Taking V to the power of N: Dividing both sides by R N K : Dividing both sides by s y and defining it as U HB : Dividing the numerator and denominator of the righthand side by s y 3Nþ1 : Simplifying gives: Multiplying the numerator and the denominator by s w s y : Inserting the coefficient ''N N ð sw sy À1Þ Nþ1 ð sw sy Þ 3Nþ1 '' into the parenthesis and replacing ''1/N'' by ''m'':  used dimensionless forms of HB model parameters of P HB and U HB : Therefore: Next, for the PL model, the following relationship holds for the pressure gradient in laminar flow: Therefore, the shear stress in the PL model s w;PL is: Dividing both sides by s y : Configuring the dimensionless ''P PL ¼ s w;pl s y ,'' we have: Or: Knowing the dimensionless U PL : Using equations Eqs. (54) and (55): Or, U PL is:  made an assumption here (without mentioning) that K HB ¼ K PL and N ¼ n, and therefore U HB ¼ U PL : Thus: # N Therefore, we can write: where e ¼ ðP HB À 1Þ 3 P HB They defined DP HÀP as the difference between P HB and P PL as follows: Combining Eq. (58) and Eq. (60):  took limitations from Eq. (61) for the following cases: Case (a) Flow Initiation (s w ! s y orP HB ! 1) First, we know that the limit of e is: In short: Case (b) Very High Flow Rates (s w ! ¥ orP HB ! ¥ ) First, we know that the limit of e is: lim The first part of the above limitation has the indeterminate form. First, we convert the indeterminate form from ''1 Â 0'' to '' 0 0 '' so that we can later use the L'Hopital's rule of limitation (Thomas and Weir 2004): Now, the L'Hopital's rule of limitation is applied by taking derivation from both the numerator and the denominator: The derivative de dP HB is found as: Therefore, using the above derivative, the limit is found as: lim P HB !1 P HB 1 À e N À Á ¼ lim Cancelling out P HB 2 , knowing already that lim P HB !1 e ¼ 1, and expanding the terms in the parentheses: Therefore, lim P HB !1 DP HÀP is: In brief: Finally, using the limits of DP HB found in cases (a) and (b),  found an approximation equation structure which meets the limits. The equation structure had some constants (n, r and w). Next, they performed tuning and found the constants of the equation. A summary of their equations is presented in a stepwise manner as follows: Step-1: Calculate the dimensionless pressure in PL model, P PL (using the average velocity V): Step-2: Calculate the transition parameter P T , and f 0 : Step-3: Calculate f YP : where n ¼ 0:97 À 0:1N À 0:11N 2 Step-4: Calculate P HB : Step-5: Calculate DP DL : Generalized shear rate method Metzner and Reed (1955) and Metzner (1957) introduced the ''Generalized Shear Rate'' method based on the effective diameter concept in order to expand PL models to non-Newtonian models. The mathematical background behind this model is trivial, and it is simply using analogy to the PL model. The effective diameter for non-Newtonian pipe flow is the diameter of a circular pipe that would have the identical pressure drop for flow of a Newtonian fluid with viscosity equal to the apparent viscosity and the same average velocity as the non-Newtonian flow (Reed and Pilehvari 1993). The procedure of this method is as follows: Step-1: Calculate constant G: Step-2: Calculate the effective diameter D eff : Step-3: Calculate the wall shear rate c: Step 4: Calculate the pressure gradient DP DL : where N is the HB model's flow behavior index.
Results and discussion Table 1 shows the summary of HB's model equations for pipe and annular flow, as derived in this work.

Conclusions and future work
Stepwise derivation of laminar flow equations and estimation methods of Herschel-Bulkley (HB) fluids through pipes are not available in the literature. Therefore, this work was focused on finding and presenting step-by-step derivation procedures in order to contribute to effective learning for engineering students and practitioners in addition to providing a clear example derivation guideline for future researchers to reach other more accurate non-Newtonian hydraulics models and estimation methods.
The following conclusions were taken: 1. Stepwise and detailed mathematical derivation of the HB fluid flow velocity, average velocity, relative velocity and flow rate was presented. 2. Different approximation methods of HB laminar flow pressure drop were presented in forms of user-friendly procedures applicable to field and academic use. 3. In a future work, these estimation methods will be compared in some case studies to find the most accurate one.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Funding This research work had no fund.