Abstract
Most of the real contaminant problems are defined domains that are geometrically complex and can have different boundary conditions in different areas. Therefore, it is usually difficult to find a solution analytically, so we use the approximate method to generate an approximate function. One answer to this problem is the finite element approach (FEM). This study presents a partial differential equation (PDE) simulation system that uses numerical techniques for the distribution of pollutant concentrations in groundwater in space and time. The movement of the liquid is described by the incompressible steady-state Navier-Strokes equation, while the transport of pollutants is described by the diffusion-convention equation. The variation formulation that forms the basis of FEM and MATLAB is discussed along with the selection of the abstract approximation space and the welfare of the weak formulation. The motivation for this study comes from a specific and considered water body with the discharge of factory effluents on the ground that ends up reducing the quality of groundwater. First, the fluid flow equation is solved to obtain velocity and pressure profiles. Steady-state concentration profiles were obtained for various values of diffusion coefficient (), baseline, and input concentrations. The results showed that decreasing the diffusion coefficient increased the number of pollutants for convective transport and decreased the number of pollutants that diffused from the entrance. Although groundwater is not completely safe, it is concluded that experimental studies are necessary decision-making basis for water resource protection, especially in water pollution emergencies.
1. Introduction
In this paper, a model based on a stationary incompressible 2-dimensional Navier-Stokes equation for fluid flow and a 2-dimensional convection-diffusion equation for polluting water transport (one type of pollutant) is then formulated and the numerical model solved for different boundary conditions and, finally, simulations to predict the movement and concentration of the pollutant in regions occupied via of porous media at different positions and times. The driving motivation for this work comes from a specific and considered industrial problem in the transport and dispersion processes of polluted particles that end up reducing the quality of groundwater [1]. In the literature, water is a basic need of life, used domestically in industry, recreation, and irrigation purposes. Groundwater is found in underground aquifers and provides approximately 97% of the world’s consumable water [2]. In Uganda, many communities, especially the rural poor, depend on streams and swamps for supply of water [3]. In recent years, the emphasis for policymakers and researchers has shifted from water accessibility to water quality problems. Environmental pollution problems require prompt action to prevent the reduction of water quality. The issue of water quality is an important factor in sustainable human development [4]. Groundwater is polluted due to anthropogenic activities on the land, which include the use of agrochemicals and poor disposal of domestic and industrial waste. Increasing trends in the use of fertilizers and pesticides in agricultural production systems and disposal of industrial waste are the reasons for contamination of groundwater [4]. A study on risk factors contributing to microbiological contamination of near subsurface water in Kampala indicated that water from shallow protected springs was polluted [5]. Pollution level increased in the rainy season due to storm water runoff and infiltrated into the groundwater [6]. The impact of pollution on shallow groundwater in Bwaise III Parish, Kampala [7], showed contamination of groundwater; in addition, the pollution levels of water resources in Pece Wetland, Gulu Municipality, indicated that domestic water sources were contaminated [8]. Modelling fluid flow problems falls in the branch of applied mathematics called computational fluid dynamics (CFD), which uses tools from numerical analysis, fluid dynamics, and computer science. Discretization techniques which include boundary element, finite difference, finite volume, and finite element methods have been globally utilized in solving solute transport and distribution of fluids [9]. Problems solved by these discretization techniques include multiphase and single-phase fluid flows in porous medium [10–12], flow in a thermal labyrinth [13], pollute transport in the atmosphere [14, 15], pollute transport in rivers [16, 17], and pollute the groundwater [2]. FEM is a numerical analysis technique which employs the concept of piecewise “approximation” to approximate partial solutions of differential equations. It is a powerful discretization tool which can accurately discretize a domain of any size or shape. The Navier-Stokes equations are fundamental models in fluid dynamics [18] that describe motions of fluids in . These equations solve velocity vector and pressure for any position and time [19]. Many groundwater problems have been modelled by Darcy’s law [1, 12, 20–22], among others. Henry Darcy discovered that water flow through a porous medium was governed by equation (1).
where is the total discharge through the medium, is the permeability of the medium, is the cross-sectional area of the flow, is the pressure gradient, is the viscosity, and is the length over which pressure drop occurs [21].
In this study, fluid flow is modelled by the Navier-Stokes equations instead of the usual Darcy’s law because, under certain assumptions, Darcy’s law is derived from the Navier-Stokes equations, [23]. The momentum balance (equation (2)) can be written as
where by: is the material time derivative of a given particle, τ is the stress tensor, is the external body force per unit mass, and is the drag force/unit volume exerted on the fluid. The momentum balance (equation (1)) for the flow of water in the saturated zone reduces to the well-known Darcy’s law under certain conditions [23]. Moreover, Darcy’s law is reliable for values of . At very low Reynolds’ number, the circulation generated can be neglected, but as Reynolds’ number increases, it becomes important and a noticeable contribution to the resistance to flow. The fact is that the circulation zone becomes smaller at even higher Reynolds’ number. Pollutants generally dissolve and are carried by infiltrating rain water into unsaturated soil above the water table. Pollutants then enter the saturated zone and migrate in the direction of groundwater flow (Figure 1).
2. Materials and Methods
Human health is harmed by suspended particles in water. Particle transportation and distribution of matter suspended in water is associated with fluid motion and turbulence. CFD is the most appropriate modelling approach for studying particle spatial distribution in a particular domain [9]. These methods enable complex and abstract mathematical models or theories to be easily visualized through the use of simulations. In the formulation of groundwater models, errors can arise from conceptual problems due to excluding relevant or including irrelevant physical process in the model or using an inapplicable model for the groundwater solute transport problem and numerical errors like truncation and round-off errors and using wrong input data [1]. The CFD model is determined by the nature of the physical process to be simulated, the objective of the study, and the available resources [24]. Equation parameters can easily be adjusted to observe how these changes will affect the model [25]. Because the mathematical model is the result of coupling multiple physical fields, which restricts the choice of function spaces when using the finite element method, an adequate choice of function spaces for pressure and velocity is required [26]. Mesh generation depends on the complexity of the geometry of the domain and can either be structured, block-structured, or unstructured. For example, in a one-dimensional case, if the computational domain is divided into equal subintervals, we generate a structured mesh, in the case of a complex domain, that requires a fully unstructured mesh. A good and reliable numerical model should be able to simulate a transport phenomenon accurately and suppress the instability arising out of discretization in the computational domain [27]. Groundwater flow and solute transport are described by second-order partial differential equations that can be parabolic, elliptic, or hyperbolic. This classification is essential because the choice of a numerical method should be best suited to the nature of the PDE. For instance, a diffusion-convection process dominated by the convective term approximates a hyperbolic type of equation, while a diffusion-convection process dominated by the diffuse term approximates a parabolic type of equation. FEM is a universal discretization tool for partial differential equations. Its advantage is that it can easily be defined on a complicated geometry and can improve the quality of the numerical solution by fine-tuning the model grid. The mesh need not be uniform; it can be made finer in areas of the domain where greater accuracy is needed. Its downside, however, is that the programming of FEM is more complicated and needs standard software packages. The finite difference method (FDM) is simpler and easier to program, but it is only applicable to a regular domain [28]. The extract from Konikow and Reilly [28] is an illustration of the finite difference and finite element methods in discretizing a domain of a field. The rectangular grid of the FDM fails to cover the entire domain of the field, while the triangular grid of the FEM can easily be constructed as shown in Figure 2.
(a)
(b)
(c)
When using the FDM technique, the domain is divided into uniform subareas and the PDE’s are then replaced by approximating functions which are written in terms of finite subspaces. The time derivative is also approximated by a finite difference expression obtained using a Taylor series expansion [23]. When formulating the finite element method, the variational direct energy balance or weighted residual approaches can be used. The weighted residual and variational approaches are the most commonly used approaches for groundwater models [1]. When using the variational approach under FEM, the appropriate PDE together with its boundary conditions and initial conditions are replaced with a corresponding variational problem. In this approach, the unknown function is approximated throughout the solution domain as in Equation (3)
where are suitable bases or shape functions, coefficients are unknowns, and is the number of nodes.
A finite element Newton method for the solution of steady-state Navier-Stokes equations for 2-dimensional incompressible flows was discussed [29]. In this work, a weak variational formulation of the problem formulation and an unequal order interpolation for pressure and velocity were adopted. A Newton method was used for the nonlinear system of coupled equations written in an incremental form and the Jacobian linear system was solved using a direct algorithm. They explained that the Newton iterative method was a suitable technique because of its efficiency since only a few iterations were sufficient for convergence to a very accurate steady solution provided the initial guess was not chosen too far from the solution. They also explained the use of the incremental form of the Newton iteration, which fully exploited the quadratic nature of the nonlinearity of the Navier-Stokes system and thereafter obtained an algorithm which was optimal both in the imposition of the boundary conditions and with respect to the cancellation of numerical errors. For solutions at high Reynolds’ numbers, the Newton method may fail when the Stokes flow is too far from the nonlinear solution. They then resorted to a continuation stepwise manner by computing different intermediate solutions for smaller values of the Reynolds’ numbers. The Laplace transformation technique has been commonly used to obtain the solution of the diffusion-convection equation because of its simplicity compared to other methods [30]. Other analytical solutions to the convection-dispersion solute transport equation are discussed [31]. A groundwater contamination, diffusion, and spreading model can be solved by Fourier transforms [2]. A similar model could have been done without assuming that all covering parameters are constant because this is highly unlikely in reality. For example, soil conductivity, which measures how fast a fluid can be transported through the soil, is a highly variable quantity. Other models of pollution in river type aquatic systems can be found [4]. To prevent groundwater contamination in both shallow and deep aquifers [32], cut-off walls are often used. Discacciati and Quarteroni [12] utilized a model that couples the Navier-Stokes and Darcy equations and modelled the filtration of incompressible fluids through a porous medium [33]. In this work, the motion of incompressible fluids was described by the Navier-Stokes equations while the filtration process was described by Darcy’s equation. A computational domain divided into two parts was considered, one for the fluid and the other for the porous media, and went ahead to discuss the coupling conditions (Beavers-Joseph-Saffman) including their mathematical justification via a homogenization theory. Interface conditions were also introduced because of the multidomain structure of the model. In this regard, a much better model for pollute transport in a groundwater contaminant problem would require the coupling of the Navier-Stokes equations with a groundwater flow equation in the porous media together with a diffusion-convection equation for pollute transport in the two regions.
2.1. A PDE Model for Groundwater Flow
To achieve the objectives of groundwater flow and pollute transport problem, a groundwater flow equation based on an incompressible steady 2-dimensional Navier-Stokes equation and a 2-dimensional convection-diffusion equation describing pollute transport is constructed. The Navier-Stokes equations are derived on the basis of the conservation principles of mass (continuity) and linear momentum (Newton’s second law of motion). The first law of thermodynamics (energy) is not considered because of the assumption that heat transfer in groundwater is negligible. Finally, the boundary conditions both Dirichlet and Neumann are stated, together with reasonable assumptions. The density of water is assumed constant, that is, it is not affected by variations in concentration. The pollute transport equation is derived based on a Fickian model which assumes that the waste is transported mainly due to the concentration gradient and that the dispersive flux occurs in a direction from higher concentrations towards lower concentrations [34]. In this work, the PDE model for pollute transport in groundwater was formulated by first deriving an equation for groundwater flow using Navier-Stokes equations.
2.2. Derivation of the Navier-Stokes Equation
The flow of water is described using steady incompressible Navier-Stokes equations. These equations are derived on the basis of the conservation principles of mass (continuity) and linear momentum. Let describe the domain covered by the flowing water. A fluid particle (Figure 3) is considered to be moving through at time .
The Lagrangian or material coordinate is . It is the observed trajectory of a specific fluid particle as the fluid flows while is a Eulerian coordinate because the fluid particle is observed through a fixed region of interest [35].
Also, let be a subset of . Define the function such that det for all the change of the particle’s position is described by equation (4).
Define to be a vector field on at time , the mass density, the mass density of the reference frame the force density, the surface force density also known as traction or contact forces, and the outward vector to . By the conservation of mass property, the total mass remains constant as shown in equation (5).
Since the total mass remains constant, the rate of change of mass in is given by
To simplify the term on the left-hand side of equation (6), Reynolds’ transport (Theorem 1) leading to equation (7) was used.
Theorem 1. For an arbitrary single-valued scalar function with continuous derivates and an arbitrary control region with boundary outward-pointing unit-normal , and the local velocity of the fluid across the following integral relation holds: where is an element of length on
Source: Reinstra and Hirschberg [36] and Powers [37].
In simple terms, the rate of change of a quantity, , in is a result of the rate at which the quantity is being produced in plus the net amount that crosses the boundary Reynolds’ transport theorem is used to translate the integral conservation laws into differential conservation laws. This helps in formulation of the basic conservation laws of fluid dynamics. To simplify calculations, the boundary integral on the right-hand side of equation (7) is transformed to the domain integral using Gauss divergence (Theorem 2).
Theorem 2. Lets a 2-dimensional domain with boundary ds an element of length (arc length) on and associated unit outward normal vector n with a vector field on then source: [36]. Gauss’ divergence theorem 1.1 shows the relationship between a domain and surface integral. Applying the Gauss divergence theorem, the right-hand side of equation (7) can be written as Substituting equations (8) into (7) yields equation (9) as Applying equation (6) to equation (9), Since the integral over is arbitrary, equation (10) can be written as
Equation (11) is known as the continuity equation. Applying the incompressibility condition to equation (11), the continuity equation now becomes . The conservation of momentum equation was derived from Newton’s second law of motion. From , where is the force action on the particle, is the mass, is the acceleration, and is the velocity of the particle. The total force acting on the subdomain is the sum of the field force and the surface force .
By Newton’s second law of motion, where is an element of length on . Considering the action of the body or field forces and the normal component of the surface forces where is the fluid stress tensor, the formula is written as
Considering the component of the vectors , and where denotes the row of the matrix , equation (12) is rewritten as
Using the divergence theorem, the boundary integral on the right-hand side of equation (13) is written as
Moreover, by Theorem 1, the term on the left-hand side of equation (13) can be written as
And by the incompressibility condition, the following is obtained: substituting equations (14) and (16) into equation (13),
Since the integral over is arbitrary, equation (17) becomes
where
Introducing constitutive relations into equation (18), and ; where is the volume viscosity, is coefficient of shear viscosity, and we set the deformation tensor in is the small strain tensor [38]. Substituting for , and for in equation (18),
where is the Laplacian of . Applying the incompressibility condition to equation (20),
Therefore, equation (21) is obtained as Substituting equation (22) into equation (18),
Assuming the mass density to be constant, equation (23) can be written as
where is the kinetic viscosity.
To make an analysis of the physical problem when the equation parameters of equation (24) change, scaling of the Navier-Stokes equation (24) was carried out. Normalization of these scales resulted in the formulation of dimensionless groups such as Reynolds’, Froude’s, Euler’s, Weber’s, Prandtl’s, and Mach’s numbers which represent the relative importance of various parts of the Navier-Stokes equations and are also key in determining the flow regimes. Scaling also helped to reduce the number of variables. For instance, the combined effect of both viscosity and density can be determined by one dimensionless variable called Reynolds’ number.
2.3. Nondimensionalization of Navier-Stokes Equations
Let
where , , and are characteristics speed, length, and time, respectively. From expressions in (25), we obtain
In equation (26), the convective termbecomes
The pressure term becomes
The viscous term in (29) becomes Substituting equations (26)–(30) into equation (24), we obtain
Also let , then
Introducing a nondimensional source term and substituting equations (34) and (35) into equation (33) and further substituting and multiplying , we obtain
where is the Reynolds number. It is a dimensionless ratio of inertial to viscous forces used to determine the flow regime. Dropping the tilde notation and assuming steady state, equations (36) and (28) form
And where equation (36) represents the dimensionless steady incompressible Navier-Stokes equations.
Pollute transport is due to the combined effect of diffusion and advection in a fluid. The PDE describing it as derived from the principle of conservation of mass using Fick’s law [30]. Underground water flows through the pores (space between soil particles) in the soil. Some soils are more compact than others. The ease with which water flows through the soil depends on the soil porosity denoted by such that
2.4. Diffusion-Convection Equation
Consider a porous domain and with denoting the density or concentration of the pollute in mass per unit area, denotes the flux, that is, mass per unit length per unit time crossing the boundary the rate at which the pollutant is increasing or reducing in due to the source term. The total mass of pollute in is and the amount of pollute that crosses the boundary is where is an element of length on and is a unit outer normal vector to . The pollutant crosses the boundary in three ways: (1)Advection: this is due to the bulk flow of water molecules where pollute particles are carried along the streamlines. It is given by (concentration velocity)(2)Molecular diffusion: this is caused by the random motion and collusions of pollutant particles. This occurs even when the fluid mass is at rest(3)Mechanical dispersion: this is due to the velocity variations caused by the different flow paths that the groundwater takes
The total increase or decrease of the pollutant in due to the source terms is given by By the mass balance principle, the rate of change of the total mass of the pollutant in equals the net rate of pollutant mass that flows through the boundary plus the net rate of increase or decrease of pollutant mass in ; this is summarized as
Thus,
By the divergence theorem, the first term on the right-hand side of equation (40) can be written so that equation (41) becomes
Applying the Reynolds transport theorem, Gauss divergence theorem, and the incompressibility condition, the derivative term on the left-hand side of equation (42) is written as
which further simplifies to
Since is arbitrary, equation (44) implies
The contribution to flux from the effect of diffusion (molecular) and solubility (mechanical dispersion) according to Fick’s law:
where is the diffusion coefficient. The equation means that the flux is proportional to the negative gradient of concentration and the pollution will diffuse from a region of high to low concentration. The total flux is the sum of advective, molecular, and mechanical dispersions [34]. Therefore,
where is the convection velocity also known as advection velocity, substituting equation (47) into equation (45) and assuming the diffusivity is constant, and by the incompressibility condition we obtain
which simplifies to
where is the average velocity vector of groundwater obtained from solving equation (37) and is the source term. Combining equations (37) and (49), the complete model equations which are a set of highly nonlinear partial differential is as follows;
With initial condition and boundary conditions prescribed by
where and represent the domain boundary prescribed by Dirichlet and Neumann conditions, respectively, and .
Also, it is expressed as
where , and are known data to be prescribed later. Dirichletian boundary conditions are called essential because they are imposed explicitly on the solution, that is, on the test space, while Neumann boundary conditions are called natural because they are implicitly defined in the variationally weak formulations in which the model is based on the following assumptions: (a)It is known a priority regions where groundwater flows(b)Groundwater flows at a rate which depends on the permeability of the soil(c)No chemical reactions occur during the interaction of pollutants with water
The existence of solutions to equation (50) is discussed below.
2.5. Existence of Solutions to the Model Equations
The following definitions are used to describe the functional spaces.
Definition 3. The space of real valued functions that are square integral on with respect to the Lebesgue measure is denoted by the relationship The space is equipped with the norm and scalar products in [39]
Definition 4. The Sobolev space is denoted by the expressions in [39] This space is equipped with the scalar product in that induces the norm in [39]
2.6. Functional Spaces
Let be test spaces for velocity, pressure, and pollute concentration, respectively, as described below;
The velocity solution space and concentration solution space are also described as below:
2.7. Weak Formulations
The weak formulations of the model equations are derived for this purpose from Theorem 5.
Theorem 5. (Green’s formula) let be an open bounded regular set of class If and u are functions of , then they satisfy where is the outward unit normal to and [39]
The weak form of the system equation (52) is as follows.
To find such that equation (65) is expressed as
By applying Green’s formula on the viscous and pressure terms of the momentum balance equation, and on the diffusive term of the convection-diffusion equation, the following is obtained:
In defining the bilinear forms as in expressions by and the trilinear form by
Equation (66) is expressed as
The existence and uniqueness of solutions to problems in equations (71) and (73) are briefly discussed.
2.8. Existence and Uniqueness of Solutions with the Weak Formulation
The existence and uniqueness of solutions to the weak formulation (equation (71)) is guaranteed by the Lax-Milgram theorem [39]. Alternatively, from the weak formulation (equation (71)), the divergence of free subspace when introduced can be given by and the weak formulation (equation (71)) with the assumption that can be formulated as follows:
Lemma 6. Let be a solution to a problem in equation (75), the there exists a unique such that is a solution of problem (equation (71)) [40]. The existence of the solution in Eequation (75) can be established from the continuity and coercivity arguments which guarantee well-posedness of the problem [40]. For this purpose, the following lemmas and theorems are stated.
Lemma 7. The trilinear form is continuous.
Proof. The proof follows from the Cauchy-Schwarz inequality, Sobolev embedding theorem, and Holder’s inequality, from which the following can be written as [41].
Lemma 8. The trilinear form has the following properties: (a)(b) [41]
Theorem 9. There exists a weak solution to the Navier-Stokes problem (equation (36)) and a constant such that
Proof. The proof which follows from the Galerkin method involves choosing a finite dimension subspace of , and writing the solution as a linear combination of the functions that form a basis for since the sequence is bounded in (result from the Galerkin equations), applying the results of Lemma 8, Poincare’s inequality, and the Bolzano Weierstrass theorem, there exists a weakly convergent subsequence such that Using compactness and embedding argument, it can be concluded that is a weak solution to Navier-Stokes problem (equation (36)) [42].
The uniqueness of the solution is discussed in Theorem 10 below.
Theorem 10. There exists a constant such that the solution of the Navier-Stokes problem (equation (49)) is unique, if . Proof: [42]. Remarks on Theorem 10 imply that the uniqueness of the solution to the Navier-Stokes Problem (equation (36)) is guaranteed for sufficiently small data.
To establish the existence and uniqueness of solutions to the weak formulation (Equation (73)), Galerkin’s method which involves constructing solutions of certain finite dimensional approximations to equation (73) and then passing to limits can be used [43] in which the procedure involves constructing approximate solutions, deriving energy estimates for the approximate solutions, and then convergence of approximate solutions to a solution.
2.9. Numerical Formulations
A numerical model for the groundwater pollute transport problem is derived by discretizing the governing PDE’s by the finite element method (FEM) and in time using the implicit Euler finite difference method.
2.10. Approximation with Mixed Finite Element Method
The ability to use a mesh of finite elements to accurately discretize a domain of any size or shape and the possibility of constructing higher order approximating polynomials for greater accuracy makes the finite element method a powerful tool in computational fluid dynamics. FEM is a numerical technique for finding approximate the solutions and by
where and are shape functions. It also uses the idea that integrals in the weak form (equation (65)) can be broken up into a sum of integrals over an arbitrary collection of disjoint subdomains whose union is the original domain. This implies that the problem can be treated locally and then assembled to the global system.
Below is an algorithm for the implementation of the FEM used in the present work.
| ||||||||||||||
The discrete weak formulation is now obtained below.
Assume is polygonal so that it can be tessellated with a set of triangles. We now define a triangulation.
Definition 11. Let be polygonal. A triangular mesh or a triangulation of is a set of (nondegenerate) simplices such that
and is the maximum diameter of triangle [39] since the spaces , and are infinite dimensional; the approximation following finite dimensional subspace is in
which depend on the discretization parameter . The discrete weak formulation for equation (83) is formulated as follows:
Find ,
such that
Defining bilinear forms
By the expressions in (87)
and the trilinear form
by
equations (84) can then be expressed in
If is a set of vector-valued-basis functions, and and scaler basis functions that span the finite dimensional space , and , respectively, such that then, the discrete velocity pressure and concentration components are written as a linear combination of their respective basis functions as follows: where , and are the number of unknowns for velocity, pressure, and concentration, respectively, is the number of nodal points where velocity or concentration is known on the boundary , and and are scalar coefficients for the velocity and pressure basis functions, respectively, while is the time-dependent coefficient for the concentration basis function. Substituting the expansions in equations (94) and (95) into equations (84), together with and , we obtain equation (97) as and
Then, equations (97) and (98) can be rewritten in or in block matrix form as in
where for and .
The right-hand vectors are expressed in
2.11. Discretization of the Diffusion-Convection Equation
The diffusion-convection equation is discretized in space using finite elements and in time using an implicit finite difference scheme. In this scheme, a system of linear equations is solved to calculate the values as functions of the previous values at each time step. Substituting the expansion in equation (96) into equation (84) together with the following equation (96) is obtained:
Next, discretize in time by substituting for in equation (106), which gives
Equation (108) can be rewritten as shown or in a matrix form as
where in
for
The right-hand vector in (112) is expressed as
The Navier-Stokes system of equation (100) is first solved, and the obtained solution is utilized to solve the system of the diffusion-convection equation (110). A potential source of instability of solutions to the diffusion-convection equation is now discussed. Methods that are applied to solve problems where no convection is present may totally fail when applied to convection-dominated problems. [44]. For instance, when the diffusivity coefficient of the diffusion-convection equation is very small, the diffusion-convection equation becomes convective dominated and for such a case, an artificial diffusive term is added through a process known as upwinding to obtain stable solutions.
2.12. Streamline Upwind Petrov-Galerkin Method
Solutions to the convection-diffusion equation in equation (66) may develop spurious oscillations (if convection dominated) unless the exact solution happens to be globally smooth [45]. An in-depth analysis on the deficiencies of the classical Galerkin approach in obtaining solutions of convection-dominated transport problems is discussed in [46] The finite element method requires the addition of and an extra term that balances the diffusive terms to obtain stable numerical solutions. A popular and efficient remedy is to augment the Galerkin finite element method formulation of the diffusion-convection equation in equation (83) by an extra term that adds artificial diffusion. This extra term is evaluated over the interior of and is a function of the residual in (113) of the PDE in
This consistently stabilized method can take the form of
where is the interior of , the number of interior points in , and is the residual of the diffusion-convection equation (114) and defined in
is an operator applied to the test function and is the stabilization parameter and defined as the artificial diffusion coefficient in equation (117). The choice of in equation (115), describes what is known as the streamline upwind Petrov-Galerkin (SPUG) finite element method.
The word streamline implies that the stabilization parameter is selected in such a way that the artificial diffusion term added is in the direction of flow and not perpendicular to it because convective transport occurs along the streamlines, [40]. The idea of adding diffusion along the flow lines helps to avoid what is called crosswind diffusion [46]. The development of a general theory for selecting optimally is still an area of research. In equation (115), the stabilization term contains a parameter whose value significantly influences the quality of the numerical solution. An optimal choice of for convection-dominated equations is usually not known because most of the existing procedures are based on somewhat heuristic procedures, [47]. Detailed discussions on these stabilization schemes are provided [24, 44, 46, 47–54].
2.13. Choice of the Finite Element Spaces , , and
For stability, the velocity is approximated using piecewise quadratic functions while the pressure is estimated using piecewise linear functions [41]. Concentration is also approximated using piecewise linear functions. The discrete spaces are defined in and in
Approximating polynomials are chosen in such a way that the velocity is computed at the vertices of the triangle and the midpoints of the edges, while pressure and concentration are computed at the triangle vertices only as shown in Figures 4 and 5. The polynomials are built such that they are described on one node and vanish with the others. The quantity of nodes related to a triangular detail is identical to the quantity of nearby ranges of freedom. Convergence houses are ruled with the aid of using the steadiness concerns expressed in the ellipticity requirement of the inf-sup circumstances of Brezzi and Babuska that is a generalization of the coercivity circumstance implied in the Lax-Milgram theorem [41]. For the Navier-Stokes equations, one such preference of a strong and convergent family of finite detail spaces used for approximating velocity and pressure, which can be inf-sup stable is the Taylor-Hood detail used in this study work, [40, 41].
2.14. Solving the Discrete Algebraic System
The algebraic system in equation (117) is nonlinear and hence requires an iterative method with a linearized problem being solved at each step. Starting with an appropriately chosen initial guess (the solution of the Stokes flow), the next guess is obtained by solving the resulting linear system. The process is repeated until the difference between the solutions obtained at iteration step and is greater than a given value of tolerance. A summary of the above steps can be seen in the algorithm (Figure 6).
Similarly, using an initial guess for the diffusion-convection equation, we solve for in
The above algorithms are implemented in FreeFem and some MATLAB codes in addition are utilized.
3. Results and Discussions
Consider a hypothetical factory located next to a water body as shown in Figure 7. The factory discharges the pollution on part of the land marked . Groundwater flows in the direction to where is the boundary between the water body and the groundwater domain while is the downstream boundary.
A 2-D cut of in the plane where represents the horizontal section of and represents the vertical section of is presented. Figure 8 shows a discretized domain of size to and to Circular obstacles are inserted to simulate a real groundwater flow process. The domain is now described in detail as follows: (a)is the inlet boundary(b)is the upstream boundary with a baseline concentration of pollute . is the interface between the water body and the computational domain (c)is the downstream boundary that is far away such that it can be assumed that the flux is zero(d)The flux on is proportional to the difference between the concentration of , and the concentration in , (e)Diffusivity coefficient is a constant
The following boundary and initial conditions for the diffusion-convection equation (121) are prescribed in and the following boundary conditions. For velocity in the Navier-Stokes equations are prescribed in expressions
In addition, no-slip boundary conditions are prescribed on the obstacles in and no boundary condition in prescribed on which implies a homogenous Neumann condition on In the next section, the Navier-Stokes and diffusion-convection equations are solved to determine the impact of various flow and concentration parameters.
3.1. Velocity and Pressure Profile
For Reynolds’ number, , the velocity profile shown in Figure 9 shows that the maximum velocity of 15.3 occurs along the line , which is zero for and . This is because of the prescribed inflow velocity, which is parabolic and the no-slip conditions imposed on the boundaries. Data in Figure 9 also shows that the velocity increases in regions where the obstacles are widely spaced. The color bar shows the variations in the magnitude of velocity with the color variations from top to bottom in the order of decreasing magnitude. The pressure distribution shown in Figure 10 shows that pressure decreases in the direction of groundwater flow. It has a maximum of 282.36 at . It was reduced from 253.34 at to 166.29 at , further reduced to 50.123 at 9.6, and finally to 0 when . This is in agreement with the physical interpretation that a fluid flows from a region of higher pressure to a region of lower pressure. Comparing Figures 9 and 10, it is seen that the velocity is highest in regions of higher pressure and lowest in regions of lower pressure.
3.2. Concentration Profile
Prescribing the following boundary and initial conditions in expressions in
For the diffusion-convection equation, the transport equation is solved using the velocity obtained from solving the Navier-Stokes equation for . The concentration profile for diffusivity coefficient is studied. After 5 time steps, a steady-state solution is reached. Figure 11 shows its concentration profile. The concentration on reduces from 100 to 6.3239 at between and . It further reduces to 0.81352 at between and . The baseline concentration reduces from 10 on to 6.03239 at and . It further reduces to 0.81325 at . The concentration elsewhere is 0.
3.3. Impact of Increasing Inflow and Baseline Concentrations
Keeping the baseline concentration constant and increasing the inflow concentration, the transport equation is solved for the boundary and initial conditions prescribed by the expressions in and. At steady state, the concentration profile is shown in Figure 12.
Diffusion of molecules depends entirely on movement from regions of high concentration to regions of low concentration. That is, diffusion occurs down the concentration gradient of that molecule. The greater the difference in concentration, the faster the molecule descends along the concentration gradient. If the difference in concentration is not large, the molecules do not move quickly and the diffusion rate decreases. Figure 12 shows that the concentration on reduced from 500 to 31.524 at between and . It was then reduced to 3.9664 at between and . The baseline concentration reduces from 10 on to 3.9664 at and . The concentration elsewhere is 0. Figures 12 and 13 show that increasing the inflow concentration while keeping the baseline concentration constant increases the amount of the pollutant diffusing from and reduces the amount of pollution due to convective transport from Keeping the inflow concentration constant and increasing the baseline concentration, the transport equation is solved for the boundary and initial conditions prescribed by expressions in and . At steady state, the concentration profile is shown in Figure 14.
Figure 14 shows that the concentration on reduces from 100 to 6.03239 at between and . It then reduces to 0.81352 at between and . The baseline concentration reduces from 50 on to 11.834 at . It then reduces to 0.81352 at and . Elsewhere, the concentration is 0. Figures 14 and 15 show that increasing the baseline concentration while keeping the inflow concentration constant increases the amount of pollutant due to convective transport from . The amount of pollutant diffusing from is constant. This is because the velocity of groundwater at the boundary is zero.
3.4. Impact of Decreasing the Diffusivity Coefficient
Using the boundary and initial conditions prescribed in equation (124) and for the Navier-Stokes equations, the transport equation is solved for different diffusivity coefficients.
Figure 16 shows the steady state concentration for described earlier. Figure 17 shows the steady state concentration for . The concentration on reduces from 100 to 6.5495 at between and . It then reduced to 1.0524 at between and . The baseline concentration reduces from 10 on to 6.5495 at and . It was reduced to 1.0524 at . The concentration elsewhere is 0.
Figure 18 shows the steady state concentration profile for . The concentration on reduces from 100 to 3.8776 at between and , and the baseline concentration reduces from 10 to to 3.8776 at and . The concentration elsewhere is 0. Figure 19 shows the concentration profile for , the concentration on reduced from 100 to 1.8236 at between and , and the baseline concentration reduces from 10 on to 1.8236 at . The concentration elsewhere is 0. From Figures 16–19, it is seen that decreasing the diffusivity coefficient increases the number of pollutants due to convective transport from the boundary and decreases the amount of pollutant that diffuses from the boundary
3.5. Comparison of Various Scenarios for the Management of Polluted Computational Domains for Pollutant Transport Problems and Their Impact on Groundwater Quality
In this study, we present simulation results for four hypothetical case studies. Figure 11 shows the results of the concentration profiles and diffusion coefficients (initial and final concentrations in this scenario) of the aquifer for four cases, including the scenario investigated. The velocity and pressure profile [34] concentration profiles for increasing input and base concentrations, decreasing diffusion coefficient and concentration profile for diffusivity coefficient . For the three cases, the movement of contaminants can be expressed as an increase in the input concentration, and while maintaining a constant reference concentration, the number of contaminants diffused from increases, and the number of contaminants due to convective movement from decreases. In addition, it can be seen that the decrease in diffusion coefficient increases the number of pollutants due to convection motion from the boundary and decreases the number of pollutants diffused from the boundary . Concentrations elsewhere are zero, indicating the effect of boundary conditions on groundwater contamination. The same result as above can be achieved by modifying the hypothetical scenario for the contamination site in terms of different concentrations and diffusion coefficients, probing the location of the contaminated site, and keeping 5-step scenario. For the hypothetical scenario, the results show that the groundwater quality can be controlled by the dividing wall in shallow aquifers, but not in deep aquifers. However, lining contaminated areas is an effective way to prevent groundwater contamination in both shallow and deep aquifers [32]. In addition, the use of linings in polluted areas is superior to other methods of collecting and treating wetland wastewater to protect aquifers from contamination. This method can be used taking into account the conditions of the geotechnical characteristics of the soil, the type of lining, the pollution load in these areas, and the type of element.
4. Conclusions and Recommendations
The idea of a depth-dependent initial filtration coefficient is consistent with the following conceptual model: Given a constant flow rate, the larger particles that are more likely to be trapped are more trapped in the porous medium, so most of the larger particles are retained at the top of the conceptual domain as relatively small particles, with the small particles being less likely to be trapped. This disproportionate absorption leads to an increase in the concentration of fine particles in the migratory particle population and the deeper parts of the porous medium have a greater ability to accommodate the particles, but the colloidal population becomes less sticky as the migration distance increases. Finally, particles that are less likely to be captured are more likely to be found in wastewater by transporting longer filters or passing through the medium. As a result, the initial filtration coefficient decreases along the porous medium. Therefore, groundwater is not entirely harmless, like the simulation results in Figures 8–19. An increase in the input and base concentrations, as well as a decrease in the diffusivity coefficient , increases the extent of groundwater contamination. Building a model that takes into account the processes involved in the transport of soil pollutants is complicated. In this model, for example, soil porosity was treated as a constant to facilitate calculation, analysis, and programming, but it is actually a highly variable quantity. This model can be improved by taking into account a 3D transitory equation from Navier-Stokes, adding a variable coefficient of porosity and diffusivity of the soil, increasing the dimensions of the calculation area and the number of obstacles to obtain a better approximation of the geometric structure of the water underground. Industries or others that dispose of toxic substances, such as chemicals, oil, or pharmaceuticals, must treat them before disposal to reduce the level of groundwater contamination. Recommendations of dumping sites should be located far away from water bodies or direction of groundwater flow. To further check the presented simulation results, experimental studies will be conducted.
Data Availability
No raw data were used in support of this study
Disclosure
The authors acknowledge that this paper was partly based on a study by Isaac Enyogoi (IE) (2017).
Conflicts of Interest
The authors declare no conflict of interest
Acknowledgments
The authors appreciate the following various departments: civil and environmental engineering, mathematics at Kyambogo University and Makerere University for the contributions made in this write-up, the International Science Programme (ISP) for contributing to the development of active and sustainable environments for higher education especially in field of applied mathematics. Special thanks to Jacob Nyende (JN) (PhD)-(PI), Henry Kasumba (HK) (PhD), Savannah Nuwagaba (SN) (PhD), William Wanasolo (WW) (PhD), and John Mango (JM) (Prof., PhD) for the guidance and commitment in this research study.