Numerical modelling of pollution transport in Tom River.pdf Shallow water equations (SWE) are a widely used approach to modelling geophysical flows that combines acceptable computational cost and precision of the results obtained. It has wide application to modelling problems from fast-developing flows in dam breaks and tsunami waves to relatively slow flows in river estuaries and transport of sediments in lowland rivers. SWE [1] are obtained by integrating RANS equations by depth and correctly describe a flow in conditions of the depth being much less than the horizontal dimensions of an investigated area. In averaging, the distribution of pressure is hydrostatic, and flow characteristics are assumed to vary little with depth. In this approach, gravity, bottom and free surface shear stresses, and the Coriolis force are the main forces that drive the flow. Because of the assumptions above, SWE are widely used for modelling the atmosphere and the flow in natural basins. SWE are used to model flow in lakes and seas when the depth-dependent processes do not play a significant role. For example a series of works [2-4] illustrate the application of SWE to modelling tidal flow and tsunami waves in tidal regions of an ocean. Other works [5-7] show the application of this approach to modelling flows with significant wind shear on the surface of the Sea of Azov and the Caspian Sea. Applications of SWE to modelling flow in rivers and estuaries are extensive [8-10, 3, 11] and are not limited to two-dimensional formulations. 2D SWE are used to compute the flow in river sections that are very long (from dozens of meters to several kilometers) and have an irregular shape that significantly affects the flow [8, 11]. The geometry of the river is especially important in solving such problems as the transport of pollution in the flow, bed deformation [12, 13], and sediment transport [4]. Nevertheless, applications of 1 The reported study was funded by Russian Foundation for Basic Research under research project N 18-3100386. Numerical modelling of pollution transport in Tom River 49 2D equations to problems of bed deformation and sediment transport are not common due to the very large time scales of these processes. Two of the approaches to modelling processes in a large spatial or time scale are 1D formulations [14-16] and combinations of 1D and 2D models [17, 18] in which 1D computation of the whole river is adjoined with precise 2D computations of the areas of particular interest. 1D solvers can provide accurate information in terms of total sediment load, though they fail at predicting the bed erosion and calculating the final shape of the cross section. One of the methods to construct the 1D river stream model from 2D depth-averaged equations is integrating them by river width, which gives equations an extra integral pa-rameter-an area of water section. Another approach is omitting the terms that describe the variation of the parameters with one of the spatial coordinates. 1D models are also used when variation of the flow parameters from the river width is negligible compared to longitudinal variation. Very fast-developing flows (such as dam-break problems) or flows in long industrial channels with smooth walls and a flat bottom are examples of flows with such parameters [19]. Solving the fully 3D flow field with specific treatment of the free surface and bottom boundary condition is limited to short river sections adjacent to hydraulic facilities and small-scale flows in laboratory statements [9, 20-22]. In this work, 2D SWE are considered and used as a preferable approach that combines acceptable computational cost for computing long river sections (dozens of kilometers) with the precision of results obtained that is essential for solving such problems as evaluating the anthropogenic changes in the river bed, modelling ice movement and local flooding in spring, and detailing the distribution of pollutants discharged into the river with wastewater. In many cases, the turbulence of the flow is not accounted for at all (the second derivatives in momentum equations are omitted) [10] or the turbulent viscosity is set as a constant. This approach gives accurate results for cases without recirculation zones or where turbulence is considered mainly to account for energy losses. An example of a flow with such features is coastal waves. But when turbulent mixing is significant, varying turbulence characteristics should be defined at each point of the flow. To close the SWE, a parabolic eddy viscosity model, a modified mixing length model, and a series of modifications of classical two-parameter difference models such as k-ε [23, 24], k-ω [25] have been developed. The difference in velocity fields obtained for two natural rivers by a depth-averaged solver with standard [26], non-equilibrium [27], and RNG [28] versions of the k-ε model is not significant [29]. Turbulent viscosity defined with non-equation models significantly differs from that obtained with the standard, non-equilibrium, and RNG k-ε models, which give very close values. In [24, 30] it is shown that the difference between standard models and the Chu and Barbarutsi modification [31] becomes obvious only for flows mainly driven by bottom friction. Simple turbulence models such as the parabolic eddy viscosity model and the mixing length model that are widely used with depth-averaged equations [8, 32] due to their simplicity tend to under-predict the turbulent eddy viscosity values for flows where turbulence is mostly two-dimensional, as it is in river flows [11]. The aim of the work presented in this article is improving the mathematical model and the numerical method constructed in [33] for turbulent river flow computations by adding the modified Streeter-Phelps [34] model of self-cleaning mechanism and illustrating the capabilities of the model computing several test-cases of pollutant transport in Tom river. 50 V.V. Churuksaeva, A.V. Starchenko Problem Statement Steady turbulent flow of the viscous incompressible liquid in an open river bed is considered. The river bed has a complex geometry. Assuming that pressure distribution is hydrostatic and water depth is much less than the horizontal size of the area investigated, this flow can be described with 2D SWE: ∂(hu) _^d(hv) = 0 ∂x ду d(hu'-) + d(huv) = -gh д(zb + h) + 1 d(hτxx ) + 1 d(ЬТху ) + (txz )s - (Txz ')b - F ∂x ду ∂x P∂X p∂v P kx" ∂(hUv) + ∂(hV2) =- ∂( Zb + h) + 1 d (h hTyx) + £ d (h Tyy) + (t - (t yz )b - f ∂x ду ду P ∂x P ду P у у where h(X,у) is the water depth; и (x,у),v (x,у) are the depth-averaged horizontal velocities; zb(X, у) is the bed elevation; P is the density; g= 9.81m / s2 is the gravity acceleration; τxx, Tx^ , тух , туу are the depth-averaged components of the viscous stresses and Reynolds stresses tensor; (τxz )s, (τ)s, (τxz )b, (τ)b are corresponding wind stress and bottom friction; and Fx , Fy are the depth-averaged Coriolis force components. The components of the Coriolis force are defined as follows: -4π 4π Fx =----hv sin φ, Fχ =----Ьи sin φ. time у time Here φ is the geographical latitude, time is the length of a day in seconds. In the cases studied, wind stresses (τ x^z )s, (τ yz )s are omitted because their influence is not relevant compared to the effects produced by the slope and bottom friction terms. Turbulence Modelling A Boussinesq hypothesis is used to connect Reynolds stresses with components of the strain velocity tensor by defining effective viscosity (ν+Vt). = P(V+V- )^ау ÷jX I; τχ = 2p(v+v' )∂X ■ 1k ∙ = 2p(v+v' )ду ■ 3- (txz )b = cf I^U; (τyz )b = cf l^v; cf = hn3. Here k is the depth-averaged turbulence kinetic energy; ν is the kinematic viscosity; Vt is the eddy viscosity; Cf is the bottom friction coefficient; and n is the Manning coefficient. In order to account for production, transport, and dissipation of turbulence in the river flow, a depth-averaged high Reynolds (Vt » ν) k - ε model is used in this work. This model was constructed by Rastogi and Rodi [26] from the original k - ε model proposed by Launder and Spalding [35] and was the first differential turbulence model modified for depth-averaged equations. Numerical modelling of pollution transport in Tom River 51 The equations of the model used are: νt = cμ-ε∙ ε ∂x where Ph=νt σε = 1.3; σk = 1.0; c1 = 1.44; c2 = 1.92; cμ = 0.09; D v* v4 c2 !-- Pkv = ck h ’ pεv = cε ,2 ;Ck = I---’cε = 3.6 3'4 у cμ . h h 4cf cf Here k is the depth-averaged turbulence kinetic energy; ε is the energy dissipation; Vt is the eddy viscosity; and n is the Manning coefficient. The constants of the model are assumed to have the same values as in [23] where considerations about their values are also given. Boundary Conditions Inflow boundary conditions are obtained from free surface level and relief data. At the inlet boundary, longitudinal velocity is set to a constant obtained from empirical data for the total discharge of the flow. Because water depth is equal to zero on the river boundary, no wall functions are used and no-slip and no-flow conditions for velocity components are used at wet-dry boundaries. Turbulence kinetic energy and its dissipation fluxes are set to zero at the boundary faces. At the outflow boundary, simple gradient conditions are set for both velocity components and turbulence parameters. Pollutant Transport Modelling Water in natural basins has important self-cleaning mechanisms. Wastewater discharged into rivers is diluted by pure river water and partially precipitates. Organic matter introduced into the river with wastewater oxidises by dissolved oxygen driven by microorganisms and algae, and it decomposes from the sun's radiation. The intensity of selfcleaning processes depends on the water level of a river, flow velocity, and the intensity of mixing, which is mainly defined by turbulence. The model proposed contains the following convection-diffusion equation in order to model transport of pollutants whose velocity is equal to the velocity of the flow and whose concentration is relatively small. ∂hL ∂ ∂ + - , ∂y = -k1Lh -k3Lh. ∂t ∂x Here L (x, y, t) is the depth-averaged concentration of organic matter; k1 is the deoxygenation rate; and k3 is the rate of organic matter sedimentation. 52 V.V. Churuksaeva, A.V. Starchenko The problem considered assumes that organic matter is present in the river and that its concentration and the water temperature change slightly with water depth. The modified dissolved oxygen sag equation [41] is also introduced to the model to define biochemical oxygen demand (BOD) as a criterion of water quality: ∂hD . ∂t сх ∂ + - , Су - = k1 Lh - k2 Dh, where k2 is the reaeration rate and D(x, у, t) is the depth-averaged oxygen deficit, which is defined as the difference between the dissolved oxygen concentration at saturation and the actual dissolved oxygen concentration. In this research, the constants are 0.3 1.0 set to: k1 =-; k2 =-; k3 = 0 [36]. ττ In this approach, oxygen deficit is a criterion of both the intensity of self-cleaning processes and water quality. Both L (X, у, t )and D( x, у, t) fluxes are set to zero at the solid boundaries and are set to constants obtained from empirical data at the inflow boundary. At the outflow boundary, simple gradient conditions are used for both variables. This description of a self-cleaning process is a modification of the classical Streeter-Phelps model. In this model, organic pollution of a river is estimated by BOD and the countervailing influence of atmospheric reaeration. This model could be successfully applied for a timely evaluation of the water quality for a season of the year in several hundred kilometers of a river. Numerical method The solver uses a staggered structured mesh to discretise the spatial domain (Fig. 1). Velocity components are defined in the nodes (rhombs) placed in the midpoints of the edges of the initial mesh. All scalar characteristics of the flow are defined in the nodes (circles) of the initial mesh. Tn ∙f v P' zb , ee Os Fig. 1. Control volume and mesh stencil Discretisation of Convective and Diffusive Terms Convective fluxes in the momentum equations are approximated with a MUSCL scheme [37] that is up to third-order accurate in regions where function is monotonic. Values are interpolated from the centers of the mesh cells to the midpoints of the mesh edges with linear function as follows: Numerical modelling of pollution transport in Tom River 53 ∆X + φe = φp σe , ue >0; ∆X - φe = φe σe , ue ≤0; +δ-δwhere σe^ φ(θe), σe-=-^Φ(θee^') are the slopes that are limited by the ∆X∆X C 2 + θ function φ(θe). φ(θ) is defined as φ(θ) = max 0,minl 2θ^3- numerical scheme that satisfies conditions of Harten's theorem [38]. In this formulation, δ θe =δe(δw ≠ 0); δe = φe -φp. δw An upwind scheme is used to approximate convective terms of the turbulence model equations. - - T.) 'V - φ(θ 2θ.2+^-2) , which gives a Discretisation of the Source Terms ∂(zb+h) ∂(zb+h) Bed slope source terms -gh-------and - gh------- in the momentum ∂X∂y equations and turbulence generation due to gradients of horizontal velocities in the equations of the k-ε because of the construction of the mesh and to obtain second-order spatial accuracy. h ∂(Zb + h) gh ° ≈ l ∂x _ e model are discretised with the central difference scheme ghJhE-^p + , hE and hP > εw > 0(in the river), ∆x 0, hE or hP

Скачать электронную версию публикации