Analysis of Nonlinear Dispersion of a Pollutant Ejected by an External Source into a Channel Flow

Abstract

This paper focuses on the transient analysis of nonlinear dispersion of a pollutant ejected by an external source into a laminar flow of an incompressible fluid in a channel. The influence of density variation with pollutant concentration is approximated according to the Boussinesq approximation, and the nonlinear governing equations of momentum and pollutant concentration are obtained. The problem is solved numerically using a semi-implicit finite difference method. Solutions are presented in graphical form and given in terms of fluid velocity, pollutant concentration, skin friction, and wall mass transfer rate for various parametric values. The model can be a useful tool for understanding the polluting situations of an improper discharge incident and evaluating the effects of decontaminating measures for the water body.

1. Introduction

Pollution of water sources like rivers and lakes, say resulting from industrial waste discharge, is a serious socioecological hazard. If the problem is not carefully controlled and monitored, large communities can be exposed to extensive health risks. Early detection of such pollution accidents, both in terms of extent and impact, is thus of major primary importance, and the subsequent requirement to take immediate corrective measure to redress the pollution problem and mitigate against its impact is equally so. In line with prevention better than cure adage, the development of accurate methods to predict the spread of a pollutant once a discharge has been detected and hence also the development of equally reliable preventive/corrective techniques is thus of paramount importance [1]. Spread of pollutants in a fluid flow depends largely on concentration coefficients [2]. These can be determined empirically for each type of pollutant. Investigations such as [3] and also the current one can help identify the pollutant physical properties (and the related mathematical parameters) likely to cause the greatest harm in the spreading the pollutant downstream. The importance of these kinds of investigations as well as the complimentary experimental works, say in large-scale water treatment and redistribution networks, thus makes them of great relevance [4–9].

It should however be remarked that to date, the literature on the transient analysis of problems consisting of buoyancy effects and nonlinear pollutant injection is still quite sparse. This then forms the crux of the current investigation. We suppose that a pollutant is introduced nonlinearly into a channel flow via an external source and proceed, under Boussinesq approximations, to investigate the transient diffusion and resultant spatial-temporal concentration.

In Section 2, the model problem is formulated and in particular the mathematical governing equations, initial and boundary conditions are outlined. We outline, in Section 3, the finite difference schemes that will be employed in the solution process, and graphical results and discussions follow in Section 4.

2. Formulation of the Problem

We consider a transient problem of fluid flow and nonlinear dispersion of pollutant in a rectangular channel as shown in Figure 1. In order to derive the governing equations, the following assumptions are made:(i)the fluid is viscous and incompressible;(ii)initially, the flow is fully developed through a rectangular channel;(iii)at time 𝑡>0, a given pollutant is injected into the flow from an external source; the viscosity of the fluid and the pollutant mass diffusivity then vary with its concentration;(iv)the influence of density variation with pollutant concentration has been considered only in the body-force term of the momentum equation and is approximated according to the Boussinesq approximation.

Figure 1 

Flow configuration and coordinate system.

Under laminar flow conditions, the problem is reduced mathematically to a transient coupled fluid flow and mass transfer problem given in one dimension as [1, 3, 8, 9] 𝜌𝜕𝑢𝜕𝑡𝜕=−𝑃𝜕𝑥+𝜕𝜕𝑦(𝜇(𝐶)𝜕𝑢𝜕𝑦)+𝜌𝑔𝛽(𝐶−𝐶₀),𝜕𝐶𝜕𝑡=𝜕𝜕𝑦(𝐷(𝐶)𝜕𝐶𝜕𝑦)+𝑆(𝐶), (2.1)(2.2) subject to the following initial and boundary conditions: 𝑢(𝑦,0)=𝑚(1−𝑦²𝑎²),𝐶(𝑦,0)=𝐶₀,𝜕𝑢𝜕𝑦(0,𝑡)=𝜕𝐶𝜕𝑦(0,𝑡)=0,for𝑡>0,𝑢(𝑎,𝑡)=0,𝐶(𝑎,𝑡)=𝐶_𝑤,for𝑡>0, (2.3)(2.4)(2.5) where 𝑢 is the fluid axial velocity, (𝑟,𝑧) are the radial and axial coordinates, respectively, 𝐶 is the pollutant concentration, 𝐶₀ is the pollutant reference concentration, 𝐶_𝑤 is the pollutant concentration at the walls, 𝑆 is the pollutant external source function, 𝑔 is the gravitational acceleration, 𝜌 is the density, 𝛽 is the concentration expansion coefficient, 𝑎 is the channel half width, and 𝑃 is the fluid pressure. As shown in (2.3), we employed a fully developed Poiseuille parabolic velocity profile as the initial condition for the flow. The pollutant concentration-dependent fluid dynamic viscosity 𝜇, mass diffusivity 𝐷, and external source are prescribed as follows: 𝜇=𝜇₀exp[𝑏₁(𝐶−𝐶₀)],𝐷=𝐷₀exp[𝑏₂(𝐶−𝐶₀)],𝑆=𝑄exp[𝑏₃(𝐶−𝐶₀)], (2.6) where 𝜇₀,𝐷₀,𝑏₁,𝑏₂, and 𝑏₃ are the viscosity coefficient, mass diffusivity coefficient, viscosity variation parameter, mass diffusivity variation parameter, and the pollutant external source variation parameter, respectively. The following dimensionless variables are introduced:𝑦=𝑦𝑎,𝑥=𝑥𝑎,𝑤=𝑢𝑎𝜈,𝜙=𝐶−𝐶₀𝐶_𝑤−C₀𝜈,𝑡=𝑡𝑎²,𝐷=𝐷𝐷₀,𝜇=𝜇𝜇₀,𝑃=𝑃𝑎²𝜌𝜈²,𝐾=−𝜕𝑃𝜇𝜕𝑥,𝜈=⁰𝜌,𝛼=𝑏₁(𝐶_𝑤−𝐶₀),𝛾=𝑏₂(𝐶_𝑤−𝐶₀),𝑛=𝑏₃(𝐶_𝑤−𝐶₀𝑎),𝜆=²𝑄(𝐶_𝑤−𝐶₀)𝜈,Gc=𝑔𝛽𝑎³(𝐶_𝑤−𝐶₀)𝜈²𝜈,Sc=𝐷₀. (2.7) Substituting (2.7) into (2.1)–(2.6), we obtain the following dimensionless coupled governing equations: Re𝜕𝑤𝜕𝜕𝑡=𝐾+𝜕𝑦(𝑒^𝛼𝜙𝜕𝑤𝜕𝑦)+Gc𝜙,𝜕𝜙=1𝜕𝑡𝜕Sc𝜕𝑦(𝑒^𝛾𝜙𝜕𝜙𝜕𝑦)+𝜆𝑒^𝑛𝜙,𝑤(𝑟,0)=𝑚(1−𝑦²),𝜙(𝑦,0)=0,𝜕𝑤𝜕𝑦(0,𝑡)=𝜕𝜙𝜕𝑦(0,𝑡)=0,for𝑡>0𝑤(1,𝑡)=0,𝜃(1,𝑡)=1,for𝑡>0, (2.8) where 𝜆 is the pollutant external source parameter, Gc is the solutal Grashof number, 𝐾 is the axial pressure gradient parameter, Sc is the Schmidt number, 𝛼 is the viscosity variation parameter, 𝛾 is the mass diffusivity variation parameter, and 𝑛 is the pollutant external source variation parameter. The dimensionless shear stress (𝐶_𝑓) and the mass transfer rate (Sh) at the channel wall are given by 𝐶_𝑓=−𝜕𝑤|𝜕𝑦_𝑦=1,Sh=−𝜕𝜙|𝜕𝑦₌₁. (2.9) In the following section, (2.8) is solved numerically and the skin-friction together with the wall mass transfer rate given in (2.9) are computed.

3. Numerical Solution

Our numerical algorithm is based on the semi-implicit finite difference scheme given in [10] for the isothermal viscoelastic case. As in [11], we extend the algorithm to the temperature equation and take the implicit terms at the intermediate time level (𝑁+𝜉) where 0≤𝜉≤1. The algorithm employed in [11] uses 𝜉=1/2; we will however follow the formulation in [10] and thus take 𝜉=1 in this article so that we can use larger time steps. The discretization of the governing equations is based on a linear Cartesian mesh and uniform grid on which finite-differences are taken. We approximate both the second and first spatial derivatives with second-order central differences. The equations corresponding to the first and last grid point are modified to incorporate the boundary conditions. The semi-implicit scheme for the velocity component reads (𝑤^(𝑁+1)−𝑤^(𝑁))Δ𝑡=𝐾+𝑒^𝛼𝜙𝑁𝑤^(𝑁+𝜉)_𝑦𝑦+[Gc𝜙+(𝑒^𝛼𝜙)_𝑦𝑤_𝑦]^(𝑁). (3.1) The equation for 𝑤^(𝑁+1) then becomes −𝑟₁𝑤^(𝑁+1)_𝑗−1+(1+2𝑟₁)𝑤_𝑗^(𝑁+1)−𝑟₁𝑤^(𝑁+1)_𝑗+1=explicitterms, (3.2) where 𝑟₁=exp(𝛼𝜙^𝑁)Δ𝑡/Δ𝑦². The solution procedure for 𝑤^(𝑁+1) thus reduces to inversion of tri-diagonal matrices which is an advantage over a full implicit scheme. The semi-implicit integration scheme for the concentration equation is similar to that for the velocity component. Unmixed second partial derivatives of the concentration are treated implicitly: Sc(𝜙^(𝑁+1)−𝜙^(𝑁))Δ𝑡=𝑒^𝛾𝜙𝑁𝜙^(𝑁+𝜉)_𝑦𝑦+[(𝑒^𝛾𝜙)_𝑦𝜙_𝑦+Sc𝜆𝑒^𝑛𝜙]^(𝑁) (3.3) The equation for 𝜙^(𝑁+1) thus becomes −𝑟₂𝜙^(𝑁+1)_𝑗−1+(Sc+2𝑟₂)𝜙_𝑗^(𝑁+1)−𝑟₂𝜙^(𝑁+1)_𝑗+1=explicitterms, (3.4) where 𝑟₂=exp(𝛾𝜙^𝑁)Δ𝑡/Δ𝑦². The solution procedure again reduces to inversion of tri-diagonal matrices. The schemes (3.2) and (3.4) were checked for consistency. For 𝜉=1, these are first-order accurate in time but second order in space. The schemes in [11] have 𝜉=1/2 which improves the accuracy in time to second order. We use 𝜉=1 here so that we are free to choose larger time steps and still converge to the steady solutions. The algorithm was also tested for both spatial and temporal convergences and shown to be independent of both mesh size and time step size.

4. Results and Discussion

Unless otherwise stated, we employ the parameter values: Gc=0.1,Sc=0.6,𝛼=0.1,𝛾=0.1,𝜆=0.5,𝑛=0.1,𝑚=0.1,𝐾=1,Δ𝑟=0.02,Δ𝑡=0.005, and 𝑡=50.

4.1. Transient Solutions

Figures 2 and 3 show the time development of the velocity and pollutant concentration profiles. The velocity profile eventually settles to a steady parabolic flow which is qualitatively similar to the initial profiles but quantitatively higher than the initial. This is expected since the injection of pollutant into the flow acts as a momentum source and hence naturally increases the flow velocity at subsequent times after the injection. Since the injection ceases at some point, we also naturally expect the velocities to eventually settle to steady state as illustrated in Figure 2.

Figure 2 

Developing velocity profiles.

Figure 3 

Developing concentration profiles.

The observed increase in the pollutant concentration with time shown in Figure 3 is a trivial consequence of the fact that pollutant is injected into the flow at time >0. The cessation of the injection process at some future time leads to the eventual steady concentration profiles displayed.

4.2. Dependence on Flow Parameters

As already noted, the pollutant injection is connected to momentum source terms and hence leads to increased steady flow velocities. Increase of the pollutant injection parameter, 𝜆, thus expectedly (indirectly) increases flow velocity as illustrated in Figure 4.

Figure 4 

Effects of $\lambda$ on velocity.

The increases pollutant concentration however also increases fluid viscosity and hence the competing effects of higher viscosity and higher pollutant injection result in the marginal velocity increases shown in Figure 4.

Figures 3 and 5 show that at a fixed value of parameter lambda, the pollutant concentration in the flow increases with time until it reaches a saturation state. Further injection of pollutant into the flow does not lead to any appreciable increase in saturation state concentration. We should however remark (as will be illustrated shortly) that large values of 𝜆 may not allow for the attainment of steady solutions but rather lead to finite time blowup phenomena.

Figure 5 

Effects of $\lambda$ on concentration.

For large values of 𝜆, the possibility of blowup phenomena is clearly illustrated in Figure 6. Physically, this results primarily from flow blockage due to excessive discharge of pollutant into the channel flow.

Figure 6 

Blowup of concentration, $n=1$.

In Figure 6, we compute and plot (for each value of 𝜆) the corresponding pollutant maximum concentration (𝜙_max) in steady state (𝑡=10) until a value of 𝜆 is reached (around 𝜆<0.4 for the current parameter values) at which no (steady) solutions exist but instead blowup phenomena is observed.

We next look at the effects of the solutal Grashof number Gc on both the flow velocity and pollutant concentration. Clearly Gc is connected to momentum source terms and thus flow velocities should directly increase with increasing Gc. This is illustrated in Figure 7.

Figure 7 

Effects of $\text{Gc}$ on velocity.

On the other hand, since the flow velocity, 𝑤, is uncoupled from the pollutant concentration, 𝜙, then changes in flow velocities are not expected to have any influence on the concentration. The solutal Grashof number only appears in the momentum equation and its increase will (even though it leads to increased velocities) thus not have any influence on the pollutant concentration as illustrated in Figure 8.

Figure 8 

Effects of $\text{Gc}$ on concentration.

The Schmidt number of various chemical substances which may act as pollutant varies depending on the chemical species. The values of Schmidt number (Sc) are chosen for hydrogen (Sc=0.22), water vapour (Sc=0.62), ammonia (Sc=0.78), and Propyl Benzene (Sc=2.62). The minor increases noted in the flow velocity, Figure 9 are, as explained before, due to competing effects of increased Boussinesq source terms and increased fluid viscosity, both due to increased Sc and hence increased 𝜙.

Figure 9 

Effects of $\text{Sc}$ on velocity.

Figure 9 shows that pollutant substances with higher Schmidt numbers lead to correspondingly higher pollutant concentrations within the flow.

Increasing the parameter 𝛼 leads to corresponding increases in the fluid viscosity and hence leads to reduced flow velocities. This is well illustrated in Figure 11.

As with the Grashof number, the fluid viscosity only plays a role in the momentum equations and since the velocity is uncoupled from the concentration equation, any changes in the viscosity (even though they result in changes in the velocity field) will have no effect on the pollutant concentration as shown in Figure 12.

The diffusion parameter 𝛾 manifests itself appreciably in the pollutant concentration but barely plays a role in the flow velocity as illustrated in Figures 13 and 14.

Due to limited fluid motion close to the wall, due to viscous effects, increased diffusion values would not be able to redistribute pollutant concentration around the flow hence the observed increases in pollutant concentration close to the walls. On the other hand, higher diffusion values within the main flow will help distribute the pollutant around the bulk flow hence leading to the reduced pollutant concentrations depicted in Figure 14.

The parameter 𝑛 plays an almost similar role to 𝜆 and hence the effect on the velocity and concentration with varying values of 𝑛 can be similarly explained. In particular, increasing 𝑛 directly increases the pollutant source strength and hence leads to increased pollutant concentration as shown in Figure 16. The resultant increase in 𝜙 leads to the aforementioned competing effects of increased Boussinesq source terms and increased fluid viscosity and hence changes (increase) in the velocity filed in response to changes (increase) in 𝑛 are quite marginal as illustrated in Figure 15.

4.3. Shear Stress and Mass Transfer Rate

According to Figures 4, 5, 13, and 14, the velocity increases with 𝜆 and hence the wall shear stress should also correspondingly increase with 𝜆 as shown in each graph of Figure 17. Similarly, the velocity decreases with increasing 𝛾 and hence the observed decrease in 𝐶_𝑓 with increasing 𝛾 also shown in Figure 17. A similar reasoning is applied to explain the features observed in Figure 18.

According to Figures 7, 8, 9, and 10, the velocity increases with Gc and hence the wall shear stress should correspondingly increase with Gc as shown in each graph of Figure 19. Similarly, the velocity increases (even though marginally) with increasing Sc and hence the observed increase in 𝐶_𝑓 with increasing Sc also shown in Figure 19. A similar reasoning is applied to explain the features observed in Figure 20.

Figure 10 

Effects of $\text{Sc}$ on concentration.

Figure 11 

Effects of $\alpha$ on velocity.

Figure 12 

Effects of $\alpha$ on concentration.

Figure 13 

Effects of $\gamma$ on velocity.

Figure 14 

Effects of $\gamma$ on concentration.

Figure 15 

Effects of $n$ on velocity.

Figure 16 

Effects of $n$ on concentration.

Figure 17 

Wall shear stress dependence on $\lambda$ and $\gamma$.

Figure 18 

Mass transfer dependence on $\lambda$ and $\gamma$.

Figure 19 

Wall shear stress dependence on $\text{Gc}$ and $\text{Sc}$.

Figure 20 

Mass transfer dependence on $\text{Gc}$ and $\text{Sc}$.

According to Figures 2, 3, 15, and 16, the velocity increases with 𝑡 until a steady state is reached and hence the wall shear stress should correspondingly increase with 𝑡 until also a steady state is reached when the time dependence ceases. This is shown in Figure 21. Similarly, the velocity increases (even though marginally) with increasing 𝑛 and hence the observed marginal increases in 𝐶_𝑓 with increasing 𝑛 also shown in Figure 21.

Figure 21 

Wall shear stress dependence on $t$ and $n$.

A similar reasoning is applied to explain the features observed in Figure 22 except in this case, based on the scalings on the vertical axis, the increase in Sh with increasing 𝑛 is more significant as compared to that for 𝐶_𝑓.

Figure 22 

Mass transfer dependence on $t$ and $n$.

5. Conclusion

We computationally investigate the transient dynamics of pollutant dispersion in a channel with a nonlinear waste discharge concentration subjected to one-dimensional gravity driven flow under Boussinesq approximations. We observe that there is a transient increase in the flow velocity with an increase in the pollutant injection. Larger increases in the injection are also shown to potentially lead to blockage due to excessive discharge of pollutant into the channel flow. We also similarly investigated the effects of the various other physical parameters on the flow velocity, the pollutant concentration, and hence also on the wall shear stress and wall mass transfer rates.

Acknowledgment

The authors would like to thank the National Research Foundation of South Africa Thuthuka programme for financial support.