Modeling and Analysis of Fluid Flow through A Non – Prismatic Open Channel with Application to Drilling

This paper presents the development and validation of a simplified dynamic model of a Venturi channel. The existing dynamic models on open channels are based on the open channel flow principles, which are the shallow water equations. Although there are analytical solutions available for steady state analysis, the numerical solution of these partial differential equations is challenging for dynamic flow conditions. There are many complete and detailed models and numerical methods available for open channel flows, however, these are usually computationally heavy. Hence they are not suitable for online monitoring and control applications, where fast estimations are needed. The orthogonal collocation method could be used to reduce the order of the model and could lead to simple solutions. The orthogonal collocation method has been used in many chemical engineering applications. Further, this has been used in prismatic open channel flow problems for control purposes, but no literature is published about its use for non–prismatic channels as per the author’s knowledge. The models for non–prismatic channels have more non–linearity which is interesting to study. Therefore, the possibility of using the collocation method for determining the dynamic flow rate of a non–prismatic open channel using the fluid level measurements is investigated in this paper. The reduced order model is validated by comparing the simulated test case results with a well–developed numerical scheme. Further, a Bayesian sensitivity analysis is discussed to see the effect of parameters on the output flow rate.


Background
Kick and loss detection in an oil drilling process is a key requirement to reduce the risk and cost in the drilling industry. In the drilling process, a drilling fluid is circulated from the tanks into the well through the drill bit, up the annulus and then back to the tanks, as shown in Figure 1. This drilling fluid is used to balance the pressure inside the well, to bring the drill cuttings out of the well and also to lubricate the drill bit. If the pressure inside is not well balanced, a kick or a loss can occur. When the pressure inside the well exceeds the formation pressure, the drilling fluids can seep into the formation and commonly known as a loss. When the well pressure is too low than the formation pressure, the formation fluids (air, water, oil or a combination) can flow into the well, producing a kick. A kick or a loss can be catastrophic, if not controlled. The difference of the fluid flow rate pumped in to the well and the return flow rate coming out of the well is a commonly used indicator of a kick/loss. However, the return flow rate is quite difficult to measure because it contains the drill cuttings as well as dissolved gas. A novel approach of using a Venturi rig for this purpose is proposed in previous studies done by Pirir et al. (2017); Jinasena et al. (2017); Chhantyal et al. (2016). This is a low cost and an easy to maintain solution. The flow rate through the Venturi can be estimated online with the use of mathematical models. However, the modeling and simulation of free surface flows are complex and challenging. Generally, the open channel hydraulics are often modeled by the wellknown and efficient shallow water equations which are also known as SaintVenant equations (Chow, 1959;Chaudhry, 2008;Litrico and Fromion, 2009;Chanson, 2004). These are nonlinear, hyperbolic Partial Differential Equations (PDEs) which are difficult to discretize and solve. The classical finite difference of finite volume methods are usually computationally expensive, hence less suitable for an online estimation application. Therefore, the use of a model reduction method is explored in this study, to obtain accurate and fast solutions for the shallow water equations. Orthogonal collocation method is used to reduce the PDE model into Ordinary Differential Equations (ODEs) efficiently, and the model results are validated using simulations. For the validation, a well balanced and well developed numerical scheme, which is specifically developed for the shallow water equations by Kurganov and Petrova (2007) is used. The use of collocation method in solving open channel problems has been studied previously by a lot of researchers including Georges et al. (2000); Dulhoste et al. (2004); Jinasena et al. (2017); Layton (2003), but these are limited only to prismatic channels. The Venturi is a non-prismatic channel which gives better fluid level differences along the channel length, therefore the model reduction for the non-prismatic channels is described in this study. Further, a sensitivity analysis is done to find the effect of different input, geometrical and fluid parameters on the model output.
The paper is organized as follows. The development of the mathematical model from the shallow water equations and the model reduction is described in detail in the Section 2. This is followed by the sensitivity analysis in Section 3. Then the simulation results and the discussion based on the results are stated in Section 4. Finally, the conclusions drawn by the results are summarized.

The Model
The model is the 1D shallow water equations for a non-prismatic channel, with the following assumptions (Chaudhry, 2008;Litrico and Fromion, 2009).
• The pressure distribution is hydrostatic.
• The channel bed slope is small i.e. the cosine of the angle it makes with the horizontal axis may be replaced by unity.
• The head losses in unsteady flow (due to the effect of boundary friction and turbulence) can be calculated through resistance laws analogous to those used for steady flow.
The Equations for a 1D, unsteady, non-prismatic, open channel system, is expressed as, where A(x, h, t) is the wetted cross sectional area normal to the flow, h(x, t) is the depth of flow, Q(x, t) is the volumetric flow rate and S f (Q, x, h) is the friction slope. I 1 , the first moment of area represents the hydrostatic pressure term and I 2 represents the pressure forces in the fluid volume, which occur from the longitudinal width and slope variations. g is the gravitational acceleration, t is the time and x is the distance along the flow direction (Chow, 1959;Chaudhry, 2008). β is known as the momentum correction coefficient or the Boussinesq coefficient and corresponds to the deviations of the local velocity over the mean velocity of the flow. The channel bed slope S b (x) is calculated by − ∂z ∂x (z is the absolute fluid level), where it is considered positive when sloping downwards. The friction slope S f is defined from the Gauckler-Manning-Strickler formulae as follows (Chow, 1959).
where n M is the Manning roughness coefficient and R is the hydraulic radius given by A P . Here, P is the wetted perimeter. For a channel with an isosceles trapezoidal cross section I 1 and I 2 can be found as follows, For a channel with a uniform side slope dSs dx would be zero. Therefore, the Eqs. 1 and 2 can further be simplified for an isosceles trapezoidal channel with uniform side slope and no lateral inflow rate. This simplified set of equations can be presented in a simpler form as follows, where U vector is the vector of conserved variables as follows, and S is the source term as follows, Then the mass flow rateṁ is calculated by multiplying the volumetric flow rate Q with a constant density ρ for a given drilling fluid.

Model Reduction for Non-Prismatic Channel
The shallow water equations are (as described earlier), a nonlinear hyperbolic system of PDEs. These are commonly known as very complicated systems to solve, because of the non-smooth solutions which could also contain shock and rarefaction waves, and the possible discontinuities (occurs due to discontinuous bottom topography or cross section) (Kurganov, 2018). Further, the numerical solutions could break down even for a system with smooth initial data and no discontinuities. Therefore, solving the shallow water equations in a stable and accurate manner, specially for a nonprismatic channel (with discontinuities) is a difficult task (Kurganov, 2018). There are well-balanced and well-developed numerical schemes available, which address these issues. However, due to the demand of high computational resources these schemes are not suitable for a real time estimation application. Therefore, a model reduction is applied here, in order to obtain fast and accurate enough, application oriented solutions. The model reduction is done by the use of orthogonal collocation method. The use of the orthogonal collocation method for shallow water equations for a prismatic channel is described in detail in . The construction of ordinary differential equations from orthogonal collocation on the spatial domain, for a non-prismatic channel is briefly described here.
The spatial length x of the channel is divided into n − 1 non-equidistant spaces between n points, which are known as collocation points (CP). The state vector U of the shallow water model at these specific points can be approximated by polynomial interpolation (Isaacson and Keller, 1966). In this study, the Lagrange interpolating polynomial is used for this purpose. The approximated states vector U a can be expressed as follows, Here, the subscript i denotes the i th position along the channel and L i is a weighted fraction defined as follows, where, x is the entire set of n number of positions. The derivatives of these approximated states can be derived as follows, where L ij is the element at i th row and j th column of the matrix L , Similarly, the other functions of x, βQ 2 A , I 1 , h and W can be approximated using the same method and the derivatives can also be found. The equations expanded from the original Eqs. 6 can be approximated with these approximated functions as follows, and Considering the functions at the selected CPs, the approximated value will be the same as the functional value, since the approximation of the function was done at these particular points. Therefore, the approximated Eqs. 14 and 15 would be equal to zero and can be expressed as follows after re-arranging, and The number of CPs in this study are selected as 2, 3 and 4, and the position of these points are selected according to the shifted Legendre polynomials . The corresponding matrices for L for 2, 3 and 4 CPs are stated below, respectively, where L is the length of the channel. The geometry of the channel that is used for the simulations is shown in Figure 2.

Kurganov-Petrova (KP) Scheme for Validation
To compare the results from the collocation model, a well developed, semi-discrete, second order and central upwind scheme, which has specifically been used to solve Saint-Venant systems (Bernstein et al., 2016;Bollermann et al., 2013) is used. With the use of the KP scheme by Kurganov and Petrova (2007), the PDE model is discretized in space into a set of ODEs. For a control volume as shown in Figure 3, the ODEs are written as follows, S is the average value of the source term calculated using the average values of the conserved variables. H j± 1 2 represent the fluxes flowing into the cell (minus sign) and out of the cell (plus sign) respectively. Assuming that there are no changes in the bed slope the fluxes are given by, can be calculated as the largest and smallest eigenvalues of the Jacobian of the system.
Here h d is the hydraulic depth and is equal to the crosssectional area of the flow divided by the top width of the flow (T ) that is exposed to the atmospheric pressure (h d = A/T ).

Bayesian Sensitivity Analysis
A sensitivity analysis is a generic method used to evaluate models. This is a way of decomposing the input uncertainty, where we can determine the parameters which are the most influential on the model output. This is done by analyzing the changes in the model output values that happen due to modest changes in the model input values. There are several techniques available for sensitivity analysis, where most approaches determine the effect of changes in one parameter when the others are not changing. However, the combined effect of multiple parameters can also be found. The sensitivities are determined as a percentage variation of the output for a percentage change of each input and/or parameters.
The errors and approximations in the input data measurements, parameter values, model structure and the numerical schemes (model solution algorithm) are the sources of uncertainty in this system. This model highly depends on the values of the parameters, therefore, only the uncertainty in parameters were studied and the analysis of the other sources of uncertainty is out of scope of this study. For this study, all the parameters were tested for sensitivity on the model output flow rate. All the geometrical parameters can be measured directly for a particular system. But the parameters which depend on the fluid properties (n M and β) cannot be measured always. Therefore a further analysis was done to quantify the effect of uncertainties in these two parameters on the model output and to understand the relationship between these two parameters and the model output.
A Bayesian analysis is used here. The Bayesian method will find the posterior probability with the use of a chosen prior and a likelihood. According to the Bayes theorem, The probability of Y given X (posterior p (Y |X)) is found by the probability of X given Y (likelihood p (X|Y )), the probability of Y (prior p (Y )) and the probability of X (model evidence p (X)) (Sivia and Skilling, 2006). Also the marginal likelihood can be obtained by marginalization over the parameters (for example X) as follows, The posterior probability of the output flow rate can be obtained using marginalization for the above mentioned two parameters as follows.  (28) Here,ṁ is the output of the model i.e. mass flow rate (volumetric flow rate Q × density ρ), that we are interested in. The model is denoted by M and the input data is denoted by D, which are considered to be noise free data. a, c, b and d are the minimum and maximum values for the two parameters, respectively. These values are chosen from test simulations, where the full range of values for which the model converges is taken. Further, with use of the product rule and the independence between n M and β, the right hand side can be expanded as follows, Here, n M and β are assumed to be random variables with Gaussian distributions of known means and standard deviations. The mean value is decided based on the test simulations, where it will give a certain desired flow rate as the output. The standard deviation is selected to cover the entire range of values that are possible for the simulated flow condition and model convergence. Further, both the data and the model are deterministic, therefore p (ṁ|n M , β, M, D) can be written as a delta function as follows, which is zero everywhere except at M (D, n M , β) =ṁ.

Results and Discussion
The models are simulated for the trapezoidal nonprismatic channel geometry using MATLAB (9.2). The initial conditions for the simulation are selected carefully. Runge-Kutta 4 th order time integrator is used with a fixed time step. The design of the channel geometry (see Figure 2) ensures trans-critical flow conditions for all conditions at the upstream. For example, if the upstream conditions are sub-critical, the flow becomes critical or super-critical at the throat. And if the upstream is at super-critical conditions, the flow at the throat becomes sub-critical. The trans-critical flows allow the possibility of using one or more boundary conditions. Hence, both of these trans-critical flow conditions are simulated with the inlet fluid level (converted in to the wetted cross sectional area) as the boundary condition for each different model. The different models are the model with KP scheme, and the models with different numbers of CPs. The results and the discussion for these followed by the sensitivity study results are stated below. The position of the collocation points along the channel were selected in a way that always the last point would lie in the middle of the throat, and the rest of the points are at upstream. Both the flow rate and the fluid levels were calculated at all points, except the boundary fluid level.

Comparison of Different Models
The collocation models were simulated for one boundary condition at the upstream, which is the wetted cross sectional area at the inlet.
Step responses were given to check the model responses for changes with time. Figure 4 shows the results obtained for the mass flow rates for the super-critical to sub-critical flow condition using each model. The mass flow rates are calculated from the volumetric flow rates for a given drilling fluid density. The results for the collocation models are compared with the KP scheme in Figure 4. It is worth mentioning that the KP scheme is difficult to implement for super-critical upstream conditions with fluid level as a boundary. Therefore, KP scheme is used with the upstream flow rate as the boundary and the resulting fluid level at upstream is used as the boundary for the collocation models. The mass flow rates for the super-critical to subcritical, trans-critical conditions are well comparable in all the models. These are relatively easier simulation conditions because the interference from waves (oscillations in fluid levels which propagates to the mass flow rate calculations) is less significant due to high velocities of the flow. However, more interesting flow behaviors can be observed with the trans-critical condition from sub-critical to super-critical, which is shown in Figure 5. Here, the KP scheme is also simulated taking the wetted cross sectional area at the upstream as the boundary.
From Figure 5, it can be clearly seen that the model with 3 and 4 CPs is more comparative with the KP results than the model with 2 CPs. Further, the results with 4 CP model is closer to the KP results than the results with 3 CP model. This is expected because the more the CPs, the better the approximation is. However, the number of oscillations that occur due to step changes is higher when the number of CPs are higher. But the time taken to reach steady state is shorter with the increase of number of CPs. Further, the amplitude of the oscillations are bigger, with 3 CP model. This is due to the Runge's phenomenon, which is the occurrence of oscillations when using polynomial interpolation over equi-spaced interpolation points. The 4 CPs are not equi-spaced, therefore, this effect is not seen with 4 CP model.
The shallow water equations, which are used in this study, is presented in the conservative form, so that it preserves the mass conservation along the entire spatial domain at steady state. This mass conservation can be seen in the collocation models, since the mass flow rates at all the CPs are the same. The 2 CP model was the quickest to simulate and the 3 and 4 CP models take about 10-15 times more simulation time than the 2 CP model. KP scheme usually takes 40-50 times more simulation time than the CP models. The results from the 2 CP model for subcritical upstream conditions is unsatisfactory. However, the use of 3 or 4 CPs for any flow condition, is enough to obtain similar results to KP scheme. Therefore, in an online estimation application the use of 3 or 4 CP model would give sufficiently accurate results with less computational time. This means that with the 3 or 4 CPs model it is possible to obtain a fast estimate of the flow rate in real-time. Although there are oscillations in the presence of an abrupt step change in a sub-critical flow condition will eventually reach to the steady state in a relatively short period of time.
As a summary, the reduced order model was able to capture the general fluid flow behaviors. Appearance of transient oscillations can be observed in both the models (reduced order model as well as the KP scheme), when a step change is given to the input. The CP models (especially the 2 CP model) show prolonged oscillations before settling down to a steady state value in the presence of sharp step changes in the input/inlet boundary conditions (at t = 2500s in Figure 4). Such abrupt step changes in the input is simulated here as this can be expected to occur during a kick or loss in drilling system. When the stepping is done gradually throughout a longer time period (ramping, at t = 2000s in Figure 4), the flow comes to a steady state value quickly, with less oscillations. This behav-ior is expected during the pipe connection procedure in drilling, where the mud pump is stopped and started gradually.

Sensitivity Analysis Results
The sensitivity analysis was done for all the models (KP scheme and collocation with 2, 3 and 4 CPs) for sub-critical upstream flow conditions. The inlet fluid level and all of the fluid and geometrical parameters are tested for sensitivity of the flow rate.
A spider plot, which shows the sensitivity of the input parameters on the output flow rate for each model is shown in Figure 6. For all the models, the channel bed slope has the highest sensitivity. In addition, the inlet fluid level and the bottom widths of the channel are also sensitive parameters. However, the calculations are done using the nominal fluid levels for bed slope changes. But in reality, when the channel bed slope changes, the fluid levels in the channel also change accordingly, thus may reduce the sensitivity to the flow rate. This is also true for the bottom width W 1 . The change of fluid level due to the changes in the bed slope and the bottom width W 1 , are not taken into account in the calculations, for simplicity of comparison. Nonetheless, the effect of the channel bed slope should be taken into consideration in a real fluid flow scenario due to its high sensitivity with the flow rate and the fluid flow conditions. A small increase in the channel bed slope of the Venturi channel that is taken into consideration can change the subcritical upstream conditions into super-critical conditions. Therefore, in this study, the bed slope is changed only within a small range, due to the limitation of the flow conditions based on the used geometry.
Most of the parameters have a bigger impact on the 2 CP model than on the other models.
Even though the effect of geometrical parameters (W 1 , W t , S s and L) on the 3 and 4 CP models are comparatively lower, the flow rate can still be affected significantly. In order to reduce the effect on the magnitude of the flow rate, the geometrical parameters can be measured accurately for each geometry. The fluid properties such as friction factor n M and β are tunable parameters, hence can be estimated together with the flow rate, in the context of adaptive estimation. However, the input (Ain) and the channel bed slope may change with time, hence have to be measured with a good accuracy and precision due to their high sensitivity.
For an estimation application, an analysis of the sensitivity of the estimable parameters on the output flow rate will be informative in advance. Therefore, the change of flow rate due to perturbations in the two estimable parameters (n M and β) is analyzed with the use of the Bayesian sensitivity analysis described in Section 3. In order to determine the marginal posterior distribution of the mass flow rate shown by Eq. 30, a large number of independent samples (10000 pairs), representing n M and β, were generated from their respective Gaussian distributions. Then for each pair, depending on the model, the corresponding steady state mass flow rates were calculated. The histogram of n M , β andṁ can be considered as an estimate for the joint posterior distribution of these triples, conditional on the appropriate model (p(n M , β,ṁ|M )). Similarly, the marginal posterior distribution of the mass flow rate (p(ṁ|M )) can be approximated by considering only the histogram of the mass flow rate of the samples. These posterior distributions for each model are shown in Figure 7. For the sake of comparison, the histograms are vertically scaled so that their highest value is unity. Further, a Gaussian distribution is fitted to the marginal posterior histogram of the output flow rate (shown by red curve in Figure 7). The two parameters for the models had to be tuned to ensure the model convergence for sub-critical to super-critical flow conditions. Therefore, in some of the models (especially in 3 and 4 CP models) the parameter values may not be meaningful. It is shown that the used Gauckler-Manning-Strickler formulae with a constant n M is suitable only for uniform flows (Chow, 1959). During the transition the flow is not uniform. Further, it is suggested that the n M depends on the hydraulic radius and the Froude number, hence could not be a constant (Tullis, 2012). However, it is also shown that at super-critical flow conditions this is relatively constant (Tullis, 2012). This might be the reason that the collocation models had to be tuned for sub-critical flow conditions.
The mean (µ) and the standard deviation (σ) of the estimated mass flow rate distributions together with the values of the chosen parameter distributions for each model are tabulated in Table 1. For comparison, the percentage deviation is shown in Table 2. The full range of possible parameter values for n M and β are not  Because of this the distribution could be truncated towards the tails, even though it is approximated with the Gaussian distribution. As shown in the tables, the 2 CP model has a significant change in flow rate, due to the high sensitivity of the parameters. Further, it can be seen from the tabulated percentage changes, that the perturbations in the two parameters have caused less impact on the flow rate, when the number of CPs are 3 or 4. Since there is no significant improvement with the 4 CP model than the 3 CP model, it is a trade off between the slightly better accuracy and the slightly higher computation time.
However, the dominance of the correlation of the two parameters with the output flow rate will strictly depend on the choice of these two parameter values, which then indirectly affects the output flow rate itself. These correlations between the parameters and the flow rate can clearly be seen from the scatter plots obtained from the simulations. The scatter plots and the contours drawn on their densities are shown in Figure 8, for all the models. The 3 and 4 CP models show a dominant correlation of flow rate with n M , while the KP model shows more dominance with β. This could be due to the choice of the mean values of the parameters for a desired output flow rate. Since the parameters are tuned in order to get a desired flow rate, we cannot decide about the accuracy between the models based on the parameter values. However, we can have an idea about the desired flow rate and the precision. Further, the 3 CP results with β (in Figure 8b) are truncated due to the crossing into the sub-critical region at the throat as oppose to the super-critical region. This further indicates that a careful selection of the values for these two parameters by means of estimation/adaptation will be needed for an accurate flow rate estimation. This study will be carried out in the future.
Further, the expected relationship of both n M and β with the flow rate is generally a negative relationship, based on the model structure. This cannot be seen with n M and the flow rate with the 2 CP model as shown in Figure 8, in fact it shows a complete opposite relationship, which is difficult to interpret according to physical laws. Since the 2 CP model is highly sensitive to most of the parameters, this might have occurred due to improper domination of other parameters over n M . With 2 CP model the CPs fall on the boundaries, and the information of the states inside the channel between the boundaries is not captured by the model. 2 CP model is faster and simpler but less accurate and less physically sound than the other models. This further implies that the model simplification comes with a cost on accuracy and physical explanation. Although these models are not suitable for estimating the friction factor, it is (3 CP or 4 CP model) accurate and fast enough to estimate the flow rate based on the level measurement.
However, to obtain a better understanding of the model, other uncertainty sources such as the structural uncertainty and numerical uncertainty may also need to be explored. The structural uncertainty which occurs due to the improper approximations of the real system and the numerical uncertainty which occurs due to the discretization and model reduction could be studied in future with a comparison to a real system. Further, a proper selection of the tunable parameters (n M and β) should be done in order to ensure a good estimation of the flow rate. Due to the high sensitivity of the bed slope and the inlet fluid level in all the models, a careful consideration of the uncertainty in these two measurements should be taken into account in a real world application. The knowledge and results obtained from this study will be beneficial for such studies and applications in the future.

Conclusions
This paper highlights a reduced order open channel model which can be used to estimate the flow rate in an open channel. The model is developed to be used for a Venturi channel in the drain back flow line during an oil well drilling to estimate the flow rate of the return fluid. The simplified model is simulated for a non-prismatic open channel Venturi with a throat. The reduced order model is a set of ODEs hence faster and less complex than the shallow water equations. The reduced order models are validated using simulations from a well-developed finite volume method. Two trans-critical conditions (super-critical to sub-critical and sub-critical to super-critical) are simulated using the upstream fluid level as the boundary condition. Three or four collocation points seem to be sufficient to obtain a good accuracy for the determination of the flow rate in both the flow conditions. Further, a parameter sensitivity analysis of all the models is also done to evaluate the reduced order models. This model structure is well suited to be used in process control and state estimation algorithms, where the state estimates have to be computed in real-time. The proposed solution of using a Venturi channel for on-line estimation of return flow rate during drilling seems possible together with the use of a reduced order collocation model.