High-Resolution Numerical Model for Shallow Water Flows and Pollutant Diffusions

Applicd Mathemaiss and MechanicsPublished by Shanghai University,( English Edition, Vol 23，No 7, Jul 2002)Shanghai，ChinaArticle ID: 0253-4827(2002)07-0741-07HIGH-RESOLUTION NUMERICAL MODEL FOR SHALLOWWATER FLOWS AND POLLUTANT DIFFUSIONSWANG Jia-scng (王嘉松)'，HE You-sheng (何友声)2(1. School of Mechanic Engineering, Shanghai Jiaotong University,Shanghai 200030, P R China;2. Department of Engineering Mechanics , Shanghai Jiaotong University,Shanghai 200030, P R China)(Contributed by HE You-sheng)Abstract: A finite- volume high resolution numerical model for coupling the shallow waterflows and pollutant diffusions was presented based on using a hybrid TVD scheme in spacediscretization and a Runge- Kutta method in time discretization. Numerical simulations formodelling dam break，enlarging open channel flow and pollutant dispersion were implementedand compared with experimental data or other published computations . The validation of thismethod shows that it can not only deal with the problem involving discontinuities and wnsteadyflows，but also solve the general shallow water flows and pollutant diffusions.Key words: shallow water flow; pollutant diffusion; TVD scheme ; finite-volume method;numerical simulationCLC numbers: 035; TV131.4Document code: AIntroductionShallow water flow exists very widely in hydraulic engineering ，Polutant diffusion that ismainly controlled by the shallow water flow has brought increasingly attention in environmentalengineering. The two- dimensional model for describing the flow characteristics has been usedextensively in engineering and the three- dimensional model has been greatly developed. But thetwo- dimensional model would continue to be widely applied in engineering considering the statusin engineering application and quality of computer. However， both the computations of theshallow water flow and the pollutant diffusion are necessary to be further developed. Forexample, the tidal flow and the pollutant dispersion at an estuary with strong upsurge require moreeffective scheme to simulate them than the conventional method. Natural water field with arbitrarycomplex boundaries requires the computational model having a flexibility to treat the complextopology of boundaries . Sewage emissions of large area demand a uniform model to describe theflow so as not to differentiate the computational domain as a near- and a far-filed. In these+ Received date: 2001-06-08; Revised date: 200中国煤化工、Foundation item: the National Natural Science FBiographies: WANG Jia song (1909-), Asscciate ProYHcNMHG@st.cdu.cn);IE You sheng (1931 - ), Professor , Academician of Chinese Academy of Engineering741742WANG Jia-song and HE You- shengsituations it is necessary to develop a high accurate and more flexible computational model. Thetotal variation diminishing (TVD) scheme is used very widely in aerodynamics because of itsnumerical performance of high accuracy， good shock capturing ability and strong stability .Recently the finite- volume method is being applied extensively to the shallow water equations forits advantages of flexibility of geometry， simplicity of a finite-difference method and lesscomputational effort than the finite-element method. The purpose of this paper is using a high-resolution method based on the finite- volume method and the TVD scheme for solving shallowwater flow and environmental hydrodynamics problems. As we know, the water flow impacts onthe pollutant concentration distribution and diffusion, whereas pollutant diffusions do not affect onthe water flows if the density of pollutants is not considered . Therefore coupling or uncoupling tosolve the shallow water equations or the pollutant concentration equations is optional. In thispaper，the two equations were coupled to be solved for convenience to apply the TVD scheme1 Governing EquationsThe governing equations for describing the shallow water flows and pollutant dispersions arethe shallow water equations and the pollutant diffusion equations . The coupling and conservativeform can be written asQ: + [r(Q)], + [G(Q)], = s(Q)，(1)whereQ = (h,q,q,C)",(2)F(Q) = (q。,q验/h + 0.5gh2 ,qxg,/h,qC/h)",(3)G(Q) = (g,,9xg,/h,q;/h + 0.5gh2 ,q,C/h)",4)s(Q) = (S,S2,Sg,S4)",5)Si=0,(6Sz =--(2v。(7)pS3 =-520,影)+(。.(”号; +器)小(8)S4=二(K.)+元(K,)+ Se,(9)where Q,F(Q),G(Q) and S( Q) denote the conserved physical vector, flux vectors in the 必and r- directions, and the source terms, respectively; h is water flow depth (m); g is gravityacceleration, qr，qy，, Thx， Tby，K, and K, are components of discharge per unit width (m2 /s),bottom shear stress and concentration dispersion coefficients， along 第and y- directions ，respectively; C is the pollutant averaged concentration (mg/1); v。= v + v is the effectiveviscosity stress; u is the turbulence viscosity and Sc is the source term .It is assumed that the turbulence viscosity is dominated from the shear stress, then it can beexpressed as!:中国煤化工以=6(10)MHCNMHG.where K is the von Karman constant ( =0.41). u. is the shear speed and given byNumerical Model for Shallow Water Flows743(11)where Tb is the bottom shear stressTb=√坛+坛,，(12)andTby = pghSg(13)The friction Sgx，Sjy can be obtained by the Manning formula_n2qa√q: + q_n2q,√q+qSfx =h 10/3Sp =(14)where n is the Manning coarse coefficient.Based on the experiential formulael11 ，the diffusion coefficients are used asKz = 5.93↓gn1 qx 1/hl/6,(15)K, = 5.93Vgn 1 q, 1/hV6.(16)2 Numerical SchemeConsider any quadrilateral element i, the integral form of Eg.(1) for the domain Q2 ofinterest and the boundary an can be written asedA+jao(F,G). ndl =Jf。_SdA,(17)。J0where A is the area of domain Q，dl denotes the arc length of the boundary日∩,n is an unitoutward vector normal to the boundary a0, which is constituted of four sides . The second term ofabove equation can be given by(F,G). ndl = 2/(Q)t*,(18)where k is the number of the sides, where * is the length of the side h, and fk denotes the outernormal flux vector at the side k .The vector Q is assumed to be a constant over each element. Eq . (17) is further discretizedas follows:孔=-元21由+ s(0.).(19)Let the right hand side terms of above equation be Res;, i.e.,(Q)), = Res;.(20)A two-step Runge- Kutta method is used to discretize (20)，then the second-order accuracyin time can be achieved byQ" = Q{n) + OtRes;(Q!"),Q?= 0! + OtRes;(Q{"),(21)= 0.5(Q{n)中国煤化工The term fh in Eq.(19) can be interpolatedTYHCNMHGTVD scheme .However，as the traditional TVD scheme is a finite-smunC SUINIIN Gilul . unsolved flowvariables are arranged at nodes， satellite elements and their topological relationships are744WANG Jia-song and HE You-shengpresented for convenience to utilize and extend the TVD scheme. All the unsolved variablesare arranged at cells of elements，thus a finite- volume TVD scheme can be obtainedi2). Thedeferential term in the source term is discretized by utilizing a finite-volume form based uponcenter-difference. So the obtained numerical model is of second order accuracy in time andspace. As the TVD scheme is defined and aimed to solve the conservational hyperbolicequations and the source term cannot be ignored in the present model, sometimes may bepositive- negative. Therefore the negative and explicit coefficient terms in the source term aretreated as implicit terms and then removed to the left hand of the descretized equation forkeeping it stable. The detailed approaches can be seen in Refs.[2-4].3 Numerical Examples3.1 Diffraction of dam-break waves around a circular cylinderThe hybrid type of TVD scheme was assessed for modeling a 1-D dam- break wave under thecase of no bed slope and bottom friction. The predicted results of water depth were in exactlyagreement with the theoretical solutions. Being used to the 2-D problem, the scheme was verifiedby comparing the simulations with the experimental data and good agreement was shown. Yangand Hsu-5」 simulated the unsteady diffraction flow around a circular cylinder with the TVDscheme and very satisfied results were obtained. Wave propagation due to dam breaking iscomputed herein with the same geometrical conditions as given by Yang and Hsu. The radius ofthe cylinder is 0.5 m. The dam is located at a distance of 1.2 m behind the cylinder. The waterdepth is set as 10 m at upstream and 2 m at downstream, respectively. A 361 x 101 circular meshis used. In this computation, the concentration equations and no friction are considered， so thesource is not included. Fig.1 show the water surface contours at times t = 0.16s and i = 0.2safter the dam failure. It can be seen that a bow-type shock takes place with the bore- solidreflection behind the cylinder, and a complex structure is formed with the bore- bore interaction ahead of the cylinder .2.01.5-1.s1.00.5-0.- 0.5--0.5一1.0-- 1.0-1.5--2.0+- 2.022一-20x/m中国煤化工(a)↑=0.16sFig.1 Constant of water depth for dam-brMYHCNMHGoreakNumerical Model for Shallow W ater Flows7453.2 Enlarging open channel flowA recirculation flow exists widely in civil and hydraulic engineering， such as in theneighborhood of seawall, port and gulf， around a slice or a bydraulic structure. There is alsorecirculation flow when pollutants are discharged into river or sea. The recirculation characteristicsof enlarging open channel flow are always used to validate turbulence models as a typical testexample for flow of cross-section variation .Generally the reattachment length of the recirculating length of op is about seven times of thebackward step H from experiments while the standard k-ε model calculates a value of5. 6H°J . Inthis example the geometry and mesh are given as follows : the backward step width of the channeIH = 0.5m，the width of the inlet section is3.0H, the length of the channel is 28.0H.4 000grids are used， The computed velocity field is1.5shown in Fig.2. The recirculating length is about是1.05.6H ~ 6.0H which is close to the result with the0.5Fstandard k-ε turbulent model. It is shown that thepresent model has a considerable accuracy .Moreover, since it is not necessary to solve the kand ε equations and determine the solid boundaryFig.2 Velocity of enlargingconditions with the wall function, the taken CPUopen channeItime is far less than that with the k-ε models .3.3 Side discharge into open channelThe sewage is discharged into rivers or sea generally with the type of water jets . Accordingto the locations of waste disposal, there are usually two types to waste disposal: side dischargeand bottom discharge. When the waste is ejected into water bodies ，the environmental waterbodies will be wrapped and taken in by the fluid upstream or downstream around discharge portand formed a mixture center. The environmental fluid will be removed by the jets and the jetsthemselves will also be curved by the environmental fluid . Since the interaction between the jetsand environmental fluid are confined by the solid boundary of riverbank, a recirculation zone willbe formed surrounding the bank and the discharge port. The size and the variation of the zone areaffected by the velocity ratio of Vo/ Uo (where Vo is the jet velocity, Uo is the mean velocity ofthe main channel)， shape of bank，flow characteristics of environmental fluid, etc. It cannot besolved with the commonly theoretical analysis method , but always with the numerical simulation .The present study considers a following case: the width and the mean velocity of the mainchannel with smooth bed are B = 1.82 mand Uo = 0.1 m/s, respectively; The side discharge isat a right- angle to the channel with the width and jet velocity ofb = 0.225 mand V = 4.0Uo，respectively; The water depth is h = 0 .06. The concentration of pollutant at discharge port isCo = 1.0 mg/l. .The boundary conditions are used as follows: 1) at unstream.”= h197, =0, c=0(inlet of main channel); qx =0, q, = hVo, C =中国煤化工at downstream :0φ/I + V_dφ/On = 0(whereφ = qx,qy,h,C;n ilFY'RCN M H Gmal to the ouletboundary; Vn is the velocity component of the outward direction normal to the outlet boundary);746WANG Jia-song and HE You-sheng3) at solids: the tangential and normal fluxes are given as zero .Select the length of the channel is 1 .08m at upstream of the discharge port and 2. 52m atdownstream. Bottom slope is ignored. With a 90 x 40 nonuniformly distributed grid， theproposed model was applied to predicted the flow patterns. The results of velocity vector fieldsand concentration distribution when the flow is steady ( after 44. 87s} under a specified totalresidual are shown in Fig .3.20F 0.062 0.125 0.187 0.250 0.312 0.375 0.437.1.0.8)10，11，121314150.500 0.562 0.625 0.687 0.750 0.812 0.875 0.9370.661.0F0.40.50.21.52.53%/mx/mFig.3 Velocity field and concentration distribution for pollutant diffusion of side discharge caseThe size of the main recirculation zone behind the discharge port is characterized by thelength L and lateral extent H of recirculating eddy. It was found experimentally that the size andshape of the recirculating eddy could be described by a single relation, i.e., the momentum fluxratioM = R2b/B (R = Vg/U)l9). In the present case M is about 0.2. Ref.[7] gave theexperimental results and Ref.[8 ] presented numerical results by using the depth-averaged k-Eturbulent model, Power-law scheme, and the approximately Riemann solvers with SIMPLECmethod. A comparison of these results is given in Table 1.Table 1 Comparison of characteristics of recirculation zoneexperimental2-D model3-D model9]parameterdata(7]PresentRoe'Power-lawl81 QUICKHybnidL/B1~1.70.91.04).711.088 .0.56H/B≈0.20. 1420.1850.1320.23.1It can be seen that the results are different with various models and numerical methods . Thepredictive ability using the proposed model and numerical method is better than solving the 2-Dk-E turbulent model with the Power-law method and solving the 3-D k-ε turbulent model with theHybrid method, but not as good as solving the 2-D kε turbulent model with the Roe’s methodand solving the 3-D k-ε turbulent model with the QUICK scheme. Thus both the numericalmethod and the mathematical model are very important. Since the present model is actually a“zero-equation” turbulence model that is much simpler and the computational efficiency is betterthan the general h~ε turbulence model, the results are hish accnrate and reasnnabie .中国煤化工4 ConclusionCNMHGA hybrid type TVD scheme and finite-volume method was introduced to solve the shallowNumerical Model for Shallow Water Flows747water equations and pollutant diffusion equations in this paper. A high-resolution finite-volumenumerical model for coupling the two equations was developed. Taking several typical flows forshallow water and water environment problems as examples, this paper modeled dam-break wavediffraction around a circular cylinder, enlarging open channel flow and pollutant dispersion forwaste side discharge case. The results were compared with experimental data and other numericalsimulations and good agreement was achieved . This work shows that the presented model can notonly solve the problem involving discontinuities and unsteady flows， but also well computegenera! shallow water flows and pollutant diffusions. It can provide a high accurate, good stableand universally applying computational method. The further work is to extend it to practicalestuary or gulf for more broad engineering problems .References :[1] Elder J W. The dispersion of marked fluid in turbulent shear flow[J]. J Fluid Mech, 1959, (5):546 - 560.2] WANG Jia-song, NI Han-gen. A high-resolution finite- volume method for shallow water equations{J]. J Hydrodyramics, Ser B ,2000,(1): 35-41.[3] WANG Jia-song, NI Han-gen, HE You- sheng. Finite diference TVD scheme for computation ofdam-break problems[J]. ASCE J Hydr Eng ,2000, 126(4) :253 - 262.[4] WANG Jia-song . Finite-difference and finite volume high-resolution methods for shallow water flowand pollutant diffusion[ R] . Post Doc Report of Shanghai Jiaotong University ,2000 ,8. (in Chinese)5] YangJ Y, Hsu C A. Computations of frce surface flows, part 2: 2D unsteady bore diffraction[J] .J Hydr Res , 1993 ,31(3) :403 -412.[6] CHEN Jing .ren. Turbulence Models and Finite- Analysis Method[ M]. Shanghai: Shanghai JiaotongUniversity Press, 1989. ( in Chinese)[7] McGuirkJ J, Rodi W. A depth- averaged mathematical model for the near field of side dischargesinto open channei flow[J]. J Fluid Mech , 1978 ,86(4):761- 781.[8] Yc J. McCorquodale J A. Depth- averaged hydrodynamics model in curvilinear collocated grid[J].ASCE J Hydr Eng , 1997,123(3) :380 - 388.9] Demuren A O, Rodi W. Side discharge into open channels: Mathematical model[J]. ASCE J HydrEng , 1982, 109( 12):1707 - 1722.中国煤化工MYHCNMHG