ADAPTIVE FINITE ELEMENT METHOD FOR ANALYSIS OF POLLUTANT DISPERSION IN SHALLOW WATER

Applied Mathematics and MechanicsPublished by Shanghai University,( English Edition) Vol. 26 No.12 Dec.2005Shanghai，ChinaCEditorial Comite of Appl. Math. Mech. ,ISSN 0253-4827Article D: 0253-4827(2005)12-1574-11ADAPTIVE FINITE ELEMENT METHOD FOR ANALYSIS 0FPOLLUTANT DISPERSION IN SHALLOW WATER *Somboon Otarawanna，Pramote Dechaumphai“(Mechanical Engineering Department, Chulalongkor University,Bangkok 10330, Thailand)(Communicated by ZHOU Zhne-wei)Abstract: A fnite element method for analysis of pollutant dispersion in shallow wateris presented. The analysis is divided into two parts:(1) computation of the velocity flowfieid and watr surface elevation, and (2) computation of the pollutant conceatrationfield from the dispersion model. The method was combined with an adaptive meshingtechnique to increase the solution accuracy ，as well as to reduce the computational timeand computer memory. The finite element formulation and the computer programs werevalidated by several examples that have known solutions. In addition, the capability ofthe combined method was demonstrated by analyzing pollutant dispersion in ChaoPhraya River near the gulf of Thailand.key words: shallow water; pollutant disersion; adaptive meshing technique; finiteelement methodChinese Library Classification: 0242. 21; 0368Document code: A2000 Mathematics Subject Classification: 76B15IntroductionNowadays，both the industrial and urban zones in Thailand have increased rapidly. Thedischarge of thermally or chemically polluted water from power stations , industrial plants,and households into rivers has become a threat to water resources. Authorities now requireproof that the environmental impact of a planned discharge will not exceed a certain level,and plant designers must keep the impact below the specified level. For this reason, bothauthorities and plant designers have strong interest in reliable methods for predicting thedistribution of pollutants resulting from a given discharge into a river.The behavior of pollutant dispersion in shallow water is govemed by the conservation ofmass and momentum, and the pollutant transport equation. The analysis may be consideredas a two-dimensional depth-averaged problem by assuming uniform wvlitiee over the depth中国煤化工with their values equal to the depth- averaged velociing differentialTYHCNMHG* Received Jun. 10 ,2004; Revised Aug. 08 ,2005Project supported by the Thailand Research Foundation (TRF)Coresponding author P. Dechaumphai , Professor, Doctor, E mail :fmepdc@ eng. chula. ac. th1574Analysis of Pollutant Dispersion in Shallow Water1575equations are coupled and nonlinear, and thus cannot be solved by analytical methodsespecially for complex flow geometry. Several computational methods have been proposed inthe past. These include the fnite difference method{1-53，the finite volume methodl6.,] andthe finite element method'8 -12. The finite element method is widely used currently becauseit can handle complex geometries effectivelyThe accuracy of solution by the finite element method mainly depends on element sizes.High solution accuracy is obtained if small clustered elements are used in the model.However , the computational time and computer memory are increased if a large number ofelements is used. Adaptive meshing technique5.141 can be applied to increase the analysissolution accuracy, and to reduce the computational time and memory. Such technique placessmall elements in the region of large change in the solution gradients to capture accuratesolution, while locating coarse elements in other regions where the solutions are nearlyuniform.The paper starts by explaining the finite element formulatioo and the correspondingsolution procedure that leads to the development of computer programs. The basic ideabehind the adaptive meshing technique is then described. Finally , the derived finite elementequations and the developed computer programs are validated using simple examples thathave known solutions prior to applying to solve more complex problems.1 Flow Model1.1 Governing equationsThe goveming equations that explain the flow behavior of shallow water flow can bederived by averaging the mass and momentum conservation equations in two-dimensionalover the depth. These equations area(Hu), a( Hu)= 0,(1a)3xdy(u些+o)=((a + 0rx)gvu+D,-u,(1b)( x+ayC2Haσ,g√u+o(uO+o9)= (0hx+(1c)~asadxCH”,where H is the total water depth, u and v are the depth-averaged velocity components in x andy-directions, respectively; g is the gravitational acceleration; and C is the Chezy frictioncoefficient. The stress components σ,σ,,Tzy and Ty are defined byσ,=2vau- g5,(2a), arσy =2vv - g5,(2b)中国煤化工τx=T,=v;YHCNMHG(2c)and whereζ is the elevation of the water surface over the mean surface level as shown inFig. 1, v is the eddy viscosity coefficient.The differential equations, Eqs.(1a) -(1c)，are to be solved with appropriate1576 .Somboon Otarawanna and Pramote Dechaumphaizboundary conditions which are either speci-Datumfied depth-averaged velocity componentsalong edge S,，过xu = u(x,y)，0 = u(x,y)，(3a,b)or surface tractions along edge Sq,之业T。=o,l+Tx m,(4a)T, =txl +σ, m，(4b)Fig.1 The notationa of shallow waterwhere I and m are the direction cosines of theproblemunit vector normal to the boundary edge.1.2 Finite element formulationThe basic unknowns for the sallow water flow problem cresponding to the continuityEq.(1a) and the two momeatum equations (1b) - (1c) are the depth-averaged velocitycomponents u,1 and the water surface elevation 5. The six aode triangular element suggestedin Ref. [15] is used in this study. The element assumes a quadratic interpolation for thevelocity component distributions and linear interpolation for the water surface elevationdistribution according to their highest derivative orders in the diferential Bqs.(1a) - (1c) asu(x,y) =Npup，0(x,y) = Npop，(x,y) = H,5w，(5a,b,c)whereβ = 1,2,.,6;μ = 1,2,3; Ng and H are the element interpolation functions for thevelocity and water surface elevation, respecitvely.To derive the finite element equations , the method of weighted residuals!I7] is applied tothe momentum Eqs.(1b) -(1c) and the continuity Eq.(1a)，[M[(uw, +ow,) (.+.)dA =0,(6a)C2H(6b)[M[(wo. +)-(.+o,)+8 °]d=O, .H[(Hu), +(Hv) ,dA =0,(6c)where A is the element area. Applying Gauss's theorem[16]to Eqs.(6a) -(6c) forgenerating the element boundary integrals, leads to the finite element equations which can bewritten in tensor form asKepyxUp4, +KaprYOqU, -Hwxp +SgprTup +SqrDg +CxUg =Qx，(7a)KsyxUgD, +KqpyrUpDy - Hwy5p +SqpmHp +SqprDg +Cqy"g = Qo,(7b)Jpr(Sn +h)4p +Jwpr(Gp +h)op-Rqs5μ = Rh，(7c)where the coefficients in these equations are defined by中国煤化工Ruw-s_N.M.d. Kw,= IN.NN.HCNMHCA, (8a,b,0)Hw =g,N,HdA，Sepm = 2w{N.,N,dA +wN.,N,dA,(8d,e)Se =2)N,N.,d4，Sapr =v[N,dNa,dA,(8f,g)Analysis of Pollutant Dispersion in Shallow Water1577Smr=unN.N.cd+2N.N,ds，C。 =是0[N&NgdA， (8h,i)C2(ζ +h).Qo ==L[N.T.dS, Qo=HIN.T,ds, Jmr= [H.H.dA， (8j,k,)pJJpar = [H,NH,dA，R. = {H,HVdS.(8m,n)1.3 Computational procedureThe nonlinear finite element equations, as shown in Eqs.(7a) - (7c), are solved byNewton-Raphson iteration method. The method requires writing the unbalanced values in theform,F2 =KqprUgu, +KarUgHy -HwiGp +Sqp= ug+Sqpr Ug +CqUp-Q。，(9a)Fa = KqpyzUgO, +KqprDp", -Hwr5u +Sgqmr4p +SqmrUp +Cp"g - lo,(9b)F, = Jpn(5p +hq)up +Jmpn(5 +hy)vg -Rm(5μ +h,). .(9c)Then application of the method leads to a set of algebraic equations with incrementalunknowns of the form(Gp=)(6*) (Luq)(6x6)> (-Hwx) (6x3)]pr(Aup) (6x1)] p(Fa2)(6x1)](Lqpx)(<6x6) (Gpr)(6x<6) (-Hwp)[OM](x2[ aM])dal ax[N]{H},[{C]I] [wjemsw[G]u[ M)2u}@)[{W|q.dS .(22)or[C]{@} +[[Kv] +[Kc]]{@} = {Q},(23)where the coefficients in this equation are defined by[C] = {W}[N]dA,(24a)[Kv] = {W}[u v][B]dA,(24b)[Kc]{{O}[M]-c[N] {H};{w>{ aM]ix)[ A1])JA+ [{B)I\]- (wumiw[GM](H)[ M)u}.(24c){Q} = [{W|q.dS.(24d)2.3Time discretizationThe explicit recurence relations are applied for time integration of Bq.(23). Theapplication leads toJ[]@}.= (齿[c] +[K,] +[K]J])18]. +10}.(25)which can be solved directly for all nodal values of the pollutant coacentration 0 in the flowdomain.3 ExamplesIn this section, three exarnples are presented. The first two examples are used tovalidate the finite element formulation derived and the computer program developed for flowfield calculation. The validity of pollutant dispersion中国煤化工ast example.3.1 Flow in rectangular canal with variationHCNMHGThe first example presents the analysis of flow behavior in rectangular canal withvariation in bottom profile as illustrated in Fig. 2. Figure 3 shows the finite element modeland the boundary conditions with the maximum inlet velocity of 1.5 m/s, eddy viscosity1580Somboon Otarawanna and Pramote Dechaumphaicoefficientv = 166. 67 m2/s, Chezy coefficientC = 50 m'^/s, and gravitational accelerationg =10 m/s*. This finite element mesh consists of 357 nodes and 160 elements. Figure 4shows close agreement between the computed solution and that of the penalty finite elementmethod in Ref. [9].10 000 m2000m5000m.4000m16000mu≈0,v=0not to scale=0↑l not 10 scale一4mu=0,v= 0Fig.2Problem statement of flow irFig.3Finite element model for flow inrectangular canal with varia-rectangular canal with variation intion in bottom proflebottom profilenot t0 scaleC=50 m!12/s C= 10 m12/sDDDDDDDDD-Computed●Penalty FEM9]Fig.4Comparison of velocities forFig.5Problem staterment of flow in rectangularflow in rectangular canal withcanal with variation in bottom frictionvariation in bottom profle3.2 Flow in rectangular canal with variation in bottom frictionThe geometry and the flow properties of this exanple are identical to those of the firstexample, except for the mean depth of 10 m everywhere and the Chezy coefficient of10 m'"/s in the shaded zone of Fig. 5. Figure 6 shows the boundary conditions and the finiteelement mesh which consists of 187 nodes and 80 elements. The comparison of the computedvelocity distribution with the solution presented in Ref, [9] using penalty finite element isshown in Fig. 7.u=0,v=0.ζ=0.一太DDDBRDDDD中国煤化工inalty FEM[9IMYHCNMHGFig.6kig.7bottom frictionvariation in bottom frictionAnalysis of Pollutant Dispersion in Shallow Water15813.3 Propagation of pollutant through open rectangular channelTo verifty the fnite element formulation and the computer program developed forpollutant dispersion analysis, the concentration evolution of pollutant along an open channelin PFig.8 is studied. The finite element mesh and the initial and boundary conditions areshown in Fig.9. The mesh contains 50 nodes and 72 clements and taking a time interval △t= 0.1s. The transient solution for timest = 6,12,18,24 ,30 and 36s is presented foruniformn velocity U = 0.05 m/s over the domain and dispersion cofficientD = 0.01 m/s.Figure 10 shows the transient solution of pollutant concentration along the x-axis. Thecomputed result is in very close agreement with the exact solution given in Ref. [12].d01n= 0Propagation of pollutantwith U= 0.05 m/s日=1'Initial conditionO= 0 at all points s=>3m一d0/0n=0Fig.8 Problem statement of pollutantFig.9 Finite element model for pollutantpropagation through opea rec-propagation through open rectang-tangular channelular channelFlow in Chao Phraya River1.0.36s30 s24s日0.5| Computed》18sρ !kmExct[12]1●12s0.0y6s0.5Gulf of Thailandx/mFig.10 Comparison of concentration forFig. 11 Computational domain ofpollutant propagation throughCbao Phraya riveropen rectangular cbannel4 Application to Chao Phraya River4.1 Flow in riverThe geometry of the Chao Phraya River is show中国煤化工hows the finiteelement model and the boundary conditions with the:MYHCNMHG of 1.5 m/s,eddy viscosity coefficient v = 15m2/s，Chezy coefficient C = 50m^/s，and thegravitational acelerationg = 9. 81 m/s'. This initial finite element mesh consists of 4111nodes and 1 948 elements.1582Somboon Otarawanna and Pramote DechaumphaiThe numerical solution obtained from the initial mesh is then used to construct thesecond adaptive mesh as described in Section 1. 4. The second adaptive mesh consisting of4 153 nodes and 1 954 elements is shown in Fig. 13. The figure shows smaller elements aregenerated in the regions where large change in velocity gradients occurs. At the same time,larger elements are generated in other regions where the velocity is nearly uniform. With thissecond adaptive mesh, the entire procedure is repeated again to generate the third adaptivemesh with 3411 nodes and 1 594 elements as shown in Fig. 14. The corresponding flowsolution and its detail are shown in Figs. 15 and 16, respectively.Inlet parabolic velocity profile, 11=γ=0V↑1u1=v-0、'= 0ζ=0Fig. 12 lnitial finite element mesh and boundaryFig.13Second finite element mesh forconditions for flow in Chao Phraya riverflow in Chao Phraya iverFig. 14 Third finite element mesh for flowFig. 15 Predicted velocity distribution forin Chao Phraya riverflow in Chao Phraya river4.2 Dispersion in riverContamination due to pollutant discharged from an industrial plant is studied. Fig. 17shows the boundary conditions with the initial condition of no polluant concentrationthroughout the computational domain. The dispersion中国煤化f0 m2/s and thetime interval has the value Ot = 100 s. The final adaTHCNMHGnodel as shownin Fig. 14 is used as the finite element mesh for the dopersou analysis. rigure 18 shows thecomputed concentration contours in the river at three hours after the plant disposal. Detail ofdistribution of pollutant concentration near the plant is also shown in Fig. 19.Analysis of Pollutant Dispersion in Shallow Water15830=00=1Industrial planta日.a=0!km./a0_on=0an0_0Fig.16Detail of predicted velocityFig. 17Computational domain and bounarydistribution for flow in Chaoconditions for pollutant dispersion inPhraya riverChao Phraya river0.06250.1250.Fig. 18 Predicted distribution of pollutant con-Fig.19 Detail of predicted distributioncentration for pollutant dispersion iof pollutant concenfration forChao Phraya river at three hours afterpollutant dispersion in Chaothe plant disposalPhraya river at three hours after5 Concluding RemarksThis paper presents the finite element method for analysis of pollutant dispersion inshallow water. The finite element formulation and its computational procedure is firstdescribed. The corresponding finite element equations area deriverd and the correspondingcomputer programs that can be executed on a sta中国煤化工have beendeveloped. The finite element method is combinedTYHcNMHGgtechniqueinorder to improve the flow solution accuracy. The adaptive meshing technique generates anentirely new mesh based on the solution obtained from a previous mesh. The new meshconsists of clustered elements in regions with large changes in the velocity gradients to1584Somboon Otarawanna and Pramote Dechaumphaiprovide higher solution accuracy. Elsewhere, coarse elements are generated to reduce thecomputational time and computer memory. The results in this paper have demonstrated thecapability of the combined method for the prediction of pollutant dispersion behaviors.References:[ 1 ] Rastogi A K,Rodi W. Predictions of heat and mass transfer in open channels[J].J Hydr EngDiy-ASCE, 1978 ,104(HY3) :397 - 420.[ 2] Vreugdenhill C B ,Wijbenga J H. A computation of flow patterns in rivers[J].J Hydr EngDiv-ASCE , 1982 ,108( HY11) :1296 - 1310.[ 3] Demuren A 0, Rodi W. Calculation of flow and pollutant dispersion in meandering channels[J].J FPluid Mech, 1986 ,172:63 -92.[ 4 ] Molls T, Chaudhry M H. Depth-averaged open channel flow model[J]. J Hydraul Eng-ASCE, 1995 ,121(6) :453 -465. .[ 5 ] Borthwick A G L, Akponasa G A. Reservoir flow prediction by contravariant shallow waterequations[J]. J Hydraul Eng ASCE, 1997 ,123(5):432 -439.[ 6 ] ZhouJ G, Goodwill I M. A finite volume method for steady state 2D shallow water flows[J]. Internat J Numer Method H,1997 ,7(1):4-23.[ 7 ] Yu L,Righetto A M. Depth- averaged turbulence尼w model and applications[ J]. Advances inEngineering Sofware ,2001 ,32(5) :375 -394.[ 8 ] Zienkiewicz 0 C,Heinrich J C. A unified treatment of steady-state shallow water and two-di-mensional Navier-Stokes equations-finite element approach[J]. Comput Method Appl M,1979 ,17/18 :673 - 698.[9 ] Cochet J F ,Dhatt D,Hubct G ,Touzot G. River and estuary flows by a new penalty finite ele~ment[ A]. In:Kawai T (ed). Pinite Element Flow Analyis[ C]. University of Tokyo Press,Tokyo, 1982 ,563 - 570.[10] Leclerc M, Bellemare J, Dumas G, Dbatt G. A fnite element model of estuarian and riverflows with moving boundaries[ J]. Adv Water Resour , 1990,13(4):158 - 168.[11] Sbrestha P L. An integrated model suite for sediment pollutant transport in shallow lakes[J].Advances in Engineering Software , 1996 ,27(3) :201 -212.[12] Tabuenca P, Vila J,Cardona J, Samartin A. Finite element simulation of dispersion in the Bayof Santander J]. Advances in Engineering Sofware ,1997 ,28(5) :313 -332.[13] Peraire J, Vahjdati M ,Morgan K, Zienkiewicz 0 C. Adaptive remeshing for cormpressible flowcomputation[ J].J Comput Phys，1987 ,72(2) :449 - 466.[14 ] Dechaumphai P , Morgan K. Transient thermal-structural analysis using adaptive unstructuredremeshing and mesh movement [ A]. In: Thermal Structures and Materials for High -SpeedFlight[ C]. American Instute of Aeronautics and Astonautics , Washington D C,1992 ,205 -228.[15] Dechaumphai P. Adaptive finite element technique for heat transfer problems{ J]. Joumnal ofEnergy ,Heat and Mass Transfer, 1995 ,17(2) :87 -94.[16] Dechaumphai P. Finite Element Method for Comoutational Fluid Dvnamics [ M]. Chula-longkom University Press , Bangkok ,2002.中国煤化工TYHCNMHG