One-dimensional model of turbulent flow of non-Newtonian drilling mud in non-prismatic channels

One-dimensional model of non-Newtonian turbulent flow in a non-prismatic channel is challenging due to the difficulty of accurately accounting for flow properties in the 1-D model. In this study, we model the 1-D Saint–Venant system of shallow water equations for water-based drilling mud (non-Newtonian) in open Venturi channels for steady and transient conditions. Numerically, the friction force acting on a fluid in a control volume can be subdivided, in the 1-D drilling mud modelling and shallow water equations, into two terms: external friction and internal friction. External friction is due to the wall boundary effect. Internal friction is due to the non-Newtonian viscous effect. The internal friction term can be modelled using pure non-Newtonian viscosity models, and the external friction term using Newtonian wall friction models. Experiments were carried out using a water-based drilling fluid in an open Venturi channel. Density, viscosity, flow depth, and flow rate were experimentally measured. The developed approach used to solve the 1-D non-Newtonian turbulence model in this study can be used for flow estimation in oil well return flow.


Introduction
One-dimensional prediction of non-Newtonian turbulent effect is more challenging than 2-D and 3-D shallow water flow prediction. One-dimensional models are, however, considerably more economical. There can be two types of friction assumptions in non-Newtonian fluids: internal friction due to viscous effect and external friction due to channel boundaries (Jin and Fread 1997, Fread 1988, 1993. External friction from channel walls can, in 1-D modelling, be calculated from the Darcy-Weisbach equation, the Chezy equation, and the Manning formula (Manning 1891;Chow 1959;Akan 2006;Abdo et al. 2018). This is similar to the Newtonian flow friction force from the walls. The Manning formula is the most widely used (Rahman and Chaudhry 1997;Sanders and Iahr 2001;Agu et al. 2017;Welahettige et al. 2018). The open-channel flow friction factor can be expressed as being equivalent to the pipe flow friction factor, pipe diameter being replaced by four times the openchannel hydraulic radius for Newtonian flow (Chow 1959;Akan 2006;Alderman and Haldenwang 2007). f = 16∕Re * is widely used for the rectangular-channel friction factor for a fully developed non-Newtonian laminar flow. R e * is here a generalization of the Reynolds number (Kozicki and Tiu 1967;Burger et al. 2010). There are in general two types of laminar flow regimes in open-channel flow when Re < 500 , and there are small flow depth and small flow velocities: subcritical laminar and supercritical laminar (Chow 1959). Laminar flows are, however, not significant in large-scale flow applications such as oil well return open-channel flow. Turbulent flow is, however, easily propagated due to high flow rates, wall friction, shape of the channel, and viscous forces.
Internal friction from the non-Newtonian fluid flow in open channels can be modelled using pure non-Newtonian flow models such as the power law (Kozicki and Tiu 1967) and the Herschel-Bulkley model (Jin and Fread 1999;Haldenwang 2003). A number of non-Newtonian turbulent open-channel friction factors have been reviewed by Alderman and Haldenwang (2007). According to the dip phenomenon (Stearns 1883), maximum velocity in open channels takes place below the free surface in narrow channels with an aspect ratio of l∕h < 5 (Sarma et al. 1983;Yang et al. 2004;Bonakdari et al. 2008;Absi 2011). Where the bed is rough, the curvature of the velocity distribution increases due to the weak secondary motion from the lateral solid walls, transporting low momentum fluid to the central section (Nezu et al. 1994). The dip phenomenon is not widely used in 1-D modelling due to the difficulty of the formulation. Rectangular channels are very common in 1-D shallow water equation modelling. Trapezoidal open channels are less widely modelled because the trapezoidal shape of the cross section increases the complexity of the equations. Mozaffari et al. (2015) and Liu et al. (2017) have studied time-dependent properties of non-Newtonian fluid.
Supercritical and subcritical flow regimes can, where the channel is horizontal, be observed occurring simultaneously before and after the Venturi region (Welahettige et al. 2017a, b). Higher-order Godunov-type numerical schemes are recommended for solving the open-channel conservation equations due to the following issues: unsteady hydraulic jump propagation, to avoid negative flow depth (reduce numerical viscosity) and maintain stability at dry or near-dry conditions (Sanders and Iahr 2001;Kurganov and Petrova 2007). The flux-limitercentred (FLIC) scheme using the source term splitting method is well suited to solving 1-D shallow water equations (Welahettige et al. 2018). The FLIC scheme is used to calculate the interface fluxes, the lower-order flux and higher-order flux being combined using a flux limiter function. The higher-order flux comes from the Richtmyer scheme, and the lower-order flux from the first-order-centred (FORCE) scheme, which is a combination of the Lax-Friedrichs and the Richtmyer schemes (Toro 2009).
A large number of studies have been conducted into drilling mud pipe flow (Alderman et al. 1988;Bailey and Peden 2000;Maglione et al. 2000;Piroozian et al. 2012;Livescu 2012;Aslannezhad et al. 2016). There are, in contrast, fewer published studies on drilling mud flow in open channels. The primary objective of this research paper is therefore to validate the 1-D numerical model for drilling mud in open non-prismatic channels flow using experimental results. The developed models will be used in the future for well return flow estimation. Model accuracy depends on the validity of the assumptions. Pure non-Newtonian models are combined with the turbulence models to provide a source term for the centred total variation diminishing (TVD) scheme. Viscosity, density, flow depth, and flow rates are measured experimentally at the laboratory scale for a water-based drilling mud.

Numerical schemes
The shear rate variation from the bottom wall to the free surface can be formulated in 3-D and 2-D models by the velocity gradient ̇z x = u∕ z where z ≤ h . Here, ̇z x is the shear rate in the x-direction perpendicular to the z-direction. It is, however, challenging to include the shear rate into a 1-D model. Three main shear stresses apply in 3-D fluid flow in the x-direction: a linear elongation deformation ( xx ) and two shear linear deformations ( zx and yx ). Linear elongation deformation can, assuming incompressible liquid properties, be neglected (Versteeg and Malalasekera 2007). Bottom surface velocity becomes zero under the no-slip condition, and shear stress from the air is negligible at the free surface. There are two strong velocity gradients for Newtonian fully developed turbulent open-channel flow: the inner region, which is approximately 20% of the flow depth, and the outer region (Bonakdari et al. 2008). The dip correction factor can be neglected if the aspect ratio is > 5 and the velocity profile is similar to the log law for smooth walls (Yang et al. 2004;Absi 2011). The shear stress effect from the side wall from wide, open channels is smaller than the shear stress from the bottom walls, yx ≪ zx . Therefore, yx can be neglected for 1-D (Longo et al. 2016). The average shear stress and shear rate for 1-D models, based on the assumptions made, can be considered to be ≈ zx , and ̇≈̇z x , respectively. The shear stress correlated to the power law (PL), the Herschel-Bulkley (HB), and the Pierre Carreau (PC) model can be given as The 1-D shallow water equations need to include a modification for the contraction and expansion region of open Venturi channels, to avoid artificial accelerations (Fread 1993;Sanders and Iahr 2001;Welahettige et al. 2018). An additional term has been suggested for the general shallow water equations in our previous study to accurately take into account the non-prismatic effect of the channel walls (Welahettige et al. 2018). The new term is a function of the flow depth and the variation of channel bottom width, k g 2 h 2 g b∕ x . For prismatic channels ( b∕ x = 0 ), the shallow water equations are converted to the general shallow water equations. In this study, we, however, attempt to formulate the non-Newtonian friction slope as a separate term for 1-D shallow water equations, S i . One-dimensional shallow water equations in Saint-Venant's form for locally trapezoidal channels and for non-Newtonian fluid can therefore be presented as: Unlike rectangular channels, the cross-sectional area ( A ) is not a linear function of flow depth ( h ) (Fig. 1). The relation between the flow depth and the cross-sectional area can be expressed as (Welahettige et al. 2018) k g 1 is the ratio between the gravity height of the crosssectional area and the flow depth in the trapezoidal shape. This helps calculate an accurate hydrostatic pressure (Welahettige et al. 2018). k g 2 is the ratio between the gravity height of the sidewall cross-sectional area and the flow depth. The affected sidewall area is approximately of rectangular shape. Therefore, k g 2 ≈ 0.5 . However, k g 1 is a function of b and h . Therefore, using a constant value for k g 1 is not valid. It is always higher than 0.5 (Welahettige et al. 2018): According to the turbulent pipe flow of a non-Newtonian fluid, shear stress can be formulated as a combination of the effect of dynamic viscosity ( ) and eddy momentum kinematic Douglas et al. 2001;Chhabra and Richardson 2011). This approach can be used for open-channel non-Newtonian turbulent flow. In turbulent flow, wall shear stress is predominant in the laminar region and turbulent eddies are predominant in the turbulent core. Laminar sublayer thickness is, in general, very small in turbulent flow (Versteeg and Malalasekera 2007). The laminar sublayer effect can therefore essentially be described by a pure laminar rheological model, and the turbulent core effect can essentially be described by the Newtonian turbulence model. The Manning formula adds wall friction by considering the hydraulic radius of the channel for the numerical model of the 1-D shallow water equations (Chow 1959). According to our previous study, Manning's turbulence friction model produces good results for open-channel turbulent water flow (Welahettige et al. 2018). The friction slope for the turbulent The hydraulic radius for trapezoidal channels is The Reynolds number for pipe flow Re = VD∕ can be converted into an open-channel Reynolds number by replacing D with 4R h . D is here the pipe diameter, and is the effective viscosity of a non-Newtonian fluid. According to the Hagen-Poiseuille equation, 8V∕D is the shear rate at the wall for Newtonian or non-Newtonian pipe flow (Chhabra and Richardson 2011). For the open channel, the shear rate for non-Newtonian flow is ̇≈ 2V∕R h (Haldenwang 2003). We use this shear rate to describe the shear stress. The shear rate ̇≈ 2V∕R h is valid only for laminar flow. The turbulent effect is, however, included in Manning's formula. Here, we assume that the laminar region is dominated by internal friction and that the turbulent region is dominated by external friction. This assumption is valid in the open channel due to the higher shear rates at the wall boundary and lower shear rates at the free surface. Kozicki and Tiu's (1967) power-law-based Reynolds number, Zhan and Ren's, Albulanga's and Naik's (Alderman and Haldenwang 2007) Bingham plastic-based Reynolds numbers, and Slatter's (1995) and Haldenwang's (2003) Herschel-Bulkley-based Reynolds numbers are widely used for open-channel non-Newtonian flow. Using the same approach, Reynolds numbers for open-channel flow can be derived from the power law, the Herschel-Bulkley, and the Carreau viscosity models. Effective viscosity = ∕̇ is taken from Eqs. (1)-(3), for laminar region channel flow, then the friction force due to the non-Newtonian viscous effect is F i = b + 2k 2 k g 1 h Δx . This is based on the assumption that average shear stress applies to the gravity height of the flow depth in a control volume. The internal friction slope can be introduced as a function of the internal friction factor,S i = A gf i . The dimensionless non-Newtonian friction factor can be introduced as (11) For a rectangular channel, the non-Newtonian friction factor then becomes f i = ∕( gh) where k 1 = 0 and k 2 = 0 . Jin and Fread (1997) derived a similar non-Newtonian friction factor for mud fluid in a rectangular open-channel flow. Non-Newtonian friction factors for the power law, Herschel-Bulkley, and Carreau fluids can be derived as follows: Equations (4)-(5) are solved using the FLIC scheme and Runge-Kutta fourth-order explicit scheme, for rectangular channels of water (without internal friction slope) by Welahettige et al. (2018). The FLIC scheme and Runge-Kutta fourthorder explicit scheme are used for solving the advection term and the source terms, respectively. This method also extends to the 1-D turbulent non-Newtonian fluid. In this study, we implement the FLIC scheme and Runge-Kutta fourth-order scheme for the turbulent non-Newtonian fluid in a trapezoidal shaped channel. Flow rate Q can be given as Q = AV , where the average velocity across the cross section is considered to be u ≈ V . The area perpendicular to the flow direction is a function of time, flow depth, and spatial domain A = A(t, h, x) , and the average velocity is a function of time and spatial domain V = V(t, x) . The pure advection Eq. (15) is solved with conserved variables u 1 = A and u 2 = AV . For continuous bottom topography channels, the bottom-width variation effect is highlighted in the TVD solving method used here. This can be compared with the conventional centred TVD solving method, 1∕Δx R ( ) (Welahettige et al. 2018). The pure advection term (advection flux and wall reflection effect) is solved first using the centred TVD method. Here, (13) (14) For PL, HB, and PC models, the source terms are , respectively, with their friction slopes, the internal friction slopes being S iPL = A gf iPL , S iHB = A gf iHB , and S iPC = A gf iPC . The wall reflection effect R ( ) m j is solved here as an advection term, due to the simple numerical calculation of the b∕ x term, and to minimize numerical diffusion (Welahettige et al. 2018). One advantage of using the centred TVD scheme is to avoid the strictly hyperbolic requirement of the partial differential equations. m is here the time index, m ∈ {1, 2, … , N} . j is the node index in the spatial grid, j ∈ {1, 2, … , l} . The source terms (gravity effect, external friction, and internal friction) are then solved using an ordinary differential equation (ODE) solver, the explicit Runge-Kutta fourth-order method (Toro 2009). The initial condition for the ODE solver is the solution from the centred TVD scheme.k g 1 u 1 , h u 1 , R h u 1 , S e u 1 , u 2 , S iPL u 1 , u 2 , S iHB u 1 , u 2 , and S iPC u 1 , u 2 can be derived in terms of conserved variables by substituting u 1 and u 2 with A and AV.
The apparent viscosity of the shear rate is 500 1 s −1 , 500 , which is a constant. Here, k s is the roughness height and is 15 µm in this study, which is similar to the value for steel walls. Fread (1988;Jin and Fread 1997) has derived a friction slope due to internal viscous dissipation with the rheological properties of the power law equation and a yield stress (similar to the Herschel-Bulkley model). This includes a semi-empirical velocity profile. According to Fread's model, the internal friction factor for trapezoidal channels can be presented as The source for the Fread model is Here, l is the free surface width. For trapezoidal channels, l becomes b + 2k 1 h.
In the sequel, PL, HB, PC, Haldenwang, and Fread models are solved by using the FLIC scheme and Runge-Kutta fourth-order method. The only differences are the source terms. We call the models as PL, HB, PC, and Fread where they are combined with Manning's friction.

Venturi rig
The complete flow loop of the rig contains a mud-mixing tank, a mud-circulating pump, an open Venturi channel, and a mud return tank. See Figs. 2 and 3. The sensing instruments in the setup are a Coriolis mass flow meter, pressure transmitters, temperature transmitters, and ultrasonic level transmitters. Chhantyal et al. (2017) and Agu et al. (2017) also conducted experiments using the same experimental setup. Level transmitters are located along the central axis of the channel and can be moved along the central axis. The accuracy of the Rosemount ultrasonic 3107 level transmitters is ± 2.5 mm for a measured distance of less than 1 m (Welahettige et al. 2017b). The accuracy of temperature transmitter is ± 0.19 °C at 20 °C. The accuracy of Coriolis mass flow meter is ± 0.1%. All the experimental values presented in this paper are averaged values of level sensor readings taken throughout a period of 5 min at each location. The channel inclination can be changed. A negative channel inclination ( angle) indicates a downward direction. The dimensions of the trapezoidal channel are shown in Fig. 3; the main flow direction is in the x-direction.

Viscosity and density measurements
The water-based drilling mud used for the experiments contained potassium carbonate as a densifying agent and xanthan gum as viscosifier (Chhantyal 2018). The drilling mud viscosity and density were measured using an Anton Paar MCR 101 rheometer and an Anton Paar DMA 4500 density meter. At the beginning of the experiments, the viscosity meter reached a constant room temperature, 25 °C, within 20 min. The constant temperature was maintained throughout the experimental period. For a fixed shear rate value, 40 measuring points were taken within 800 s. The averages of 40 measuring points were considered in this study. The standard deviation of viscosity is small, 1 × 10 −5 , within the 40 measuring points. This indicates that the rheometer reaches a steady state. The combined uncertainty of viscosity is 0.015 mPa s, which is calculated based on the quantifying uncertainty in analytical measurement (QUAM) method (Ellison et al. 2000). An external Anton Paar Viscotherm VT 2 cooling system has a standard temperature uncertainty of 0.02 K (Idris et al. 2017). For DMA 4500, the temperature accuracy is ± 0.03 K. The uncertainty of density is determined as 0.34 kg/m 3 (Han et al. 2012). The experimentally measured viscosity data were fitted to Eqs. (1)-(3) based on nonlinear regression techniques: the power law model, the Herschel-Bulkley model, and the Carreau model as shown in Fig. 4. The curve-fitted parameters from Table 1 are used: Herschel-Bulkley model for Haldenwang, Fread, PL, HB, and PC models. The density of the drilling mud is 1336 kg m −3 . Experimental viscosity measurements are in the range of the rheometer accuracy limit where the shear rate is 100-1500 s −1 . A shear rate of less than 100 s −1 shows a significant variation in all the models. The accuracy may not be significant in the open Venturi channel flow where the shear rate is less than 100 s −1 for the fluid we used in the experiment. This is because flow regimes are turbulent and the average shear rate is higher than 100 1 s −1 . This will be further analysed by comparing PL, HB, and PC model results. The average errors between  Maglione et al.'s (2000), and Chhantyal's (2018) drilling fluid data are also taken from the literature and used for the comparison the calculated and the experimental values are 6.14%, 7.7%, and 4.4%, respectively, for the power law model, the Herschel-Bulkley model, and the Carreau model with 100 s −1 <̇< 1500 s −1 . According to the nonlinear least squares approach in MATLAB R2018a, coefficients are calculated with 95% confidence bounds. The R-squared values are above 0.95 for all the fitted models. The calculated yield stress is 0.1451 Pa based on the Herschel-Bulkley model. The yield stress is a small value for the fluid used in this study. The open-channel flow is highly turbulent and has a high Reynolds number; the small yield stress will not significantly influence the flow regimes in the open-channel flow at high turbulence level.
The drilling f luid used in this study shows shear thinning properties in the range of shear rate 100 s −1 <̇< 1500 s −1 (Fig. 4).
Real and experimental drilling fluids from the literature are used for further comparison in the open-channel flow modelling. Khodja et al.'s (2010) and Maglione et al.'s (2000) real drilling fluids rheology based on the Herschel-Bulkley model is shown in Fig. 4 and Table 1. Chhantyal (2018) used a model drilling mud and the experimental setup used in this study. The model drilling mud rheology was given in terms of the power law model.

Results and discussion
The rheological parameter and flow parameters used in the simulations are shown in Tables 1, 2, and 3.

Steady results
Haldenwang, Fread, PL, HB, and PC model results are compared with experimental flow depth along the channel at a steady state (Fig. 5). The results are steady state, reached from an unsteady condition. At the beginning, the drilling fluid enters the empty channel at a constant inlet flow rate. Steady results are achieved after 310 s. The critical flow depth ( h c ) is 40 mm at 0 ≤ x ≤ 2.95 and 3.45 ≤ x ≤ 3.7 , and bottom width and flow depth are the same in the two ranges (Welahettige et al. 2018). According to the Froude number, the flow regime is subcritical ( h > h c ) before the Venturi region and supercritical ( h < h c ) after the Venturi region (Fig. 5). All the models derived in this study (PL, HB, and PC) give similar results, which indirectly imply that the curve-fitted rheological parameters for each model act in the same way as for the drilling fluid used in this study. The fluid accumulates before the Venturi contraction when the channel is horizontal. Then, due to the channel contraction effect, a hydraulic jump moves upstream before the steady state is reached. The flow depth increases and the velocity decreases due to high energy loss, the friction models giving Table 1 Rheological parameters of the drilling muds used in this study. Figure   higher friction. The Haldenwang model shows higher friction in this case than the PL, HB, and PC models. According to Eq. (17), wall roughness height for the channel is assumed to be 15 µm, the value for steel walls. This wall roughness value was tested in our previous study (Welahettige et al. 2017b). The PL, HB, and PC models, however, used the Manning roughness factor for steel-smooth walls, which is 0.012. The average deviation from experimental results is 5% in PL, HB, and PC models at a steady state. The PL model results are further compared with the Chhantyal's (2018) experimental results (Fig. 6). The rheology of the drilling fluid is given in terms of the power law model. They used a mechanical filter to remove foam, flow depth, therefore being less influenced by foam than in our experimental results. The flow depth difference between experiment and simulation is reduced in the subcritical region when the foam is removed. The PL model gives a good prediction of the Chhantyal's (2018) experimental results. The average deviation from the experimental result is 8%. The yield stress of the fluid is comparatively small. The viscosity of the fluid used in Chhantyal's work is 0.04 Pa s at 1 s −1 shear rate, as shown in Fig. 4. The average viscosity value calculated in the simulation is 0.005 Pa s at a steady state, at an average shear rate of 450 s −1 . This implies that, in this case, the open-channel flow regimes do not reach the small shear rate ranges. The effect of low shear rates might therefore be insignificant. Figure 7 shows flow depth variation with time for step changes in the channel inlet flow rate. The pump outlet flow mass rate varies between 10 and 40 kg min −1 from the set point. At the beginning of the experiment, the channel flow rate is 470 kg min −1 at a steady state.

Unsteady results
Step changes are carried out for the set point of the pump flow rate at time t = 64 s and t = 188 s. We show here two-level sensors readings,  They are fixed at x = 2.12 m and x = 3.2 m from the inlet of the channel, above the free surface along the channel central axis. The PC model and the Haldenwang model results are compared with the experimental readings at the dynamic condition. LT-18 is located after the Venturi region. LT-15 is located before the Venturi region. Even though sudden step changes in the flow rate occur at time t = 64 s and t = 188 s, the experimental level reading gradually changes flow depth, a complete step change taking more than 40 s. This is despite the pump having a small time constant of 1.6 s. This is due to the time required for fluid to travel from the inlet of the open channel to the level sensor locations, and to unstable wave propagation. Higher flow depth is shown in LT-15 than in LT-18 due to the hydraulic jump formation travelling upstream before the Venturi contraction. The PC model results give a higher accuracy  The accuracy of the models depends on the assumptions used in model development, the curve-fitted rheological parameter, and boundary conditions. An assumption for the 1-D models is that velocity is only considered in the x-direction, V x . However, the velocity component V z is comparatively small due to free surface movement in the z-direction, which is restricted by the surface tension of the fluid and gravitational force, V x ≫ V z . Channel sidewalls are balanced with the velocity component in the y-direction, V x ≫ V y . However, contraction and expansion of side walls influence and change the direction of the flow path (Welahettige et al. 2017b). The average aspect ratio was 4 in this study. The effect of the dip phenomena was therefore assumed to be small. In this study, the average error of flow depth is 0-6%.
One of the main advantages of using a 1-D model compared to a 3-D model is less execution time. According to our 3-D CFD simulation of the same case (related to Fig. 5 for water), the 1-D model took 1 min to execute and the 3-D CFD model took more than 5 h (Welahettige et al. 2017b).

Effect of source terms
There are no direct experimental results for the derived friction slopes in this study. The internal and external friction slopes are therefore calculated for different drilling fluids. The three drilling fluids all have different densities and viscosities, the rheology of the fluids being given in the Herschel-Bulkley model. Figure 8 shows a drilling fluid flow depth comparison. The Herschel-Bulkley parameters that specify the drilling fluid rheology are shown in Table 3. According to the viscosity and shear rate curves, for a given shear rate ( > 40 s −1 ), fluid viscosity ranked from the highest to lowest is the Khodja et al.'s (2010) fluid, the Maglione et al.'s (2000) fluid, and the fluid of this study. The inlet flow rate was kept constant at 0.0056 m 3 s −1 for all the fluids. The different densities, however, imply that the inlet mass flow rate is different for each fluid. As explained above, a flow depth that is larger than the critical flow depth is subcritical and contrast supercritical. The Reynolds number is higher than 5300 throughout the channel. Flow regimes therefore become subcritical turbulent and supercritical turbulent. The Reynolds number is calculated from Eq. (10) by substituting effective viscosity for Eq. (3). The channel inclination angle is − 1.7°. Experimental flow depth and the PC model results are well matched for the entire region of the channel. The model drilling fluid used in this study creates an oblique jump at the Venturi throat (Welahettige et al. 2017b). The two other drilling fluids show a hydraulic jump formation at the quasi-steady state. The fluid used in this study, however, has a lower viscosity than the two other drilling fluids and the highest density. This causes lower energy loss and high mass flow rates. The fluid used in this study gives supercritical flow regimes throughout the channel. However, the inlet supercritical flow regimes cannot be maintained for the high viscous fluids (Khodja's and Maglione's drilling fluids), because they generate a hydraulic jump and the level rises to keep the same flow rate at the steady state. Khodja's drilling fluid has a higher density than Maglione's drilling fluid. Khodja's drilling fluid does, however, show greater movement of the hydraulic jump front in the upstream direction due to the higher viscosity of the fluid. The friction slopes used to calculate the flow depth in the drilling fluid are studied further in Fig. 9.
The calculated friction slope terms, internal friction slope ( S i ) and external friction slope ( S e ), are shown in Fig. 9. The external friction term is the highest for the lowest viscous fluid, and the internal friction is the lowest for the lowest viscous fluid (the fluid used in this study, which relates to the HB model). Flow depths are supercritical and subcritical before and after the hydraulic jumps. In the subcritical region, the internal friction slope is predominant. In the supercritical region, internal and external friction terms actively contribute to numerical calculations.

Conclusion
In this study, a 1-D non-Newtonian turbulent model for non-prismatic open-channel flow was developed based on non-Newtonian rheological models and Newtonian turbulence models. The fluid friction term can be divided into two terms: internal friction and external friction. Internal friction is due to non-Newtonian viscosity and external friction is due to wall friction. The higher-order FLIC scheme and Runge-Kutta fourth-order method were used to solve the new 1-D non-Newtonian turbulence models. The approach used to solve the 1-D non-Newtonian turbulence model in this study can be used for flow estimation in oil well return flow. The flow depth prediction error varies from 2 to 8% in this study, depending on the model's assumptions and experimental results. The internal friction term is predominant in subcritical flow because laminar flow regimes participate in improving the viscous forces, at a steady state. The external friction and internal friction terms contribute to supercritical flow regimes.