二維動邊界非恒流數(shù)值模擬在多沙河口海岸的應(yīng)用

二維動邊界非恒流數(shù)值模擬在多沙河口海岸的應(yīng)用

20世紀(jì)70年代以來，隨著科學(xué)技術(shù)的快速發(fā)展，國內(nèi)外許多河口沿海沉積物數(shù)學(xué)模型相繼出現(xiàn)，有效促進(jìn)了河口沿海沉積物研究的發(fā)展。然而，傳統(tǒng)的河口海岸數(shù)學(xué)模型主要是模擬低含沙量的河口海岸的粘土運(yùn)動，而能夠模擬多砂和沿海尺度的河口海岸數(shù)學(xué)模型很少。近年來，一些科學(xué)家在黃河河口沿岸開發(fā)了一套粘土模擬模型。為了充分體現(xiàn)黃河三角洲陸地活動的動態(tài)影響機(jī)制，并在水流、波浪和波浪的作用下建立了黃河河口生態(tài)面的三維砂勢模型，從而取得了良好的效果。該模型是黃河三角洲陸地活動的動態(tài)影響機(jī)、泥沙運(yùn)動與濕地的動態(tài)平衡性。1模型結(jié)構(gòu)1.1潮位及自轉(zhuǎn)角速度式中:U、V分別為潮流速在x及y方向的垂線平均值分量;Z為潮位;Cf為謝才系數(shù);F為柯氏系數(shù),F=2ωsinφ,式中ω為自轉(zhuǎn)角速度,φ為地理緯度;h為水深;Zb為海底起始高程.1.2非均勻沙中泥沙對度計算的影響在河口海岸地區(qū),潮流和波浪是泥沙運(yùn)動的最主要動力.“波浪掀沙和潮流輸沙”使河口海岸地區(qū)的泥沙運(yùn)動極為活躍,但也更為復(fù)雜.因此,在計算和預(yù)報河口海岸地區(qū)的泥沙運(yùn)動和沖淤變化時,僅僅考慮潮流的作用是不全面的,還應(yīng)考慮波浪的作用.二維非恒定流不平衡情況下懸移質(zhì)含沙量隨空間和時間的變化規(guī)律的偏微分方程式(考慮泥沙粒徑不均勻性時,以下各方程式中含沙量和沉速是相應(yīng)于某一粒徑組的)?(hSl)?t+?(hUSl)?x+?(hVSl)?y-εsh(?2Sl?x2+?2Sl?y2)+αωl(Sl-S*l)=0(4)?(hSl)?t+?(hUSl)?x+?(hVSl)?y?εsh(?2Sl?x2+?2Sl?y2)+αωl(Sl?S*l)=0(4)式中:l表示非均勻沙中的第l組泥沙,h為水深;U、V分別x、y方向垂線平均流速;Sl為第l組泥沙垂線平均含沙量;S*l為第l組泥沙水流挾沙力;ωl為第l組泥沙沉速;εs為泥沙擴(kuò)散系數(shù);α為恢復(fù)飽和系數(shù).為方便起見,在以下分析中省略下角標(biāo)l.設(shè)Uc\,Uw和Vc\,Vw分別表示潮流速和波浪質(zhì)點(diǎn)速度在x和y方向上的分量,則在流場任意點(diǎn)x和y方向上瞬時流速為:U=Uc+Uw(5)V=Vc+Vw(6)U=Uc+Uw(5)V=Vc+Vw(6)將上述二式代入式(4)即得到潮流和波浪共同作用下的懸移質(zhì)泥沙不平衡輸沙方程式如下:?(hS)?t+?[hS(Uc+Uw)]?x+?[hS(Vc+Vw)]?y-εsh(?2S?x2+?2S?y2)+αω(S-S*)=0(7)?(hS)?t+?[hS(Uc+Uw)]?x+?[hS(Vc+Vw)]?y?εsh(?2S?x2+?2S?y2)+αω(S?S?)=0(7)波浪對泥沙的作用在一個周期內(nèi)為往返輸移,其單向輸沙作用一般只表現(xiàn)在半個周期內(nèi).由于波浪的周期比較小,一般只有幾秒至十幾秒;而潮流的周期是半天或一天.因此,當(dāng)觀察的時間或計算的步長大于波浪周期時,就可以對上式進(jìn)行波浪周期的平均.由于潮流速與波浪質(zhì)點(diǎn)速度基本無關(guān),而波浪質(zhì)點(diǎn)速度在一個波周期內(nèi)的平均值基本為零,對式(5)在波浪周期內(nèi)取時均后簡化為:?(hS)?t+?(hSUc)?x+?(hSVc)?y-εsh(?2S?x2+?2S?y2)+αω(S-S*)=0(8)?(hS)?t+?(hSUc)?x+?(hSVc)?y?εsh(?2S?x2+?2S?y2)+αω(S?S?)=0(8)式(8)在形式上與式(4)相同,但差別在于式(8)中的挾沙能力S*為潮流和波浪共同作用下的挾沙能力.由于懸移質(zhì)的不平衡輸沙運(yùn)動引起海底地形的沖淤變化,其河床變形方程為:?(hS)?t+?(hSUc)?x+?(hSVc)?y-εsh(?2S?x2+?2S?y2)+γ′?η?t=0(9)式中:γ′為海底淤積泥沙的干容重,η為海底沖淤厚度.通過上述分析可見,在計算河口海岸大范圍泥沙沖淤變化時,只需計算潮流流速場,無需計算波浪流速場,只是在計算挾沙能力時考慮波高因子的影響,這與竇國仁的分析結(jié)果是一致的.1.3關(guān)于窄縫法的邊界在河口海岸的計算區(qū)域內(nèi),由于潮流的漲落,實(shí)際水域邊界是變動的.如何模擬變化的水域一直是比較困難的問題,本文第二作者采用窄縫法解決了這一問題,這里不詳述.窄縫法是設(shè)想在岸灘上各空間步長內(nèi)存在一條很窄的縫隙,縫內(nèi)的水和岸灘前的水相連.這相當(dāng)于把岸灘前的水域延伸到岸灘內(nèi).這樣就可以把計算邊界設(shè)在岸灘的窄縫內(nèi),成為具有一定水深的固定邊界.1.4深水波傳播速度和波長波浪衰減過程以破碎點(diǎn)為界分為兩部分,破碎點(diǎn)以外以底摩擦為主,破碎點(diǎn)以內(nèi)主要是波浪破碎.波浪在破碎前波高可以表示為:Η/Η0=ΚsΚrΚf(10)式中:H為淺水波高;H0為深水波高;Ks\,Kr\,Kf分別為淺水因子、折射因子和底摩擦因子.淺水波傳播速度C和波長L服從如下規(guī)律:根據(jù)微小振幅波理論:式中:α0為波峰線與海岸線的夾角.根據(jù)1984～1994年黃河口的實(shí)測資料得到Kf計算式如下Κf=exp(-0.022Η1/10√Δxˉh√L0)(15)式中:Δx為依據(jù)波所在位置到計算點(diǎn)的距離;ˉh為該距離內(nèi)的平均水深;H1/10、L0分別為依據(jù)波的波高和波長.破波時破波波高Hb與破波水深hb的關(guān)系如下:hb=1.28Ηb(16)破波帶內(nèi)波高衰減過程可由下式計算:Η5/2=0.455h5/2b-1.138β+1h1.5-βb(hβ+1b-hβ+1)(17)式中:β為表征破波帶波能耗散特征的經(jīng)驗(yàn)常數(shù),對淤泥質(zhì)海灘β≈1.6.1.5泥沙沉速與粒徑的關(guān)系在河口地區(qū),細(xì)顆粒泥沙在鹽水中易于發(fā)生絮凝,因此,必須考慮絮凝的作用.如以ωs和ω分別代表絮凝團(tuán)粒及單顆泥沙在水中的平均沉速,兩者之比值F稱為絮凝因子,采用黃建維根據(jù)我國河口淤泥沉降試驗(yàn)建立的關(guān)系式:F=ωsω=7×10-4d-1.950(18)式中:d50為泥沙中值粒徑,以mm計.該式的適用范圍為d50<0.02mm,大于0.02mm的泥沙基本不發(fā)生絮凝現(xiàn)象.對于細(xì)顆粒泥沙在清水中的沉速一般采用Stockes公式計算:ω=118γs-γγgd2ν(19)1.6摩阻流速gcfkj2黃河口地區(qū)的泥沙主要來源于徑流,并在潮流和波浪的共同作用下向深海輸移,因此,單向流(包括徑流和潮流)與波浪共同作用下的挾沙能力公式是否符合實(shí)際,是泥沙數(shù)學(xué)模型能否準(zhǔn)確模擬河口海岸泥沙運(yùn)動規(guī)律和沖淤變化的關(guān)鍵所在.本數(shù)學(xué)模型采用筆者建立的挾沙能力公式:S*=ρsδνΤ+BdSVmu4*ωsu2*k[(1+√gCfk)J1-√gCfkJ2]/(1-ζbζb)Ζ(20)式中:δ為單位床面面積內(nèi)的平均猝發(fā)面積;T+B為近壁區(qū)低流速帶無因次時間間隔,大量的水槽試驗(yàn)表明:T+B≈100,δ≈0.016;d為泥沙粒徑;SVm為單位床面層的極限體積含沙量;ν為水的運(yùn)動粘滯系數(shù);k為卡爾曼常數(shù);b為近底某處距床面距離,一般取ζb=b/h=0.01～0.05;J1=∫1ζb(1-ζζ)Ζdζ;J2=-∫1ζb(1-ζζ)Ζlnζdζ;Z=ωs/kωs.u*為波浪和單向流共同作用下的摩阻流速,采用Bijker的研究成果:u*=u*c[1+12(ζumU)2]12(21)其中,u*c為單向流的摩阻流速;U為垂線平均流速.ζ=akCf√g(22)Bijker給出a=0.45,本文根據(jù)1984～1994年黃河口實(shí)測資料得到a=0.18.um為床面波浪質(zhì)點(diǎn)運(yùn)動最大水平分速:um=πΗΤsinh2πhL(23)u*k為臨界摩阻流速,根據(jù)謝爾茲(Shields)曲線確定.目前,人們對水流分組挾沙能力或者挾沙力級配的機(jī)理尚不十分明了,各家采用的方法也各不相同.本次計算采用韓其為的處理方法,即:式中:pl為懸沙級配.2沙化水的計算如前所述,本數(shù)學(xué)模型中計算公式中的待定參數(shù)是采用黃河口1984～1994年實(shí)測資料和水槽試驗(yàn)資料率定的.下面采用黃河口1996年實(shí)測資料進(jìn)行驗(yàn)證計算.本次計算域?yàn)殂夂恿髀房陂T附近海區(qū)南北長20km(北京坐標(biāo)4180.0～4200.0),東西長22.5km(北京坐標(biāo)20692.5～20715.0).計算網(wǎng)格的空間步長為Δx=250m,Δy=500m,時間步長為Δt=60s.潮流、潮位和含沙量驗(yàn)證計算時間為1996年9月5～6日實(shí)測海流過程,與此相應(yīng)河流入海的流量為1560～1570m3/s,含沙量為21.7～25.0kg/m3;沖淤驗(yàn)證計算時間為1996年6月至11月;開邊界含沙量由式(20)計算求得.計算域內(nèi)各測點(diǎn)和驗(yàn)證斷面的位置如圖1所示.計算潮位和潮流速與實(shí)測值的比較如圖2和圖3所示,表明計算與實(shí)測符合良好.含沙量計算結(jié)果與實(shí)測資料的比較如圖4所示,表明計算與實(shí)測基本符合.根據(jù)1996年汛前6月和汛后10月的施測的兩次地形資料,對數(shù)學(xué)模型進(jìn)行沖淤驗(yàn)證計算.從圖5可見,計算的河口淤積形態(tài)、淤積部位和淤積量與實(shí)測值基本吻合.3根據(jù)河口沉積物流動特征及新濕地發(fā)展預(yù)測3.1泥沙模型及預(yù)測計算的初始地形為1996年10月實(shí)測的汊河流路河口及附近海域的地形,選擇利津站1990～1994年實(shí)測水文系列作為黃河入海的水沙條件,將懸移質(zhì)泥沙分五組,即<0.007mm,0.007～0.01mm,0.01～0.025mm,0.025～0.05mm,>0.05mm.在此基礎(chǔ)上采用泥沙數(shù)學(xué)模型對未來5年濕地的演進(jìn)進(jìn)行了預(yù)報計算.3.2長6周黃河三角洲泥沙數(shù)量的計算進(jìn)入黃河口的泥沙既有小于0.01mm的細(xì)沙又有大于0.025mm的粗沙,究竟多大粒徑的泥沙被輸入到深海,多大粒徑的泥沙淤積在三角洲地區(qū),多大粒徑的泥沙既有一部分輸入深海又有一部分沉積下來,其各自的數(shù)量是多少,這些問題都是黃河三角洲重大和關(guān)鍵的泥沙問題,同時又是迄今為止沒有得到很好解決的問題.此次采用泥沙數(shù)學(xué)模型進(jìn)行了計算,結(jié)果如表1所示.由表1可見,小于0.007mm粒徑組的泥沙有86.6%輸入到深海,0.007～0.01mm粒徑組的泥沙有77.0%輸入到深海,0.01～0.025mm粒徑組的泥沙有10.2%輸入到深海,0.025～0.05mm粒徑組的泥沙有4.1%輸入到深海,大于0.05mm粒徑組的泥沙全部都沉積在三角洲地區(qū).也就是說小于0.01的細(xì)顆粒泥沙大部分(不是全部)輸入到深海,少部分沉積在三角洲地區(qū);大于0.025mm的粗沙絕大部分落淤.3.3河流濕地演變黃河三角洲是世界上新生濕地增長最快的地區(qū),濕地演替與泥沙運(yùn)動密切相關(guān),濕地內(nèi)動植物分布及生態(tài)演替過程與濕地形成的大小、高程、部位的關(guān)系亦息息相關(guān).因此,對新生濕地演變的預(yù)測是十分重要的.本次計算得出了未來5年汊河流路河口海岸濕地大小、高程和部位.圖6給了5年后河口岸線的變化,5年共造陸(即陸地濕地)76.0km2,在原有黃河三角洲的基礎(chǔ)上形成一個新的小三角洲(沙嘴).圖7和表2所示則給出了新生濕地部位、高程和淤積面積,表明大于10m淤積厚度的面積為10.5km2,5～10m淤積厚度的面積為60.4km2,1～5m淤積厚度的面積為167.8km2.4泥沙輸移特性分析針對黃河河口海岸岸線變化劇烈和含沙量變幅大的特點(diǎn),開發(fā)和建立了適合黃河河口海岸應(yīng)用的平面二維動邊界非恒定水流泥沙數(shù)學(xué)模型.模型中采用筆者建立的潮流與波浪共同作用下的懸移質(zhì)挾沙能力公式,適應(yīng)黃河河口海岸含沙量變幅大的特點(diǎn),反映徑流、潮流和波浪等動力因子對泥沙輸移的影響;采用窄縫法解決黃河河口海岸岸線變化劇烈的問題;通過分粒徑組計算以反映海流對不同粒徑組泥沙的輸移特性;河口地區(qū)泥沙顆粒較細(xì),模型考慮了絮凝影響.經(jīng)潮位、潮流、含沙量和泥沙沖淤的驗(yàn)證,計算與實(shí)