Mathematical Modelling of Sewage Overflow Through Pipe-Manhole Drainage Sewer Systems Using CFD; A Case of Mbale City, Eastern Uganda

1 Department of Mathematics, Faculty of Science and Education, Busitema University, P. O. Box 236 Tororo, Uganda. 2 Department of Computer Engineering, Faculty of Engineering, Busitema University, P. O. Box 226, Busia, Uganda. 3 Department of Mathematics, Faculty of Science and Education, Busitema University, P. O. Box 236 Tororo, Uganda. * Author for Correspondence Email: ichbinrao@gmail.com or rawichi@sci.busitema.ac.ug.


INTRODUCTION
In urban areas, drainage sewer networks contain a large number of structures. In urban development history, Municipal Sewage systems were built to rapidly collect rain runoff, wastewater and sewage (Esemu et al., 2020). These network structures consist of pipes, manholes, pumping stations, force mains, and other facilities required to collect and transport wastewater (Esemu et al., 2020;Moeini and Afshar, 2017). Manhole is one of the most common features of an urban drainage sewer system. It is usually placed at a drainage network connection where a change in pipe direction, size or gradient is necessary.
Research on urban drainage pipelines focuses on hydraulics such as pipe slopes and flow rate so that sewage and faecal sludge are to be delivered efficiently (Han and Reddy, 2012). The flow inside a manhole is normally complex, highly turbulent and possibly multiphase (Beg et al., 2017). When the flow enters a manhole, it experiences a sudden expansion and contraction involving energy losses. The drainage system is usually designed to operate without surcharging. In normal condition, the sewage level in the manhole remains below the crown level of the conduit. Thus, the flow is only gravity-driven and considered as an open channel flow. During extreme rainfall times, the stormwater gets infiltrated into the sewer system and the sewage level increases, thus the system becomes surcharged making a significant contribution to the overall head losses in the network (Beg et al., 2017). The hydraulic regime within a sewer system is a function of several design parameters including pipe material, sewer function, diameter and gradient. It also suggests that a minimum gradient should be achieved for a specific diameter. As a result, some older sewers do not meet the modern requirements for self-cleaning that result in high blockages (Esemu et al., 2020;Swamee and Sharma, 2013).
Most cities around the world have sewage systems that combine sanitary and stormwater flows within the same network. This is why these networks are known as combined sewage systems (CSS) (Chen et al., 2013). This discharge to the environment, known as combined sewage overflow (CSO) contains biological and chemical contaminants that create major environmental and public health hazards. Therefore, an urban catchment without an effective sewer system may easily encounter inundation and other consequential problems. Mbale Municipality sewer system was built as a sanitary sewer, however during rainy seasons stormwater enters into the system through nonauthorized connections, by some residents, to the sewer system and open manholes hence it experiences CSO. Although the sewer network in Mbale was established as a separate sewer, there are reports of illegal stormwater connections, especially in the Namatala catchment. There are several residential properties that connect their stormwater drains to the sewerage network, which causes siltation in the pipelines and treatment ponds (MWE, 2018;Muwanga, 2015). The sewer system has only two operational convectional sewerage treatment ponds that receive and treat waste from the domestic and industrial piped networks with many of the houses not connected to the urban sewer network. It was constructed according to Bernoulli equation in the early 1950s for the supply of clean water and the collection, transportation and treatment of human excreta and greywater within the Municipality (MWE, 2018;Muwanga, 2015). This means that the flow conditions are based on the Bernoulli principle. The principle implies that an increase in the speed of a fluid occurs simultaneously with a decrease in static pressure or a decrease in the fluid's potential energy (Batchelor and Batchelor, 2000).
The sewer network has not been developed and expanded to meet the capacity of the rapidly growing population and industries in Mbale urban area and its surroundings. It is not uncommon to find overflowing sewage spilling on the streets and these pose health hazards to the residents (MWE, 2018;Muwanga, 2015). This is attributed to drainage system designs, the geographical location of the town, rapid urbanisation and industrialisation. Therefore, there is need for optimal sewer design through improving pipe diameters to manholes, slopes and pressure drop that will consequently minimize capital investment whilst ensuring a good system performance under specific design criteria.
The open-source CFD software OpenFOAM is equipped with a well-designed C++ library that allows the numerical simulation of various engineering applications (Greenshields, 2019). Through its object-oriented structure and the opensource code concept, OpenFOAM is flexible and can be adjusted to very specific environmental fluid flow problems.
This study presents a Mathematical model for the sewage overflow through pipe-manhole drainage sewer systems using the open source CFD software OpenFOAM based on interFoam solver.

LITERATURE REVIEW
The highlights of the key achievements of the National Water and Sewerage Corporation (NWSC) on sewer mains extensions per year where 23.7 km in Financial year (FY) 2015/16 and 26.14 km in FY 2016/17, 24.02 km in FY 2017/18 (Muwanga, 2015). With these figures of extension countrywide, many of the households in urban centres like Mbale city are yet to be connected to public sewer networks. Sewer networks are expected to improve sanitation services in urban areas of Uganda which is still below the National Plan. This indicates that the NWSC has work in expanding and connecting houses to the sewer system. Thus, there is need for increased household connections to the public main sewers to achieve optimized usage of sewerage systems.
An important consideration for sewer performance is the increased volume flow in sewers resulting from increasing sewage efficient sewer connection practices. Large amount of sewage is reliant on the flow of sewers as a means for its movement through the system. Thus, a reduction in flow may reduce the ability of a sewer to efficiently and effectively remove large amounts of sewage (Moeini and Afshar, 2017). Therefore, if the benefits of improving sewer networks are to be fully realised, consideration must be given to how sewer design can be updated to accommodate increased volumes of sewer flows (Esemu et al., 2020;Petit-Boix et al., 2016).
To achieve this, fluid flow equations are considered as below: Conservation of mass, Conservation of momentum where V, ρ, P, µ and F denotes fluid flow velocity, density of fluid, fluid pressure, the dynamic viscosity of fluid and body force respectively.

Sewer Networks Modelling
The solution of sewer network optimization problems requires determination of pipe diameters, flow velocity and pipe slopes (Esemu et al., 2020;Moeini and Afshar, 2017). Finding an optimal least-cost design of sewer networks is difficult because of numerous complex hydraulic and engineering constraints that lead to nonlinear equations (Swamee and Sharma, 2013). Moreover, the size of sewer networks is generally large, therefore, leading to large scale optimization problem which is generally difficult to solve (Moeini and Afshar, 2017). Sanitary sewer network, as one of the physical networks is an essential structure of any urban area. Urban areas without an effective sanitary sewer networks may encounter many problems, such as public health threats and environmental damages. Having a suitable and cost-effective sewer network is normally interpreted as finding solution for sewer network design optimization problem that minimizes infrastructural cost, without violating operational requirements. The cost of construction of a sewer system can be significantly reduced if the system layout, pipe diameters and pipe slopes can be effectively optimized (Afshar et al., 2011).
Computational Fluid Dynamics is applied to gain insight into most fluid dynamic processes and flow phenomena like in sewer lines to add value in the analysis of urban drainage sewer systems. According to Park et al. (2016), CFD presents itself as a useful tool for investigating domain space for physical system design, performance variables and for diagnosing or troubleshooting system behaviour. Beg et al. (2017) and Beg et al. (2016) conducted a CFD study on the comparison of flow hydraulics in different manhole types and numerical investigation of the flow field inside a manhole-pipe drainage system respectively.

MATERIALS AND METHODS
The study depended on combined CFD modelling and numerical methods. The numerical modelling part utilised open source CFD tools OpenFOAM and the solver interFoam was used to model the flow phenomena. In CFD modelling, Navier-Stokes equations of a one-dimensional incompressible fluid flow was used as shown in Equation (1).
A transport equation for the volume fraction of each cell in OpenFOAM solver interFoam was solved to find the shape of the free surface (air/liquid interface) given in form of conservation of linear momentum as: where α denotes the volume fraction.
According to Haghighi and Bakhshipour (2012) and Akgiray (2004), steady-state flow in sewers is described by Manning's equation to find sewer hydraulics as; where V = flow velocity, r = hydraulic radius of pipe = 4 , S = slope (sewer gradient), n = Manning's coefficient and D = diameter of pipe(sewer).
Thus, sewage flow rate given as; where A is the wetted cross-sectional area of the pipe given as, = 4

Study Area
The study area was Mbale Municipality, the central business area of Mbale city located in eastern Uganda. It was selected as the study site because of the existence of the public sewage network since the colonial times of 1950s and currently managed by National Water and Sewerage Corporation (NWSC). The area has currently a population of at least 100,000 inhabitants and the sewer system was originally constructed for a population of around 45,000 people (Muwanga, 2015). Most of the manholes which are interlinking the sewer system are overburdened with the overloaded sewage flow which has led to the sewage outflow through them.
The Municipality has two isolated sewer system networks: One system discharge into the Namatala treatment plant and the other one discharges into the Doko treatment plant. The Namatala ponds have a capacity to treat about an average flow of 1,030 m 3 per day and the Doko ponds on average have a capacity to treat about a daily flow of 1,120 m 3 .

Research Design
The study considered the following parameters flow velocity, pressure drop and pipe diameter to understand the flow through a pipe which affects the flow dynamics in a sewer system. Data for model variable were collected from documentation and some relatively small samples of sewer diameters were also collected from the NWSC offices. The CFD computations involved carrying out structured and unstructured computational mesh discretisation for Navier-Stokes equations and the VOF using cross-sectional design of sewers. The mathematical model algorithm design was done using CFD and simulations were carried out using the open source CFD software OpenFOAM.
According to Haghighi and Bakhshipour (2012), for a given layout, a feasible sewer design is defined as a set of pipe diameters, slopes and pump location indicators which satisfies all the constraints.
Referring to Afshar et al. (2016), typical constraints of sewer networks design are: i) Keeping flow velocity between minimum and maximum bounds for self-cleaning capability and preventing scouring respectively, ii) For each manhole, placing upstream pipes at appropriate higher elevations than the downstream ones and assigning the outlet pipe's diameter equal to or greater than the upstream inlet pipes, iii) Maintaining the proportional water/wastewater depth under the specified maximum value, iv) Choosing sewer diameter that is effective for delivery from the available commercial lists, maintaining the minimum buried depth to prevent damages from the traffic loads and other surface flooding like erosions.

Sewer Networks Design Optimization
A cost function for optimal sewer network design that aims to find pipe diameter and slopes which minimizes pressure drops and sewer system cost was derived based on a model proposed by Esemu et al. (2020) and Zaheri et al. (2019). Excavation and pipe installations costs are the main part of construction costs which are functions of the pipe slopes and pipe diameters (Moeini and Afshar, 2017). Sewer network design intends to minimize the construction cost of the network. The problem of sewer network design, in its general form, may be formulated as (Afshar et al., 2011); Here CT is the total construction cost of the network, N is the number of pipes in the network, Di is the diameter of the i th pipe, Xi is the average excavation depth of the i th pipe; a function of the pipe slope, hm is the manhole depth, hp and hD are the pump and drop heights respectively.
In the absence of pumps and drops, for every feasible design, the construction cost is estimated from the following equation (Afshar and Rohani, 2012).
where CT = cost function of sanitary sewer network; N = total number of sewer pipes; M = total number of manholes, Li is the length of the i th pipe, Kp is the unit cost of the i th pipe defined as a function of its diameter (Di) and average excavation depth (Xi), and Km is the cost of manhole construction as a function of manhole depth (hm).

CFD Formulation of the Sewer Network Optimal Design
In this formulation, each node of the sewer network was regarded as a computational cell, the nodal elevation of the network; decision variable of the optimization problem was considered as a cell state. The neighbourhood of each cell was defined as the pipes which are connected to each CFD cell. Instead of an ad-hoc transition rule based on the engineering judgment of the physical behaviour of the problem (Afshar et al., 2011), which was used in most CA applications to water resource problems, an alternative mathematically CFD approach is proposed here for higher efficiency.
The proposed updating rule for an arbitrary cell, j; the network node, is derived by requiring that the problem objective function of Equation (5) is minimized, with respect to the cell state, hj , while all other nodal elevations are kept constant. This leads to the following local optimization problem, defined on the neighbourhood of cell j, as: Find hj to minimize: = 1 ( 1 , 1 ) + 2 (( 2 , 2 ) + (ℎ ) ; subject to the following constraints: where subscripts, 1 and 2, refer to the pipes in the neighbourhood of cell j. It is assumed here, for simplicity of notation, that only two pipes are connected to a typical node, j, of the network. For nodes which are at the intersection of more pipes, Equation (6) should be amended to account for the cost of all pipes connected at node j.
Differentiating the objective function, with respect to the decision variable, hj , leads to the following nonlinear equation, defined as: where the first and second terms on the right-hand side consider the derivative of pipe cost, with respect to the j th decision variable (hj), since this cost is a function of diameter (D) and the average pipe cover (X). The third term on the right-hand side also considers the derivative of the manhole cost, with respect to the j th decision variable (hj).

This cost can be easily calculated if Kh(hm)j is known.
Applying Newton-Raphson method for linearization to the Equation (7) results in the updated nodal elevation as: where ∆ℎ = ℎ +1 − ℎ and ℎ +1 is the updated value of the nodal elevation, k is the nonlinear iteration index, and∆ℎ is the change in cell state value. This procedure is repeated for all cells of the network until convergence is reached. The nodal elevation equation is as; With calculated nodal elevations, pipe slopes and diameters can be calculated as defined, and thus the problem is solved.

Mathematical Model Equations
The numerical modelling part utilizes open source CFD tools OpenFOAM v2.3.0. The solver interFoam was used to model the flow phenomena. It considers the fluid system as isothermal, incompressible and immiscible two-phased flow. It uses a single set of Navier-Stokes equations for the both fluids (water and air) and additional equations to describe the free-surface (Lopes et al., 2018). The velocity at free-surface is shared by both phases. The solver deals with Reynolds averaged conservation of mass and momentum for incompressible flow (Jasak, 1996;Rusche, 2002).

Conservation of mass;
∇. = 0 Conservation of momentum; where V is the velocity vector towards Cartesian coordinate system, ρ is the density of the fluid, τ is the shear stress tensor. Here ρg is the acceleration of the fluid because of the gravity and F is the volume surface tension force acting at the interface.
The modified pressure p * (as P−ρgh in OpenFOAM code) is adopted in interFoam by removing the hydrostatic pressure (ρgh) from the pressure P, where its gradient is given as: The stress term is defined as: where, µ is the dynamic viscosity of the fluid.
The volumetric surface tension F is calculated by the continuum surface Force Model (Brackbill et al., 1992) as; where κ is the direction of the surface curvature, σ is the surface tension.
Substituting for ∇P, ∇τ and F in Equation (13) yields; InterFoam solver utilizes Volume of Fluid (VOF) method (Hirt and Nichols, 1981) to capture the free surface location. VOF utilizes an additional volume fraction indicator term α to determine the amount of each fluid amount in a computational cell volume.
= [ 1; 0 0; VOF method also utilizes an interfacial compressive term for the interface region to keep it confined at a certain region (Beg et al., 2017).
where, Vc is the compressive velocity that acts at a perpendicular direction to the interface. This is calculated as; where, Cα is a Boolean term (value is 0 or 1) which activates (Cα = 1) or deactivates (Cα = 0) the interface compressive term.

Computational Mesh
The computational meshes were generated using the mesh generator supplied with OpenFOAM, called blockMesh, generated meshes from a description specified in an input dictionary, blockMeshDict located in the constant/polyMesh directory. The simulations are simplified to two dimensional (2D) cases. In both types of simulations, the geometry is 0.5 m high for sewer pipe/closed channels and 1.0 m high for the open channel and 20 m in length. When running blockMesh for the sewer pipe (closed channel), the created mesh from the blockMeshDict was as shown in Figure (1).

Assessment of Sewer Network
Mbale Municipality existing sewerage system was built as a sanitary sewer network based on pressure head and gravity concept with only two pump stations. All the domestic and industrial sewage is transported through the isolated sewer networks to the two separate treatment plants where it undergoes natural treatment prior to being discharged to the environment. The main sewer line was constructed in the early 1950s to carry about 2000 m 3 per day of sewage for at most 45,000 people. However, with the current population increase in the city by at least 45% since the 1990s (MWE, 2018) and the demand for usage of the public sewer line, the sewage production is about double the amount it is meant to carry which has led to the sewer lines being overburdened hence the overflow through the manholes. This means, more than 800 m 3 per day of sewage is uncollected and therefore, discharged to the environment untreated.
Some of the sewer lines during rainfall seasons are infiltrated by the stormwater either through unauthorized connections of stormwater channels by the public to the sewer lines or through damaged manholes, this has always led to the increased flow in the system and blockages which causes backflow of sewage hence overflow through any opening, more so through the manholes.

CFD Model Simulations Results and Flow pattern visualisation
The standard k−ε and VOF model visualisations generated in the analysis for optimal sewer flow rate; stratified velocity, pressure drop in a sewer pipe/ channel are presented using numerical simulation profiles and plots in terms of volume fraction (α) of the water-air phase as shown in Figure (2) to Figure (8). In this case, the red colour represents water (α = 1), the blue colour represents air (α = 0) and between the two phases is the interface.

Improvement of the Sewer System Infrastructure
The most effective way to optimize a sewer network system based on gravity flow and fluid pressure is by adjusting sewer diameters, flow velocities (minimum and maximum) and slope gradients in order to achieve optimal flow rate to transport and deliver sewage to the treatment plant with minimal overflows through sewers. It is also one of the easiest ways to control flow rate. However, specifying pipe sizes, minimum velocity accurately may help to control those factors that can change over time. For example, in old metallic and reinforced concrete piping systems, as it is observed in the NWSC records, the piping diameters should have been wider in design work than was initially necessary in order to account for material friction and some corrosion and/or scaling. Minimum (selfcleansing) velocity and Maximum velocity for each sewer, where the solid particles will remain in suspension without settling at the bottom of the sewer has to be considered. Self-cleansing velocity should be maintained at least once each day during the working life of the pipe to ensure its efficiency. Maximum velocity must be a velocity that causes no scouring action or abrasion. It depends upon the material used for the construction of sewer.

Figure 8: A graph showing Velocity development of a fluid against time in an open channel
Currently, there are no continuous specific records on the sewer system flow from NWSC, which make the monitoring of the flow system difficult. Thus, the existing sewer system should be monitored using the application of the computation fluid techniques, like the use of OpenFOAM software to always have records of the sewer flow system parameters which tend to affect flow or facilitate flow to help in the identification of advanced overflows in the system. This calls for research departments on computation fluid dynamics at each centre of the NWSC offices to be established.
Currently, the sewerage coverage in the Municipality is poor with only 14% of the households with water supply having a sewerage connection. Therefore, this should be gradually expanded from the current 14% to a level that matches a decent living standard for the city by increasing the number of connections to cater for the increasing population within the city which cause increased sewage output to drastically increase the amount of collection and transportation of sewage to treatment plants. However, this also calls for improved treatment techniques to cater for the increased amount of sewage to these treatment plants. This will reduce operational costs of the Municipal and NWSC authorities thus leading to improved hygiene and health in the town.

Recommendations for Improvement of the Existing Sewer Network
To optimize by minimizing the sewage overflow through pipe-manhole sewer network system for effective fluid flow, the recommended sewer diameters and their corresponding slopes should be adhered to during network design and construction in order to achieve optimal flow rate to transport and deliver sewage to the treatment plants.
Minimum sanitary sewer slopes are established to provide a minimum cleansing velocity of 0.67 ms −1 at full flow or half-full flow for a Manning's pipe roughness (n) of 0.013. Hence to avoid adverse slopes caused by inaccurate construction or settlement, a minimum pipe slope should be considered for the sewer pipes. The recommended minimum slopes (S) for different pipe diameters (D) for sewer systems are shown in Table 1. The recommended flow velocity in sewers should be in the range [0.67ms −1 ,5.5ms −1 ] where 0.67 ms −1 is the minimum (self-cleansing) velocity and 5.5 ms −1 is the maximum velocity for each sewer (i.e. velocity at which, the solid particles will remain in suspension without settling at the bottom of the sewer). Self-cleansing velocity should be maintained at least once each day during the working life of the pipe to ensure its efficiency. Maximum Velocity must be a velocity that causes no scouring action or abrasion, which depends upon the material used for the construction of sewer pipe.