EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6181 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Simulation of Cutting Transport in a Drilling Operation by Using TVD RKDG Scheme with WENO Limiters Ilyas Khan1,5,6, Asad Rehman2,∗, Shahid Mehmood3, Saqib Zia3,∗, Wei Sin Koh4 1 Department of Mathematical Sciences, Saveetha School of Engineering, SIMATS, Chennai, Tamil Nadu, India 2 Department of Information Technology, Pothohar Campus Gujar Khan, University of Punjab, Pakistan 3 Department of Mathematics, COMSATS University Islamabad, Park Road Chak Shehzad, Islamabad, Pakistan 4 INTI International University, Persiaran Perdana BBN Putra Nilai, 71800 Nilai, Negeri Sembilan, Malaysia 5 Hourani Center for Applied Scientific Research, Al-Ahliyya Amman University, Amman, Jordan 6 Department of Mathematics, College of Science Al-Zulfi, Majmaah University, Al-Majmaah 11952, Saudi Arabia Abstract. In this article, we present a numerical investigation of two-phase flow involving liquid and solid particles (cuttings) during drilling operations in wellbores. The dynamics of this complex system are simulated using a high-resolution second-order Total Variation Diminishing (TVD) Runge-Kutta Discontinuous Galerkin (RKDG) method, enhanced with a Weighted Essentially Non-Oscillatory (WENO) limiter. Unlike traditional slope-limiter techniques, the integration of WENO methodology with the RKDG scheme significantly improves the accuracy and stability of the numerical solution by effectively capturing sharp interfaces and discontinuities without spurious oscillations. A two-phase drift-flux model is employed to account for the relative motion between the liquid and solid phases, enabling realistic simulation of the transient and dynamic behavior of cuttings transport under various well conditions. Through a series of test cases, we analyze the influence of varying rates of penetration (ROP) on critical parameters such as bottom hole pressure (BHP), cuttings concentration, and pressure evolution over time. The results demonstrate that the proposed scheme offers enhanced accuracy and computational efficiency compared to conventional approaches. 2020 Mathematics Subject Classifications: 65Nxx, 65M08, 65Z05 ∗Corresponding author. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6181 Email addresses: i.said@mu.edu.sa (I. Khan), assad013@gmail.com (A. Rehman), sm0694921@gmail.com (S. Mehmood), saqibzia81@hotmail.com (S. Zia), weisin.koh@newinti.edu.my (W. S. Koh) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 2 of 21 Key Words and Phrases: Two phase flow, rkdg, drift flux flow model, shock wave, source terms, smart microgrids 1. Introduction The transportation of cuttings in wells has been the subject of experimental investi- gation by numerous scientist since the second last decade of 20th centuary. Tomren et al. [1] were among the pioneering researchers to study cuttings transport in directional wells. The extraction of drill cuttings from the bore-hole is a pivotal step during drilling operations, especially when working on horizontal, extended-reach, or slanting wells. Lack of proper bore-hole cleaning can result in substantial issues, such as pipe obstructions, reduced drilling velocity, increased bit wear, fluid loss, and excessive torque and drag. A well-designed drilling mud and an appropriate circulation system are essential for ensur- ing the extraction of cuttings from the bore-hole. The success of the drilling operation is greatly impacted by a number of critical elements, such as the type of formation, wellbore inclination, cutting size, penetration rate, and eccentricity of the drill pipe. It is difficult to create an ideal hole cleaning plan that considers each of these aspects [2–5]. Two phase flow models are critical for predicting fluid behaviours in many different industrial applications, including power production, cooling systems, and the extraction of gas and oil. These models hold significant value for multiple objectives, including gain- ing insights into fluid dynamics within pipelines, heat exchangers, and other machinery; facilitating the design and optimization of industrial processes; improving efficiency and reducing operational costs; and mitigating safety threats such as leaks and explosions[6–8]. To explore the attributes of two phase flows, scientists have devised various models and an extensive array of numerical methods. However, regardless of significance, these models can pose numerical challenges because of the complexities arising from the interface between the both phases. Specifically, these models present various numerical obstacles, such as discontinuous solutions, shock formation, and loss of hyperbolicity. To address these numerical challenges, drift flux model was introduced by Zuber and Findly [9]. Subsequently, several authors have made improvements and modifications to the drift flux flow model, as demonstrated in articles [10–14]. In transient-state and steady-state models the primary distinction is lies in the representation of the equations of motion. The rate of change of velocity, density and temperature involve in the transient equations of motion. A pipeline flow is typically characterized as a transient process due to the changing inlet and outlet flows. Moreover, variations in ambient conditions, control set points, and temperature further contribute to this transient behavior. Therefore, having a transient flow model that accurately captures the dynamics occurring in the system proves advantageous [15]. In this article, RKDG is developed to simulate the cuttings flow in well operations by considering the two phase drift flux model. This method combines the best features of finite element methods (such as weak formulation and test function spaces) and finite volume methods (such as nonlinear limiters, numerical fluxes and smart microgrids). First I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 3 of 21 of all, the discontinuous Galerkin numerical method was developed by Reed and Hill [16] to solve the hyperbolic equations. After that this numerical scheme was extended and employed to solve parabolic, elliptic and hyperbolic types of partial differential equations [17–19]. This scheme was further modified by Cockburn and his team in [20–23]. The explicit time discretization techniques for total variation diminishing (TVD) were discussed in [24–27]. The proposed numerical scheme maintains the property of total variation bounded in the means. This scheme is capable to handle sharp shock transitions, without introducing oscillations. Further, the artificial viscosity is introduced by the numerical flux through upwinding to maintain stability in the scheme when employing piecewise-constant approximations, for detail see [20, 21, 24, 25]. This work presents a refined numerical approach for simulating two-phase drilling hydraulics by coupling the Runge-Kutta Discontinuous Galerkin (RKDG) method with WENO limiters within a drift-flux framework.To the best of our knowledge, only the Advection Upstream Splitting Method – Velocity (AUSMV) numerical scheme has been employed previously to simulate cuttings transport in well operations [28]. The incor- poration of WENO limiters [29–32] into the time-stepping algorithm enhances stability and effectively suppresses spurious oscillations. The resulting scheme accurately resolves shock phenomena and maintains sharp solution profiles with high fidelity. Through various test cases, the influence of changing rates of penetration (ROP) on bottom hole pressure (BHP), cuttings concentration, and pressure evolution is systematically analyzed. The results confirm that the proposed method outperforms conventional schemes in both ac- curacy and computational efficiency. Future studies may focus on extending the current model to incorporate non-Newtonian fluid effects and thermal variations, which are criti- cal in real-world drilling environments. Additionally, the proposed numerical scheme can be extended to compute the solutions of nonlinear partial differential equations such as [33–36]. The rest of article is structured as follows. In section 2, the mathematical form of two phase drift flux model and eigen structure is presented. In section 3, the TVD RKDG is derived for the proposed mode. In section 4, different cases are considered to analyze the effects of cuttings on bottom hole pressure and observe the effects of different parameters in the process of cutting flows. Finally, the conclusions are drawn in section 5. 2. Drift flux model for two phase flows In this section, we present the mathematical representation of the two phase drift flux model [9] and provide an explanation of its eigenvalues. The model is based on several key assumptions to simplify the complex physics of multiphase transport. Firstly, both the liquid and solid phases are assumed to share a common pressure field. Assuming that the cuttings are sufficiently small and uniformly distributed. The relative motion between phases is captured using empirical drift velocity correlations. Additionally, the model assumes isothermal conditions and neglects thermal effects on fluid or particle behavior. Drift flux model consists of three partial differential equation, two of them are liquid and cutting phase mass conservation equations and one mixture momentum equations are I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 4 of 21 given here ∂ ∂t (alρl) + ∂ ∂x (alρlul) = 0, ∂ ∂t (agρg) + ∂ ∂x (agρgug) = 0, ∂ ∂t (alρlul + agρgug) + ∂ ∂x (alρlu 2 l + agρgu 2 g + p) = s, (1) where ag and al represent the volume fractions of cutting and liquid, ug and ul denote the cutting and liquid velocities, ρg and ρl represent cutting and liquid densities, respectively. The cuttings density was assumed fixed to be 2500kg/m3. The source term is denoted by s. The mutual pressure is represented by p . s = Ff +Gf = −32umixµmix D2 + (agρg + alρl)g, (2) The diameter of pipeline from the inner side is denoted by D and g = 9.81 m/s2. Viscosity term µmix and the average velocity term umix in the frictional forces Ff are µmix = alµl + agµg and umix = agug + alul. (3) There are seven unknown variables are involved in drift flux model al, ag, ρl, ρg, ul, ug, and p. Consequently, to obtain a solution, it is necessary to establish four additional relationships, which are as follows: ρl = ρl,0 + p− pl,0 α2 l , (4) The pressure is calculated by the following relation p = pl,0 + ( alρl (1− ag) − ρl,0 ) · a2l , (5) al + ag = 1. (6) and ug = Coumix + ud, (7) here,αl represent the speed of sound in liquid phase. Additionally, ρl,0 and pl,0 are the reference liquid density and reference liquid pressure, while Co and ud represent the drift cutting velocity and a flow-dependent parameter. The specific values for these parameters are given in Table 1. 2.1. Eigenvalues calculation The drift flux model (1) is represented in the conservation laws form ∂tW+ ∂xF(W) = S(W), (8) where I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 5 of 21 Table 1: Parameter values. Description Parameters Value cutting phase viscosity µg 2.10−2 Pa.s liquid phase viscosity µl 5.10−2 Pa.s velocity of sound in cutting phase αg 1500 m/second velocity of sound in liquid phase αl 1000 m/second reference density liquid ρl,0 1000 Kg/m3 reference pressure liquid pl,0 105 Pa choke pressure pchoke 100 bars W =  alρl agρg alρlul + agρgug  , F(W) =  alρlul agρgug alρlu 2 l + agρgu 2 g + p  ,S(W) =  0 0 s , with, S(W) = (s1, s2, s3) = (0, 0, s), F(W) = (f1, f2, f3) = (alρlul, agρgug, alρlu 2 l + agρgu 2 g + p) and W = (w1, w2, w3) = (alρl, agρg, alρlul + agρgug). Now, the Eq. (8) can be expressed as ∂t  w1 w2 w3 + ∂x  w1wl w2wg w1w 2 l + w2w 2 g + p(w1, w2)  =  0 0 s  . (9) We have observed that the entire flux vector cannot be entirely expressed using the con- served variables in this system, making it impossible to precisely determine the sound velocity. The papers [37, 38] provide comprehensive explanations for approximation of the sound velocity. Although we have omitted the procedural details here, we will only consider the eigenvalues. Eigenvalues for the considered model is, λ1 = u− α, λ2 = u, (10) λ3 = u+ α. Where, α = √ p agρl(1− Coag) . (11) 3. Construction of TVD RKDG scheme for the drift flux model In this section, we obtain the two phase drift flux flow model using the RKDG scheme. This model is written in compact form as ∂tW+ ∂xF(W) = S(W), t > 0, x ∈ Ω. (12) I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 6 of 21 with initial conditions W(x, 0) = W0(x). (13) We subdivide the domain Ω into cells Ci = [xi− 1 2 , xi+ 1 2 ], where i = 1, ..., N. We denote the center of ith cell by xi and defined as 0.5 ( xi− 1 2 + xi+ 1 2 ) . The size of ith cell is defined as △xi = xi+ 1 2 − xi− 1 2 , for all i. Multiply the Eq. (12) by a smooth function τ(x) and integrate the resulting equation over the ith cell Ci. By integrating, we obtain∫ Ci ∂tW(x, t)τ(x)dx− ∫ Ci F(W(x, t))∂xτ(x)dx + F(W(xi+ 1 2 , t))τ(xi+ 1 2 )− F(W(xi− 1 2 , t))τ(xi− 1 2 ) = ∫ Ci S(W(x, t))∂xτ(x)dx, (14) and ∫ Ci W(x, 0)τ(x)dx = ∫ Ci W0(x)τ(x)dx. (15) Consider for each time t ∈ [0, tfinal], the solution approximation Wh of W belongs to discontinuous Galerkin finite element space Vh = {υ : υ|Ci ∈ Pk(Ci), i = 0, 1, 2...N}. (16) The Pk(Ci) is the polynomials space of degree at most k in cell Ci. We choose the local orthogonal basis {φl i, l = 0, 1, 2, ..., k} , named as scaled Legendre polynomials Pl(x), over cell Ci as defined below φl i(x) = Pl((x− xi)/△xi), l = 0, 1, 2, ...k. (17) Then, the approximate solution Wh(x, t) ∈ Vh is defined as Wh(x, t) = k∑ l=0 Wl i(t)φ l i(x), for x ∈ Ci, (18) where, Wl i(t) = 2l+1 △xi ∫ Ci Wh(x, t)φ l i(x)dx by using the orthogonal property∫ Ci φl iφ m i dx = { 0, l ̸= m △xi 2l+1 , l = m. Now, replace a smooth function τ(x) by the test function φl(x) ∈ Vh and flux function F(W(xi+ 1 2 , t)) by numerical flux F̂i+ 1 2 in Eqs. (14) and (15) to obtain a numerical scheme as follow∫ Ci ∂tWh(x, t)φl(x)dx− ∫ Ci F(Wh(x, t))∂xφl(x)dx+ F̂i+ 1 2 φl(x − i+ 1 2 )− F̂i− 1 2 φl(x + i− 1 2 ) = I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 7 of 21∫ Ci S(Wh(x, t))φl(x)dx, (19) with ∫ Ci Wh(x, 0)φl(x)dx = ∫ Ci W0(x)φl(x)dx. (20) The flux F̂i+ 1 2 = 𭟋 ( Wh(x − i+ 1 2 , t),Wh(x + i+ 1 2 , t) ) is any monotone numerical flux. Here, we use the Local Lax-Friedrichs as a monotone numerical flux which is defined below 𭟋 ( W− i+ 1 2 ,W+ i+ 1 2 ) = 1 2 ( F(W− i+ 1 2 ) + F(W+ i+ 1 2 )− ϑ(W+ i+ 1 2 −W− i+ 1 2 ) ) , (21) where ϑ = maxW |F′(W)| and Wh(x ± i+ 1 2 , t) are the limited values of approximate solution Wh(x, t) at the cell interface xi+ 1 2 . By using above definitions, the Eqs. (19) and (20) simplify to dWl i(t) dt = −2l + 1 △xi (∫ Ci F(Wh(x, t))∂xφl(x)dx+ F̂i+ 1 2 φl(x − i+ 1 2 )− F̂i− 1 2 φl(x + i− 1 2 )+∫ Ci S(Wh(x, t))φl(x)dx ) , (22) and Wl i(0) = 2l + 1 △xi ∫ Ci W0(x)φ l i(x)dx. (23) Next, to obtain non-oscillatory discontinuity propagation and high order accuracy by RKDG method, the limiter procedure is divided into two parts. Firstly, identify the cells, named as “troubled cells”, that may need the limiting procedure. Secondly, reconstruct the polynomial solutions in these troubled cells by using WENO reconstruction. Here, the troubled cells are identified by total variation bounded (TVB) limiter [39], which is described as Wh(x + i+ 1 2 ) = W (0) i +Wi, Wh(x + i− 1 2 ) = W (0) i −Wi, (24) where Wi = k∑ l=1 Wl iφ l(xi+ 1 2 ), Wi = − k∑ l=1 Wl iφ l(xi− 1 2 ). Now Wi and Wi are modified by the modified minmod (MM) function, as follow W mod i = MM(Wi,W (0) i+1 −W (0) i , W (0) i −W (0) i−1), I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 8 of 21 W mod i = MM(Wi,W (0) i+1 −W (0) i , W (0) i −W (0) i−1), where MM is defined as MM(a1, a2, a3) = { a1, if |a1| ≤ M(△x)2, M(a1, a2, a3), otherwise, where M is a positive constant and the minmod function M is defined as M(a1, a2, a3) = { S ·min1≤i≤n |ai|, if sign(a1) = sign(a2) = sign(a3) = S, 0, otherwise. Note that, the cell Ci is called a troubled cell when W mod i ̸= Wi and W mod i ̸= Wi. For the troubled cell, the moments Wl i are reconstructed by the WENO reconstruction which retains the cell average W (0) i . Note that, for the (k + 1)-th order accurate Pk based DG scheme, we use the Gauss Lobatto quadrature rule which is accurate at least up to the order (2k + 2) and the order of accuracy for the WENO reconstruction must be at least (2k+ 1). To achieve this, (2k+ 1) neighboring cells Ci−k, ..., Ci+k are used for reconstruction of the point values of W at the Gauss Lobatto quadrature points. Next, for the reconstruction of WENO limiter for RKDG scheme, we will follow the same approach as described in [30] and obtained the WENO approximation to the moments as follows Wl i = 2l + 1 △xi k∑ m=0 ωm ∫ Ci pm(x)φl i(x)dx l = 1, ..., k. (25) Here, pm(x) denotes the k-th degree reconstructed polynomial associated with each of the stencils τm = ∪k l=0Ci+m−l, m = 0, ..., k and ωm are the nonlinear weights defined as ωm = ωm∑ m ωm , ωm = γm (ε+Bm)2 , (26) where Bm is the smooth indicator for each cell Ci which measure the smoothness of function pm(x) in the cell Ci. The Bm is defined as follow Bm = ∑ m ∫ Ci △x2l−1 i ( ∂l ∂xl pm(x) )2 dx. (27) The remaining two terms γm and ε in Eq. (26) respectively denote the linear weights and small constant which is usually taken as 10−6. Further, the linear weights γm should satisfy the following relation∫ Ci q(x)φl i(x)dx = k∑ m=0 γm ∫ Ci pm(x)φl i(x)dx, l = 1, ..., k, (28) I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 9 of 21 where q(x) denotes the 2k-th degree reconstructed polynomial, associated with the larger stencil T = ∪k m=0τm [29]. For time discretization, we will use second order TVD Runge- Kutta method [40]. The system of ordinary differential equations (22) can be written in the compact form as follows, dWh dt = Lh(t,Wh) . (29) To approximate equation (29), one may used r-order TVD RK, W (k) h = k−1∑ l=0 [ αklW (l) h + βkl∆tLh(t n + dl∆t,W (l) h ) ] , k = 1, 2, · · · , r , (30) by using the boundary conditions W (0) h = Wn h , W (r) h = Wn+1 h . (31) The coefficients for the second-order TVD RKM are provided as follows:[20] β10 = α10 = 1 , β21 = α20 = α21 = 1 2 , d0 = 0 , d1 = 1, β20 = 0. (32) While, the coefficients for the third-order TVD RK are provided as follows β10 =α10 = 1 , β20 = 0 , α20 = 3 4 , β21 = α21 = 1 4 , α30 = 1 3 , β31 = β30 = α31 = 0 , β32 = α32 = 2 3 , d2 = 1 2 , d1 = 1 , d0 = 0 . (33) The following condition ∆t ≤ ( 1 2k + 1 ) min(∆xj) λmax , (34) is known as Courant-Friedrichs-Lewy (CFL) condition [20], where ∆t, k, ∆xj , min(∆xj) and λmax represent time step, degree of polynomial used in DG scheme, spacial cell size, the most restrictive condition over the grid and maximum speed. The CFL condition is used in the time step scheme for ensuring the stability and convergence of the numerical scheme. Here, k = 1, 2 is the order of scheme which is second and third order respectively, and λmax reprsents the maximum eigenvalue. The adaptive nature of this time step results in a decrease for scenarios involving substantial variations (large slopes) in the solution and an increase otherwise. I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 10 of 21 4. Numerical Simulation In this section, five different test cases [28] are considered to analyze the efficiency of and proposed numerical scheme.More precisely, the proposed numerical scheme is used to examine the influence of cuttings on bottom-hole pressure, observes the variations of BHP over time, and explores the impact of changing penetration rates on BHP throughout the time duration. We consider the consistent geometry across all given cases, the depth of well is taken as 2000 m and inner diameter of the well is 5 inches, while the outer diameter is 8.5 inches. For complete description, the reader is referred to [28]. Case 1: A no-slip model is considered in the this case and over various time intervals the changes in cutting concentrations are observed. Before pumping starts, the well is completely full with liquid, and the flow rates are 22 kg/s with a cutting flux of 2 kg/s at the bottom of the well. In Figure 1, at 500 s, it is evident that well lower part contains cuttings, indicating the transportation of cuttings out of the well. However, initially, the cutting concentration is zero up to a depth of 1250 m, and then it steadily increases, eventually approaching a value of 0.04, as shown in the figure. Moving on to Figure 2, the cutting has progressed higher up in the well at 1500 s. It is also observed that the concentration values are constant, indicating a stable state, specifically from 1400 m to 2000 m. In contrast to Figures 1 and 2, Figure 3 shows that the steady state of the cutting concentration is achieved at a shorter well depth at a depth of 2000 m. Constantly cutting concentration values are observed from a depth of 1150m onward. After 3000 s, the steady state is reached in Figure 4, and the cutting concentrations are consistently seen over the whole well. Figure 5 displays the final steady state conditions when the final time is set to 4000 s. The cutting concentration is measured at 0.04, indicating that the cuttings are distributed throughout the entire well. Figure 6 shows a graph of bottom hole Pressure over time in the well, highlighting how the cutting rate impacts BHP. When we first begin pumping cuttings into the well, it takes around 200 s for BHP to start rising. This rise is due to the friction within the well caused by the start of the flow rate. As the well continues to fill with cuttings, the BHP steadily increases. Case 2: In this case, we introduce a change in the mass rate of cutting, which is now 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 1: Case 1: t= 500 s, mass rate = 2 kg/s, plotting of well depth vs cuttings concentration. 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 2: Case 1: t = 1500 s, mass rate = 2 kg/s, plotting of well depth vs cuttings concentration. I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 11 of 21 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 3: Case 1: t= 2000 s, mass rate = 2 kg/s, plotting of well depth vs cuttings concentration. 0.028 0.03 0.032 0.034 0.036 0.038 0.04 0.042 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 4: Case 1: t= 3000s, mass rate = 2 kg/s, plotting of well depth vs cuttings concentration. 0.034 0.035 0.036 0.037 0.038 0.039 0.04 0.041 0.042 0.043 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 5: Case 1: t= 4000 s, mass rate = 2 kg/s, plotting of well depth vs cuttings concentration. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 192 194 196 198 200 202 204 206 208 210 212 B ot to m ho le P re ss ur e [b ar ] Figure 6: Case 1: t= 4000 s, plotting of bottom hole pressure vs time. set at 4 kg/s. We observe variations in the concentrations of the cutting for different time intervals, while keeping all other variables the same as in Case 1 . Figure 7 demonstrates that when using a higher cutting rate, the cutting concentrations are also higher, reaching a value of 0.077. This concentration is higher compared to the previous case. Similar to Case 1, the cutting concentrations initially start at zero up to a depth of 1250 m, and then steadily increase, approaching the value of 0.077. In Figure 8, we observe that the concentration of the cutting gradually increases from a depth of 400 m to 1400 m. At a level of 0.77 concentration, it reaches a steady state. This indicates that the cutting is distributed primarily in the lower part of the well. According to Figure 9, the cutting concentration reaches its steady state at a shorter hole depth at 2000 s when compared to Figures 7 and 8. At around 1150 m down, the concentration reaches its steady state. Figure 10 depicts that the cutting concentration can be distributed throughout the entire well relatively quickly. Constant cutting concentrations are observed from a depth of 250 m to 2000 m. At 4000 s, constant cutting concentrations of 0.077 are shown over the whole well in figure 11. The well’s bottom hole pressure against time plot is shown in Figure 12, which illustrates how the cutting rate impacts BHP. Like Case 1, an initial I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 12 of 21 increase in BHP occurs when pumping cuttings and liquid begins, attributed to friction. Gradual increases in BHP are observed due to the hydrostatic pressure as the cuttings are moved upward in the well. It takes around 3500 s to reach a steady condition. When the injection rates are zeroed out at 4000s, well friction goes away. The difference between the simulations in both cases is mainly due to the content of cuttings in the well. The steady state pressure observed during circulation is higher in this case, which can be attributed to the higher cutting rate and resulting higher cutting concentration. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 7: Case 2: t= 500 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concentration. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 8: Case 2: t= 1500 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concentration. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 9: Case 2: t= 2000 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concentration. 0.055 0.06 0.065 0.07 0.075 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 10: Case 2: t= 3000 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concentration. Case 3: In this case, we will simulate with a cuttings mass rate of 6 kg/s, while the end times change according to time intervals. The cuttings have penetrated a depth of 1250 m at 500s, as seen in Figure 13. However, as a result of the increased cutting rate of 6 kg/s, the concentration of cuttings is roughly 0.115 at the conclusion of the simulation, which is greater than in the prior situations. In Figure 14, the concentration of the cuttings rises steadily from 400 m to 1400 m, after which it steady state at 0.115, much as in the other instances. This points to the bottom of the well as the primary distribution point for the cuttings. The steady state of the cutting concentration is achieved at a shorter I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 13 of 21 0.066 0.068 0.07 0.072 0.074 0.076 0.078 0.08 0.082 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 11: Case 2: t= 4000 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concentration. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 190 195 200 205 210 215 220 225 B ot to m ho le P re ss ur e [b ar ] Figure 12: Case 2: t= 4000 s, mass rate = 4 kg/s, plotting of bottom hole pressure vs time. well depth compared to Figures 13 and 14, as shown in figure 15, when the end time is raised to 2000 s. At around 1150 m down, the concentration reaches its steady state. At 3000 s, as shown in Figure 16, the cuttings have been distributed throughout the entire well. Constant cutting concentrations can be observed from a depth of 250 m to 2000 m. In Figure 17, at 4000 s, the cuttings concentration is the same across the entire well, indicating the achievement of final steady state conditions. From 2 kg/s to 4 kg/s and then 6 kg/s, the cuttings concentration grew in accordance with the cuttings mass rate. In Figure 18, we can see the well’s bottomhole Pressure (BHP) plotted against time. As in earlier instances, the first rise in BHP at 200 s is caused by the friction that occurs when cuttings and liquid are pumped. As the cuttings ascend in the well, the hydrostatic pressure rises, and the BHP rises with it. A steady condition is achieved after 3500 s. When the injection rates are turned off after 4000 s, the only pressure that remains is hydrostatic pressure, as well friction goes away. Here, we see stronger pressure than in the prior instances since the cuttings concentration is higher because to the higher cuttings rate. 0 0.02 0.04 0.06 0.08 0.1 0.12 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 13: Case 3: t = 500 s, mass rate = 6 kg/s, plotting of well depth vs cuttings concentration. 0 0.02 0.04 0.06 0.08 0.1 0.12 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 14: Case 3: t = 1500 s, mass rate = 6 kg/s, plotting of well depth vs cuttings concen- tration. I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 14 of 21 0 0.02 0.04 0.06 0.08 0.1 0.12 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 15: Case 3: t = 2000 s, mass rate = 6 kg/s, plotting of well depth vs cuttings concen- tration. 0.085 0.09 0.095 0.1 0.105 0.11 0.115 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 16: Case 3: t = 3000 s, mass rate = 6 kg/s, plotting of well depth vs cuttings concen- tration. 0.098 0.1 0.102 0.104 0.106 0.108 0.11 0.112 0.114 0.116 0.118 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 17: Case 3: t = 4000 s, mass rate = 6 kg/s, plotting of well depth vs cuttings concen- tration. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 190 195 200 205 210 215 220 225 230 235 B ot to m ho le P re ss ur e [b ar ] Figure 18: Case 3: t = 4000 s, mass rate = 6 kg/s, plotting of bottom hole pressure vs time. Case 4: In this test case, we consider two different rate of penetrations: an ROP of 10 m/h and an ROP of 30 m/h, which correspond to mass rates of 0.166 kg/s and 0.498 kg/s, respectively. Figure 19 illustrates that the cutting concentrations increase gradually until reaching a constant value of 3.55 × 10−3. The cuttings front appears to be smeared out, which can be attributed to numerical diffusion. In Figure 20, we observe that near this constant value, the cutting concentrations remain steady over a time period of 4000 s. At a ROP of 10 m/h, the simulation results are shown in Figure 21. After 200 s of inactivity, the flow rates are activated, leading to a pressure rise caused by friction, after which the well is left static. On the other hand, there is much less cutting buildup since the cutting rate is far smaller. At a ROP of 30 m/h, the simulation results are shown in Figure 22. After 1500 s, you may see cuttings at the well’s bottom. At this mass rate, which corresponds to a Rate of Penetration (ROP) of 30 m/h the concentration is more noticeable. The concentration of the cuttings is 0.01 in Figure 23, which is almost three times higher than the concentration of 0.00355 in Figure 20. As the ROP grows, the cutting concentration visibly increases. Higher ROP produces a larger volume of cuttings I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 15 of 21 in a shorter time, increasing the solid phase concentration in the annular flow. This results in sharper gradients and stronger discontinuities in volume fraction, velocity, and pressure fields. Figure 24 shows that our ultimate steady-state BHP is 201.9 bar, which is higher than the 199.8 bar result from the prior simulation. Because of the increased cutting concentration, the BHP has increased. Also, between 2000 s and 4000 s, there is a change in BHP of 2.1 bar (201.9 - 198.8), which is the hydrostatic contribution of the cuttings. Additional research is required to confirm that pressure pulses are the source of the peak at 2000 s. There is a rising tendency in the maximum BHP value, which means that BHP grows as ROP does. This test case demonstrates that the WENO-limited RKDG scheme remains stable even under high rates of penetration (ROP), which typically introduce steep changes in bottom hole pressure (BHP) and cuttings concentration. 0 0.5 1 1.5 2 2.5 3 3.5 4 Cutting Concentration 10-3 -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 19: Case 4: t = 1500 s, mass rate = 0.166 kg/s, plotting of well depth vs cuttings con- centration. 2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 Cutting Concentration 10-3 -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 20: Case 4: t = 4000 s, mass rate = 0.166 kg/s, plotting of well depth vs cuttings con- centration. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 193 194 195 196 197 198 199 200 201 B ot to m ho le P re ss ur e [b ar ] Figure 21: Case 4: t = 4000 s, mass rate = 0.166 kg/s, plotting of bottom hole pressure vs time. 0 0.002 0.004 0.006 0.008 0.01 0.012 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 22: Case 4: t = 1500 s, mass rate = 0.498 kg/s, plotting of well depth vs cuttings con- centration. Case 5: In the present scenario, we discretize the well into 25 boxes and 50 boxes to examine the efficiency of proposed numerical scheme and effect of discretization on the I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 16 of 21 0.0085 0.009 0.0095 0.01 0.0105 0.011 0.0115 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] Figure 23: Case 4: t = 4000 s, mass rate = 0.498 kg/s, plotting of well depth vs cuttings con- centration. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 193 194 195 196 197 198 199 200 201 202 B ot to m ho le P re ss ur e [b ar ] Figure 24: Case 4: t = 4000 s, mass rate = 0.498 kg/s, plotting of bottom hole pressure vs time. two phase flows. Figure 25 displays the impact of increasing the number of boxes on the analysis. The transition zone between pure liquid and liquid/cuttings becomes sharper as the number of boxes increases. Clearly, steady state cutting concentration value becomes more precise with a finer grid (0.07 for 25 boxes vs 0.077 for 50 boxes). Thus, refining the grid brings us closer to the true value, which is 0.067. Figure 26 demonstrates that using 50 boxes results in less smearing of the cutting concentrations at 1500 s compared to 25 boxes. Similarly, in Figure 27, for 50 boxes and 2000 s, we observe a less smeared front. The steady state concentration values at the end differ between the two cases. Figure 28 shows steady state conditions at 3000 s, while in Figure 29, although we reach steady state, the cutting concentration values differ at 4000 s. Figure 30 illustrates the outlet rate, which is observed to be 24 kg/s initially. This is due to the pumping of liquids and cuttings from the bottom of the well. Nevertheless, the liquid flow rate diminishes to steady-state value of 22 kg/s, once the cuttings arrive at the outlet. Figure 31 displays the impact of numerical diffusion on the liquid mass rate. Initially, we observe a transit time of 2172 s for the cutting in the well. The simulation starts at 150 s. The transition is quicker and closer to the theoretical value when using 50 boxes compared to 25 boxes. A smooth curve showing a linear rise in the cutting mass rate is shown in Figure 32 from 2000 s to 2700 s. Four thousand seconds is the time required to achieve a steady state. Here, there is no difference between the cutting mass rate at the well’s top and bottom. The impact of numerical diffusion becomes apparent once again when the cuttings reach the well’s exit, when there is a rapid rise from 0 to 4 kg/s. Figure 33 demonstrates the numerical diffusion effect. It is observed that by using 50 boxes the steady state attained more rapidly as compared to 25 boxes. It takes 3000 s to achieve the cutting mass rate associated with 50 boxes and 4 kg/s. This result is better when compared with the result of 25 boxes, takes 1000 s more. This result shows that the curve display a complete transition. Figure 34 shows the effect of numerical diffusion on BHP. We can observe an increase in BHP from 201− 205 bar at 200 s, which is due to friction. As the cuttings start to distribute in the well, there is a gradual build-up of pressure. The maximum BHP is 222 bar for the solid curve and 224 bars for the dashed curve. When I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 17 of 21 the number of boxes is 50, we observe a difference of two bars in the maximum BHP value. This is because we obtain a more accurate cutting concentration value with 50 boxes. The proposed numerical scheme demonstrates improved efficiency and accuracy as the number of grid points increases, indicating good convergence behavior with mesh refinement. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] 25 boxes 50 boxes Figure 25: Case 5: t = 500 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concentration. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] 25 boxes 50 boxes Figure 26: Case 5: t = 1500 s ,mass rate = 4 kg/s, plotting of well depth vs cuttings concentration. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] 25 boxes 50 boxes Figure 27: Case 5: t = 2000 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concen- tration. 0.055 0.06 0.065 0.07 0.075 0.08 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] 25 boxes 50 boxes Figure 28: Case 5: t = 3000 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concen- tration. 5. Conclusions This study presented a high-order TVD Runge-Kutta Discontinuous Galerkin (RKDG) numerical scheme with WENO limiters to analyze the behavior of two-phase flow in drilling operations using the drift-flux model. The proposed method effectively handled the non- linear nature of transport processes, accurately captured sharp interfaces, and suppressed non-physical oscillations near strong gradients. Through a series of test cases, a clear relationship was observed between increasing rate of penetration (ROP), cuttings concen- tration, and bottom hole pressure (BHP). These results demonstrate that the developed I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 18 of 21 0.066 0.068 0.07 0.072 0.074 0.076 0.078 0.08 0.082 Cutting Concentration -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 W el ld ep th x [m ] 25 boxes 50 boxes Figure 29: Case 5: t = 4000 s, mass rate = 4 kg/s, plotting of well depth vs cuttings concen- tration. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] -5 0 5 10 15 20 25 M as sr at e of L iq ui d at o ut le t [ kg /s ] Figure 30: Case 5: t = 4000 s, mass rate = 4 kg/s, plotting of liquid mass rate out vs time. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] -5 0 5 10 15 20 25 M as sr at e of L iq ui d at o ut le t [ kg /s ] 25 boxes 50 boxes Figure 31: Case 5: t = 4000 s, mass rate = 4 kg/s, plotting of liquid mass rate out vs time. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 0 0.5 1 1.5 2 2.5 3 3.5 4 C ut tin g M as sr at e ou t [ kg /s ] Figure 32: Case 5: t = 4000 s, mass rate = 4 kg/s, plotting of cuttings mass rate out vs time. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 0 0.5 1 1.5 2 2.5 3 3.5 4 C ut tin g M as sr at e ou t [ kg /s ] 25 boxes 50 boxes Figure 33: Case 5: t = 4000 s, mass rate = 4 kg/s, plotting of cuttings mass rate out vs time. 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time [s] 190 195 200 205 210 215 220 225 B ot to m ho le P re ss ur e [b ar ] 25 boxes 50 boxes Figure 34: Case 5: t = 4000 s, mass rate = 4 kg/s, plotting of bottom hole pressure vs time. scheme can provide reliable and accurate predictions for transient drilling hydraulics, of- fering an efficient tool for optimizing operational parameters. However, the current model I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 19 of 21 does not account for thermal effects or complex particle interactions like settling, agglom- eration, or erosion. Future work can extend this study by incorporating more detailed two- fluid models, considering variable fluid properties, and including temperature-dependent behavior. References [1] P. H. Tomren, A. W. Lyoho, and J. J. Azar. Experimental study of cutting transport in directional wells. SPE Drilling Engineering, 1(01):43–56, 1986. [2] R. A. Sanchez, J. J. Azar, A. A. Bassal, and A. L. Martins. The effect of drillpipe rotation on hole cleaning during directional well drilling. In SPE/IADC Drilling Conference. OnePetro, 1997. [3] M. M. Huque, S. Butt, S. Zendehboudi, and S. Imtiaz. Systematic sensitivity anal- ysis of cuttings transport in drilling operation using computational fluid dynamics approach. Journal of Natural Gas Science and Engineering, 81:103386, 2020. [4] M. M. Huque, M. A. Rahman, S. Zendehboudi, S. Butt, and S. Imtiaz. Experimental and numerical study of cuttings transport in inclined drilling operations. Journal of Petroleum Science and Engineering, 208:109394, 2022. [5] A. Zakerian, S. Sarafraz, A. Tabzar, N. Hemmati, and S. R. Shadizadeh. Numerical modeling and simulation of drilling cutting transport in horizontal wells. Journal of Petroleum Exploration and Production Technology, 8:455–474, 2018. [6] M. Pelanti and R. J. LeVeque. High-resolution finite volume methods for dusty gas jets and plumes. SIAM Journal on Scientific Computing, 28(4):1335–1360, 2006. [7] F. Dobran, A. Neri, and G. Macedonio. Numerical simulation of collapsing volcanic columns. Journal of Geophysical Research: Solid Earth, 98(B3):4231–4259, 1993. [8] J. Dufek and G. W. Bergantz. Dynamics and deposits generated by the kos plateau tuff eruption: Controls of basal particle loss on pyroclastic flow transport. Geochem- istry, Geophysics, Geosystems, 8(12), 2007. [9] N. Zuber and J. A. Findlay. Average volumetric concentration in two-phase flow systems. Journal of Heat Transfer, 1965. [10] J. E. Romate. An approximate riemann solver for a two-phase flow model with numerically given slip relation. Computers & Fluids, 27(4):455–477, 1998. [11] K. K. Fjelde and K. H. Karlsen. High-resolution hybrid primitive–conservative upwind schemes for the drift flux model. Computers & Fluids, 31(3):335–367, 2002. [12] E. Steinar and K. K. Fjelde. Hybrid flux-splitting schemes for a two-phase flow model. Journal of Computational Physics, 175:674–701, 2002. [13] M. Baudin, C. Berthon, F. Coquel, R. Masson, and Q. H. Tran. A relaxation method for two-phase flow models with hydrodynamic closure law. Numerische Mathematik, 99:411–440, 2005. [14] S. T. Munkejord, S. Evje, and T. Fl̊atten. The multi-stage centred-scheme approach applied to a drift-flux two-phase flow model. International Journal for Numerical Methods in Fluids, 52(6):679–705, 2006. I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 20 of 21 [15] J. L. Modisette and J. P. Modisette. Transient and succession-of-steady-states pipeline flow models. In PSIG Annual Meeting. OnePetro, 2001. [16] W. H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. [17] J. S. Hesthaven and T. Warburton. Nodal Discontinuous Galerkin Methods: Algo- rithms, Analysis, and Applications, volume 54 of Springer Texts in Applied Mathe- matics. Springer, 2008. [18] J. X. Qiu and C.-W. Shu. A numerical study for the performance of the runge- kutta discontinuous galerkin method based on different numerical fluxes. Journal of Computational Physics, 212:540–565, 2006. [19] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of discon- tinuous galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis, 39:1749–1779, 2002. [20] B. Cockburn and C.-W. Shu. Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws ii: General framework. Mathematics of Computation, 52:411–435, 1989. [21] B. Cockburn and S.-Y. Lin. Tvb runge–kutta local projection discontinuous galerkin finite element method for conservation laws iii: One-dimensional systems. Journal of Computational Physics, 84:90–113, 1989. [22] B. Cockburn and C.-W. Shu. Runge–kutta discontinuous galerkin methods for convection-dominated problems. Journal of Scientific Computing, 1:173–261, 2001. [23] B. Cockburn and C.-W. Shu. The runge–kutta local projection p1-discontinuous galerkin method for scalar conservation laws. Mathematical Modelling and Numerical Analysis, 25:337–361, 1991. [24] C.-W. Shu. Total-variation-diminishing time discretizations. SIAM Journal on Sci- entific and Statistical Computing, 9:1073–1084, 1988. [25] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock- capturing schemes. Journal of Computational Physics, 77:439–471, 1988. [26] R. Y. Adrian et al. Numerical simulation on the behaviour of celular stainless-steel beam (cssb) at elevated temperature. Journal of Engineering Sciences and Technology in Physics, 19(1):25–39, 2024. [27] B. AlHadeethi et al. Numerical modeling of flow pattern with different spillway locations. Mathematical Modelling of Engineering Problems, 11(5):1219–1226, 2024. [28] M. A. Iqbal. Use of ausmv scheme for simulation of particle transport. M.sc. thesis, University of Stavanger, Norway, 2015. [29] C.-W. Shu. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. Technical Report NASA/CR-97-206253, ICASE Report No. 97-65, NASA ICASE, 1997. [30] J. X. Qiu and C.-W. Shu. Runge–kutta discontinuous galerkin method using weno limiters. Journal of Scientific Computing, 26:907–929, 2005. [31] X. Zhang and C.-W. Shu. On maximum-principle-satisfying high order schemes for scalar conservation laws. Journal of Computational Physics, 229:3091–3120, 2010. [32] X. Zhang and C.-W. Shu. Positivity preserving high order discontinuous galerkin I. Khan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6181 21 of 21 schemes for compressible euler equations with source terms. Journal of Computational Physics, 230:1238–1248, 2011. [33] E. A. Az-Zo’bi. On the reduced differential transform method and its applica- tion to the generalized burgers-huxley equation. Applied Mathematical Sciences, 8(177):8823–8831, 2014. [34] E. A. Az-Zo’bi. A reliable analytic study for higher-dimensional telegraph equation. Journal of Mathematics and Computer Science, 18(4):423–429, 2018. [35] E. A. Az-Zo’bi et al. Numeric-analytic solutions of mixed-type systems of balance laws. Applied Mathematics and Computation, 265:133–143, 2015. [36] J. A. Haider et al. Results validation by using finite volume method for the blood flow with magnetohydrodynamics and hybrid nanofluids. Modern Physics Letters B, 38(24):2450208, 2024. [37] K.K. Fjelde and K. H. Karlsen. High-resolution hybrid primitive-conservative upwind schemes for the drift-flux model. Computers & Fluids, 31(3):335–367, 2002. [38] E. Steinar and K. K. Fjelde. Hybrid flux-splitting schemes for a two-phase flow model. Journal of Computational Physics, 175(2):674–701, 2002. [39] C.-W. Shu. Tvb uniformly high-order schemes for conservation laws. Mathematics of Computation, 49(179):105–121, 1987. [40] Gottlieb and C.-W. Shu. Total variation diminishing runge-kutta schemes. Mathe- matics of Computation, 67(221):73–85, 1998.