Experiments and numerical simulations on transport of dissolved pollutants around spur dike

Li-ping CHEN*, Jun-cheng JIANG

College of Urban Construction and Safety Engineering, Nanjing University of Technology,Nanjing 210009, P. R. China

1 Introduction

A spur dike may be defined as a structure extending outward from the bank of a stream for the purpose of deflecting the current away from the bank to protect it from erosion. Flow separation around a spur dike and circumfluence downstream occur because the streamline is compressed and velocity increases (Duan et al. 2010; Cui et al. 2008; Hua et al. 2008).Pollutant transport around the spur dike is highly three-dimensional. There is no doubt that a three-dimensional hydrodynamic water quality model is needed to simulate pollutant transport around a spur dike. Mayerle et al. (1995)and Tang et al. (2007)assumed hydrostatic pressure distribution to simulate flow in the vicinity of a spur dike. The hydrostatic pressure distribution assumption is not reliable for computation of flow fields around a spur dike because the water surface curve around a spur dike changes significantly. Wu and Yuan (2007)and Ai and Jin(2008)developed non-hydrostatic models for free surface flows. In computational fluid dynamics, the VOF (volume of fluid)method is a numerical technique for tracking and locating the free surface. The method is based on the idea of the so-called fraction function. The equation of fraction function cannot be directly solved easily, since the fraction function is discontinuous. Nevertheless, some attempts have been made. The most popular approach to the equation of fraction function is the so-called geometrical reconstruction, originating in the works of Hirt and Nichols (1981). Approaches to the geometrical reconstruction have developed rapidly (Ubbink and Issa 1999; Sussman and Puckett 2000; Francois et al. 2006;Zou and Zheng 2008; Caboussat et al. 2008; Dolbow et al. 2008; Afkhami et al. 2009; Raessi et al. 2010; Larmaei and Mahdi 2010). The VOF method was used in the simulation of a spillway tunnel (Sha et al. 2006). Fu et al. (2006)used the VOF method to calculate the flow field of Acipenser sinensis’s spawning site downstream of the Gezhouba Dam in China. A three-dimensional large eddy simulation model of unsteady hydraulic jumps and pressure fluctuations was established by adopting an improved VOF free-surface tracking method (Qin et al. 2010).

In the VOF method, the geometrical reconstruction is inefficient. The basic differencing schemes are not good for interface flow. The reason is that the low-order scheme causes serious numerical diffusion around the interface and the high-order scheme causes unphysical oscillations. The normalized variable diagram (NVD)bounded differencing scheme was introduced by Leonard (1988). Based on NVD, the switch between the second-order upwind and the central and first-order upwind differencing (SOUCUP)(Zhu and Rodi 1991)and the combination of a second- and third-order interpolation profile applied in the context of the normalized variable formulation (STOIC)(Darwish 1993)differencing schemes were proposed. In these differencing schemes, the switch between the basic schemes causes an unsteady-state solution. Jasak et al. (1999)modified the normalized variable in terms of gradients of the dependent variable to propose the Gamma bounded differencing scheme. In the Gamma scheme, a smooth transition between upwind differencing and central differencing is realized, but the calculations for gradients of the dependent variable increase the computational cost.

In this study, the compressive VOF scheme was developed by introducing an artificial compressive term into the VOF equation. The derivation of the three-dimensional transient compressive pollutant transport model (CPTM)is based on the compressive VOF scheme. The difference between the CPTM and the general water quality models is that the CPTM takes into account the effect of the water-air interface on pollutant transport. Based on the NVD and the cubic equation (CE), this study developed the CE bounded differencing scheme. Without reconstructing the interface, the computational efficiency of the compressive VOF with the CE differencing scheme is very high. For the calibration and validation of the CPTM, laboratory experiments were carried out in a straight rectangular flume with a non-submerged spur dike.On the basis of comparison between the simulations and experiments, the flow structure and pollutant transport features around a non-submerged spur dike were obtained.

2 Model development

2.1 Continuity equation and momentum equation of computational cell

The VOF method offers ideas that can be used in two-phase flows. The fraction of volume occupied by water in a computational cell is usually denoted by α. If α=1, the liquid phase, i.e. water, occupies the computational cell. If α=0, the gas phase, i.e. air,occupies the computational cell. If 0＜α＜1, the computational cell is partially filled and has an air-water interface. In each computational cell, the volume fractions of all phases sum to unity. If β=1-α, the density ρ and velocity uiof the computational cell are defined as

where ραand ρβare the densities of water and air, respectively; uαiand uβiare the velocities of water and air, respectively; and the subscript i equals to 1, 2, or 3, representing the three directions in the Cartesian coordinate system.

The continuity equation of the computational cell is calculated according to the mass conservation theorem:

Utilizing Eq. (1)and Eq. (2), Eq. (3)reduces to

The momentum equations of the liquid phase and the gas phase in a computational cell are written as follows:

where μαand μβare the turbulent dynamic viscosities of water and air, respectively; pαand pβare the pressures of water and air, respectively; and giis the gravity acceleration component. It is assumed that

where p is the pressure of the computational cell.

Eq. (5)is added to Eq. (6), and the momentum equation of the computational cell is formulated as

2.2 Compressive VOF model equations

According to the mass conservation theorem, the differential equations of α in tensor form can be written as follows:

Utilizing Eq. (1)and Eq. (2), Eq. (10)becomes:

2.3 Transport equation of dissolved pollutant

The transport equation for a dissolved pollutant in water is

where C is the concentration of the dissolved pollutant, SCis a source term for the biochemical process, Dmis the turbulent dispersion coefficient, and σCis the turbulent Schmidt number relating the turbulent diffusivity of the dissolved pollutant to the eddy viscosity νt. A value of σC=1 was used in this study.

Utilizing Eq. (1)and Eq. (2), Eq. (12)is written as follows:

Eq. (4), Eq. (9), Eq. (11), and Eq. (13)constitute the system of equations describing dissolved pollutant transport.

Turbulence is accounted for by adopting the standard k-ε model (Markatos 1986). The governing equations for the turbulent kinetic energy k and the dissipation rate ε are

3 CE differencing scheme

The finite volume method (FVM)of discretization was applied. Free surface flow and the transport of dissolved pollutants are bounded; that is, the water and air at the free surface are immiscible and the dissolved pollutant in the water cannot pass through the water-air interface. Therefore, it is necessary to consider the boundedness of free surface flow and dissolved pollutant transport. The bounded differencing scheme based on CE was developed in this study.

Fig. 1 Differencing schemes in NVD

Then,

The role of the K coefficient can be seen in Fig. 1(b). The larger the value of K, the more blending will be introduced. At the same time, the transition between UD and CD will be smoother. An upper limit on K comes from the accuracy requirement. For good resolution, K should ideally be kept within the range of 0.1 ≤K ≤0.5.

4 Numerical simulation and experiment

4.1 Physical model

Simulations and a series of experiments were performed on a non-submerged spur dike in a laboratory flume. The flume was 0.4 m wide, 0.4 m deep, 14.0 m long, and had a 3 700 bed slope. The mean flow depth was 0.087 m. The spur dike was rectangular, 0.1 m long, 0.1 m high, and 0.015 m wide, and was installed in the middle of the flume. Spur dike models were angled at 60ο, 90ο, and 120οupstream in the experiments.

The computational zone was 3.0 m long, 0.4 m wide, and 0.2 m high. The coordinates of the computational and experimental zones are shown in Fig. 2. The bottom of the flume was at z = 0. The dumping point was at x = 0.1 m and y = 0.05 m on the free surface. The sampling cross sections were at x = 1.12 m and x = 1.27 m. The sampling sites at each cross section were y = 0.01 m, y = 0.09 m, and y = 0.17 m, respectively. The heights of sampling sites were all z = 0.04 m. A laser doppler velocimeter was used to measure velocity.

Fig. 2 Coordinates of computational and experimental zones

4.2 Conditions for simulations

The length and height of the computational zone were set at 3.0 m and 0.2 m, respectively.Through multiple computations, it was found that the length of the vortex behind the spur dike increased only 0.2%, while there was a 40% increase in the length of the computational zone.The increase of the height of the computational zone had no influence on the water surface curve.

The computational zone was anomalous because of the existence of the spur dike. In the mesh generation, the computational zone was divided into five blocks as shown in Fig. 3, and the structure grid was adopted. The numbers of nodes in x, y, and z directions were 120, 27, and 40, respectively.

Fig. 3 Mesh around spur dike

As for boundary conditions, the top and bottom of the computational zone were regarded as atmospheric pressure and a no-slip boundary, respectively. The inlet boundary was composed of the water flow inlet and air inlet. The water discharge was 0.003 6 m3/s at the inlet. The air velocity at the inlet was set at zero. The k and ε at the inlet boundary were evaluated according to empirical formulas. The gradients of all physical quantities at outlet cross sections were assumed to be zero.

The relative tolerance limits the relative improvement from initial to final solution. In the experiments in this study, the relative tolerances for all physical quantities were set at 10-5to force the solution to converge.

4.3 Results and analysis

For the physical model, the computer used in this study was the dual-core Intel Pentium 4 CPU (3.00 GHz). The computational times for one second of flow with the CE differencing scheme and the Gamma differencing scheme are, respectively, 568 s and 694 s. The computational efficiency of the CE differencing scheme is higher than for the Gamma differencing scheme.

The calculated velocity fields at plane z = 0.02 m around spur dikes angled at 60ο(Case 1),90ο(Case 2)and 120ο(Case 3)are shown in Fig. 4. It can be seen that there is a large vortex behind each spur dike.

Fig. 4 Velocity distribution at plane z = 0.02 m in three cases

A comparison of the vortex lengths behind the spur dike in three cases is shown in Table 1.In Case 2, the vortex length caused by the spur dike at a 90οangle is the longest. In Case 3,water flow around the spur dike with a 120οangle forms the shortest vortex.

Table 1 Comparison of vortex lengths behind spur dike

The small relative errors between the calculated and experimental data indicate that the models proposed in this study are applicable to the flow around spur dikes.

The water level around the spur dike in Case 2 for longitudinal section y = 0.05 m is shown in Fig. 5. The simulation is in good agreement with the experiment. The water level around the spur dike changes severely, and the free surface at the tip of the spur dike rises and then drops sharply.

At the cross section around the spur dike, the water level of the main flow is higher than that of the back flow. The water surface curve at the cross section x = 1.05 m in Case 2 is shown in Fig. 6. The water surface opposite of the spur dike is higher than that on the same side of the spur dike. Similarly, free surfaces around the spur dike in Case 1 and Case 3 have the same rule as in Case 2. The drops of the water level along the longitudinal profile y = 0.05 m and the maximum differences of the water level at the cross section x = 1.05 m in three cases are compared in Table 2. The changes of the water level in Case 2 are more obvious than in Case 1 and Case 3.

Fig. 5 Water surface curve along longitudinal profile y = 0.05 m

Table 2 Water level comparison

For water flow, there are three reasons for errors between calculations and experiments.The first one is that water occasionally flows over the brim and above the spur dike in experiments. However, there is no water above the spur dike in these calculations. The second one is that the computational zone around the spur dike is anomalous. Thirdly, coarse meshes in block 1 and block 2 cause errors. Subdivision of meshes around the spur dike can reduce error but increases computational cost.

Fig. 7 shows the secondary vortical flow at a cross section (x = 1.05 m)in three cases.The secondary vortical flow in Case 2 is the most significant because the streamline is severely compressed.

Fig. 7 Secondary vortical flow at cross section x = 1.05 m

4.4 Results and analysis of dissolved pollutant transport

The dissolved pollutant transport experiments were conducted under the same flow conditions as the flow velocity measurement experiments. A quantity of 50 mL of dissolved pollutant (rhodamine B)was instantaneously dumped into the flume with an initial concentration of C0= 2 g/L . A spectral photometer was used to measure the pollutant concentration.

Fig. 8 shows the results of the simulation of rhodamine B concentration transport around the spur dike in Case 1. Case 2 and Case 3 have similar patterns of concentration transport.

Fig. 8 Simulation of rhodamine B concentration transport around spur dike in Case 1

Fig. 8 shows that the front edge of the rhodamine B cloud arrives at the tip of the spur dike at two seconds and then flows around the spur dike following the main flow. There is a concentration tail behind the rhodamine B cloud. Two concentration centers arise downstream of the spur dike and a dead water zone appears at the fourth second. At the same time, the concentration tail stretches. The concentration center in the dead water zone disappears quickly because of exchange with the main flow. Between the third and sixth seconds the shape of the pollutant cloud changes continually. At the sixth second, the pollutant cloud passes through the spur dike and enters the circumfluence zone. The pollutant remains in the circumfluence zone for a long time due to the weak exchange with the main flow.

Fig. 9 is the pollutant concentration as a function of time at the cross section x0L0=12 10 in three cases, where x0is the distance between the sampling transection and the spur dike,and L0is the length of the spur dike.

As shown in Fig. 9, the pollutant concentration at y L0= 17 10 in Case 1 is lower than in Case 2 and Case 3 in the first ten seconds, because this sampling site is at the edge of the pollutant cloud. After the concentration center leaves the x0L0= 12 10 cross section, the concentration at y L0= 17 10 decays sharply in all three cases. The concentration at y L0= 110 increases slowly and then decreases. In Case 1, the concentration at y L0=9 10 decays from the maximum value and then rises significantly because of back flow. According to the rise of concentration at y L0= 9 10 in Fig. 9, less pollutant with back flow passes through y L0= 9 10 with the increase of the spur dike angle.

Fig. 9 Rhodamine B concentration as function of time at cross section x0 L0 =12 10

During dissolved pollutant transport, the errors between experiments and calculations are caused by dumping behavior. In calculations, the pollutant has no velocity relative to the water.However, the pollutant has velocity relative to the water for the sudden dumping in experiments. A new device for dumping should be designed in future studies to reduce error.

Fig. 10 is the pollutant concentration as a function of time at the cross section x0L0= 27 10 in three cases. After falling from the maximum concentration, as the spur dike angle increases, the concentration fluctuation at y L0= 17 10 tends to increase in the experiments. The reason for concentration fluctuation is that the pollutant carried by the small corner vortex is discontinuous. The maximum concentration at y L0= 110 in Case 3 is higher than in other cases. The difference between concentrations at y L0= 110 and at y L0= 9 10 increases with the spur dike angle. In Case 1, the maximum concentration is slightly higher at y L0= 9 10 than at y L0= 110. However, in Case 2, the maximum concentration at y L0= 9 10 is less than at y L0= 110. In Case 3, the maximum concentration at y L0= 9 10 is much less than at y L0= 110. Before the pollutant carried by back flow passes through y L0= 9 10, much pollutant carried by back flow enhances the concentration at y L0= 110 in Case 3.

5 Conclusions

Based on the compressive VOF scheme, the water quality model CPTM was developed.Without reconstructing the interface, the compressive VOF scheme with the bounded differencing scheme CE can effectively simulate the free surface and highly three-dimensional flow around a spur dike. Spur dikes with different angles create different flow patterns. The vortex length caused by a spur dike with a 90οangle is the longest. Water flow around the spur dike with a 120οangle forms the shortest vortex. The water level shows a steep change around the spur dike.

The different three-dimensional flow patterns around the spur dike determine the different transport features of dissolved pollutants. The decay of the pollutant is highest in the main flow. The rise and decay rates of the concentration in the circumfluence zone are proportional to the velocity of back flow.

The concentration of fluctuation behind the spur dike becomes intensive with the increase of the spur dike angle. The reason for concentration fluctuation is that the pollutant carried by the small corner vortex is discontinuous. It is concluded that the small corner vortex increases with the spur dike angle.

Afkhami, S., Zaleski, S., and Bussmann, M. 2009. A mesh-dependent model for applying dynamic contact angles to VOF simulations. Journal of Computational Physics, 228(15), 5370-5389. [doi:10.1016/j.jcp.2009.04.027]

Ai, C. F., and Jin, S. 2008. Three-dimensional non-hydrostatic model for free-surface flows with unstructured grid. Journal of Hydrodynamics, Ser. B, 20(1), 108-116. [doi:10.1016/S1001-6058(08)60035-9]

Caboussat, A., Francois, M. M., and Glowinski, R. 2008. A numerical method for interface reconstruction of triple points within a volume tracking algorithm. Mathematical and Computer Modelling, 48(11-12),1957-1971. [doi:10.1016/j.mcm.2008.05.009]

Cui, Z. F., Zhang, X. F., and Feng, X. X. 2008. Numerical simulation on scour around spur-dike by 3D turbulent model. Journal of Hydrodynamics, Ser. A, 23(1), 33-41. (in Chinese)

Darwish, M. S. 1993. A new high-resolution scheme based on the normalize variable formulation. Numerical Heat Transfer, Part B, 24(3), 353-371. [doi:10.1080/10407799308955898]

Dolbow, J., Mosso, S., Robbins, J. and Voth, T. 2008. Coupling volume-of-fluid based interface reconstructions with the extended finite element method. Computer Methods in Applied Mechanics and Engineering, 197(5), 439-447. [doi:10.1016/j.cma. 2007.08.010]

Duan, J. G., He, L., Fu, X. D., and Wang, Q. Q. 2010. Mean flow and turbulence around experimental spur dike. Advances in Water Resources, 32(12), 1717-1725. [doi:10.1016/j.advwatres.2009.09.004]

Francois, M. M., Cummins, S. J., Dendy, E. D., Kothe, D. B., Sicilian, J. M., and Williams, M. W. 2006. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. Journal of Computational Physics, 213(1), 141-173. [doi:10.1016/j.jcp.2005.08.004]

Fu, X. L., Li, D. M., and Chen, Y. B. 2006. Flow field calculation and analysis of Acipenser Sinensis’s spawning site in Gezhouba downstream. Advances in Water Science, 17(5), 700-704. (in Chinese)

Hirt, C. W., and Nichols, B. D. 1981. Volume of fluid (VOF)method for the dynamics of free boundaries.Journal of Computational Physics, 39(1), 201-225. [doi:10.1016/0021-9991(81)90145-5]

Hua, Z. L., Xing, L. H., and Gu, L. 2008. Application of a modified quick scheme to depthaveraged k-ε turbulence model based on unstructured grids. Journal of Hydrodynamics, Ser. B, 20(4), 514-523. [doi:10.1016/S1001- 6058(08)60088-8]

Jasak, H., Weller, H. G. and Gosman, A. D. 1999. High resolution NVD differencing scheme for arbitrarily unstructured meshes. International Journal for Numerical Methods in Fluids, 31(2), 431-449.[doi:10.1002/(SICI)1097-0363(19990930)31:2＜431::AID-FLD884＞3.3.CO;2-K]

Larmaei, M. M., and Mahdi, T. 2010. Simulation of shallow water waves using VOF method. Journal of Hydro-environment Research, 3(4), 208-214. [doi:10.1016/j.jher.2009.10.010]

Launder, B. E., and Spalding, D. B. 1990. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering, 81(3), 269-289.

Leonard, B. P. 1988. Simple high-accuracy resolution program for convective modeling of discontinuities.International Journal for Numerical Methods in Fluids, 8(10), 1291-1318. [doi:10.1002/fld.1650081013]

Markatos, N. C. 1986. The mathematical modeling of turbulence flows. Applied Mathematical Modelling,10(3), 190-220. [doi:10.1016/0307-904X(86)90045-4]

Mayerle, R., Wang, S. S. Y., and Toro, F. M. 1995. Verification of a three-dimensional numerical model simulation of the flow in the vicinity of spur dikes. Journal of Hydraulic Research, 33(2), 243-256.[doi:10.1080/ 00221689509498673]

Qin, L., Tian, H. J., and Zhao, H. 2010. Large eddy simulation of flow fluctuating pressure based on an improved VOF free surface track method. Journal of Hydroelectric Engineering, 29(3), 92-96. (in Chinese)

Raessi, M., Mostagnimi, J., and Bussmann, M. 2010. A volume-of-fluid interfacial flow solver with advected normals. Computers and Fluids, 39(8), 1401-1410. [doi:10.1016/j.compfluid. 04.010]

Sha, H. F., Wu, S. Q., and Chen, Z. W. 2006. 3D numerical simulation for spillway tunnel. Advances in Water Science, 17(4), 507-511. (in Chinese)

Sussman, M., and Puckett, E. G. 2000. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows. Journal of Computational Physics, 162(2), 301-337.[doi:10.1006/jcph.2000.6537]

Tang, X. L., Ding, X., and Chen, Z. C. 2007. Experimental and numerical investigations on secondary flows and sedimentations behind a spur dike. Journal of Hydrodynamics, Ser. B, 19(1), 23-29. [doi:10.1016/S1001-6058(07)60023-7]

Ubbink, O., and Issa, R. I. 1999. A method for capturing sharp fluid interfaces on arbitrary meshes. Journal of Computational Physics, 153(1), 26-50. [doi:10.1006/jcph.1999.6276]

Wu, C. H., and Yuan, H. L. 2007. Efficient non-hydrostatic modelling of surface waves interacting with structures. Applied Mathematical Modelling, 31(4), 687-699. [doi:10.1016/j.apm.2005.12.002]

Zhu, J., and Rodi, W. 1991. Low dispersion and bounded convection scheme. Computer Methods in Applied Mechanics and Engineering, 92(1), 87-96. [doi:10.1016/0045-7825(91)90199-G]

Zou, J. F., and Zheng, Y. 2008. Application of bounded and compressed VOF method to interfacial flow.Journal of Zhejiang University (Engineering Science), 42(2), 253-258. (in Chinese)