基于垂向坐標(biāo)變換的三維潮流輸沙數(shù)學(xué)模型

基于垂向坐標(biāo)變換的三維潮流輸沙數(shù)學(xué)模型

甘寧洋的開發(fā)和建設(shè)的前提和基礎(chǔ)，需要對甘寧洋的水動(dòng)力和沉積物特征有詳細(xì)的了解。伶仃洋的基本特征是海區(qū)范圍大,口門多(其西北部有四條河流匯入),島嶼多,地形不規(guī)則,流態(tài)及泥沙運(yùn)動(dòng)復(fù)雜,上下游位相不等等。因此,進(jìn)行伶仃洋海區(qū)三維潮流及懸沙場的數(shù)值模擬是完全必要的。文獻(xiàn)建立了伶仃洋三維潮流數(shù)學(xué)模型并進(jìn)行了三維潮流場的數(shù)值模擬,本文則在文獻(xiàn)的基礎(chǔ)上進(jìn)一步建立三維懸沙運(yùn)動(dòng)的數(shù)學(xué)模型,對伶仃洋的懸沙場進(jìn)行數(shù)值模擬,目的是了解其平面分布特征和空間結(jié)構(gòu)特征,為伶仃洋的港口航道工程的規(guī)劃和施工提供依據(jù)。1數(shù)學(xué)模型的構(gòu)建1.1基本方程和解決方案1.1.1懸沙對流擴(kuò)散方程三維潮流輸沙的基本方程包括潮流連續(xù)方程、動(dòng)量方程和懸沙對流擴(kuò)散方程。在f平面右手直角坐標(biāo)系(x,y,z)中,方程形式為潮流連續(xù)方程:?u?x+?v?y+?w?z=0(1)潮流動(dòng)量方程:?u?t+u?u?x+v?u?y+w?u?z=-1ρ?Ρ?x+??x(Νx?u?x)+??y(Νy?u?y)+??z(Νz?u?z)+fv(2)?v?t+u?v?x+v?v?y+w?v?z=-1ρ?Ρ?y+??x(Νx?v?x)+??y(Νy?v?y)+??z(Νz?v?z)-fu(3)?Ρ?z=-ρg(4)懸沙對流擴(kuò)散方程:?S?t+u?S?x+v?S?y+w?S?z=??x(Dx?S?x)+??y(Dy?S?y)+??z(Dz?S?z)+?(ωS)?z(5)式(1)～(5)中:x,y,z為直角坐標(biāo)系坐標(biāo)(z軸垂直向上,原點(diǎn)置于靜止海面);u,v,w分別為沿x,y,z方向的流速分量;P為靜水壓力;ρ為水體密度;g為重力加速度;ω為泥沙沉降速度;S表示水體含沙濃度;Nx、Ny、Nz分別為x、y、z方向的水體渦動(dòng)粘性系數(shù);Dx、Dy、Dz分別為x、y、z方向的懸沙擴(kuò)散系數(shù);t為時(shí)間;f為柯氏參數(shù);式(4)為靜水壓力方程,認(rèn)為水體在垂向滿足靜水壓力方程的假定或說近似在海岸河口地區(qū)被廣泛采用。1.1.2懸沙對流擴(kuò)散方程引入σ坐標(biāo)變換,即σ=z-ζΗ,經(jīng)變換后,(x,y,σ)坐標(biāo)系中的方程形式為:潮流連續(xù)方程:?ζ?t+?Ηˉu?x+?Ηˉv?y=0(6)潮流動(dòng)量方程:?uΗ?t+?u2Η?x+?uvΗ?y+?uW?σ=-gΗ?ζ?x+Η[??x(Νx?u?x)+??y(Νy?u?y)]+1Η??σ(Νz?u?σ)+fΗv(7)?vΗ?t+?uvΗ?x+?v2Η?y+?vW?σ=-gΗ?ζ?y+Η[??x(Νx?v?x)+??y(Νy?v?y)]+1Η??σ(Νz?v?σ)-fΗu(8)W=??x[Η(ˉu(σ+1)-∫σ-1udσ)]+??y[Η(ˉv(σ+1)-∫σ-1vdσ)](9)w=W+u(σ?Η?x+?ζ?x)+v(σ?Η?y+?ζ?y)+(σ?Η?t+?ζ?t)(10)懸沙對流擴(kuò)散方程:?S?t+u?S?x+v?S?y+WΗ?S?σ=??x(Dx?S?x)+??y(Dy?S?y)+1Η2??σ(Dz?S?σ)+1Η?ωS?σ(11)式中:ζ為相對于x,y坐標(biāo)平面的水面起伏;H為實(shí)際水深,H=h+ζ,h為海底到x,y平面的距離,即靜水深。由σ表達(dá)式可知,在海面處,z=ζ,σ=0;在海底處,z=-h,σ=-1,因此在(x,y,σ)坐標(biāo)系下,海面和海底都成了坐標(biāo)平面,σ介于-1～0之間,為新的垂向坐標(biāo),W為(x,y,σ)坐標(biāo)系下的垂向速度分量,ˉu、ˉv分別為x、y方向速度分量沿水深的平均值,即ˉu=∫0-1udσ,ˉv=∫0-1vdσ。1.1.3確定條件x、y、坐標(biāo)系中的確定條件如下2水面邊界條件ρΝzΗ(?u?σ,?v?σ)=(τsx,τsy),W=0?式中:τsx、τsy分別為水面風(fēng)應(yīng)力在x、y方向的分量。泥沙沉降臨界切應(yīng)力ρΝzΗ(?u?σ,?v?σ)=(τbx,τby),W=0-DzΗ?S?σ-ωS={E(τbτe-1)τb＞τe0τd≤τb≤τeωS(τbτd-1)τb＜τd}(13)式中:τd、τe、E分別為泥沙沉降臨界切應(yīng)力、沖刷臨界切應(yīng)力和侵蝕常數(shù);τbx和τby分別為床面摩擦力矢量τb在x和y向的分量:(τbx,τby)=ρgub2+vb2c2(ub,vb)(14)式中:c為Chezy系數(shù),c=1nΗ16?n為Manning糙率系數(shù);ub、vb為u、v在近水底處的值。海岸邊界流速沿岸邊界法線方向分量為零;S沿外法線方向?qū)?shù)為零。流水方向單位變量s0水位或流速分量u,v或uˉ,vˉ給定;S在流出計(jì)算域時(shí)有:?S?t+(V?n)?S?n=0(15)當(dāng)流入計(jì)算域時(shí)S=S0(16)式中:S0為已知值,V為水流流速矢量,n為水邊界外法向方向單位矢量。2.2模型建立的基本方程采用簡單的有限差分方法對式(6)～(11)進(jìn)行離散求解,在垂向上采用等間距分層,水平方向上采用矩形網(wǎng)格,將所有變量均布置在同一網(wǎng)格節(jié)點(diǎn)(i,j,k)上(i,j表示水平方向節(jié)點(diǎn),k表示第k層),時(shí)間導(dǎo)數(shù)采用前差,水平方向空間導(dǎo)數(shù)采用顯式中心差分(對于式(11)中的水平對流項(xiàng)的空間偏導(dǎo)數(shù)采用“迎風(fēng)”格式),σ向空間導(dǎo)數(shù)采用隱式中心差分,經(jīng)整理后,由式(6)的離散結(jié)果將得到求解ζ的顯式差分方程,由(7)、(8)、(11)的離散結(jié)果將分別得到沿垂向求解u、v、S的三對角方程,由(9)、(10)的離散結(jié)果可分別得到求解W和w的顯式差分方程。設(shè)Δt為時(shí)間步長,在KΔt→(K+1)Δt(K=0,1,2,3,4,……,)內(nèi),關(guān)于S的差分方程為:Ai,j,kΚ+1Si,j,k-1Κ+1+Bi,j,kΚ+1Si,j,kΚ+1+Ci,j,kΚ+1Si,j,k-1Κ+1=Di,j,kΚ(17)式中:Ai,j,kΚ+1=-Wi,j,kΚ+12Ηi,jΚ+1Δσ-(Dz)i,j,kΚ+1(Ηi,jΚ+1Δσ)2+(Dz)i,j,k+1Κ+1-(Dz)i,j,k-1Κ+14(Ηi,jΚ+1Δσ)2+ω2Ηi,jΚ+1Δσ(18)Bi,j,kΚ+1=1Δt+2(Dz)i,j,kΚ+1(Ηi,jΚ+1Δσ)2(19)Ci,j,kΚ+1=Wi,j,kΚ+12Ηi,jΚ+1Δσ-(Dz)i,j,kΚ+1(Ηi,jΚ+1Δσ)2-(Dz)i,j,k+1Κ+1-(Dz)i,j,k-1Κ+14(Ηi,jΚ+1Δσ)2-ω2Ηi,jΚ+1Δσ(20)Di,j,kΚ+1=Si,j,kΚΔt-(u?S?x)i,j,kΚ-(v?S?y)i,j,kΚ+[??x(Dx?S?x)]i,j,kΚ+[??y(Dy?S?y)]i,j,kΚ(21)式(21)中的水平擴(kuò)散項(xiàng)空間偏導(dǎo)數(shù)采用中心差分離散;為了增加模型的穩(wěn)定性和減少數(shù)值耗散,對式(21)中的水平對流項(xiàng)空間偏導(dǎo)數(shù)采用三階迎風(fēng)格式離散。以?S?x為例,離散格式如下:對于內(nèi)點(diǎn)(i=3,M-2,M為x方向網(wǎng)格點(diǎn)數(shù)):當(dāng)ui,j,k≥0時(shí):?S?x=(Si-2,j,k-6Si-1,j,k+3Si,j,k+2Si+1,j,k)/6Δx+o(Δx)3(22)當(dāng)ui,j,k≤0時(shí):?S?x=(-2Si-1,j,k-3Si,j,k+6Si+1,j,k-Si+2,j,k)/6Δx+o(Δx)3(23)對于左邊界點(diǎn)(i=1,2):?S?x=(-11Si,j,k+18Si+1,j,k-9Si+2,j,k+2Si+3,j,k)/6Δx+o(Δx)3(24)對于右邊界點(diǎn)(i=M-1,M)?S?x=(-2Si-3,j,k+9Si-2,j,k-18Si-1,j,k+11Si,j,k)/6Δx+o(Δx)3(25)式中:Δx為x方向的空間步長。2數(shù)值模擬與結(jié)果分析2.1潮過程及懸沙運(yùn)動(dòng)參數(shù)數(shù)值模擬的區(qū)域?yàn)閺幕㈤T至珠?！箦u斷面,南北距離大約70km,東西距離大約51km,模擬海區(qū)有6處開邊界,即虎門、蕉門、橫門、洪奇瀝、香港水道和珠海—大濠島斷面(見圖1),水平方向用矩形網(wǎng)格剖分計(jì)算域,空間步長Δx=Δy=320m,計(jì)算網(wǎng)格節(jié)點(diǎn)17903個(gè),垂向分11層,Δσ=0.1,泥沙計(jì)算時(shí)間步長取200s。數(shù)值模擬的資料為1991年12月19日9時(shí)至21日13日共53小時(shí)的潮過程(枯季)和1992年7月15日17時(shí)至17日21時(shí)共53小時(shí)的潮過程(洪季)的懸沙運(yùn)動(dòng)。模型中的各種參數(shù)的確定方法如下:2.1.1Dx、Dy、Dz泥沙水平擴(kuò)散系數(shù)Dx、Dy用歐拉公式確定:{Dx=5.93gΗ|u|/cDy=5.93gΗ|v|/c(26)泥沙垂向擴(kuò)散系數(shù)Dz采用修正的Munk型公式確定,在(x,y,σ)坐標(biāo)系中的表達(dá)式為:Dz={a1-kmu*σ(1+σ)Η1+σ＜0.5a1+kmu*Η/41+σ≥0.5(27)式中:u*為摩阻流速,u*=τb/ρ;km為卡門常數(shù)(km=0.4),a1為修正Munk型公式所加的一個(gè)小量,以使z=-h(即σ=-1)時(shí)Nz不為零。2.1.2τd、τe、E沉降臨界切應(yīng)力τd=0.08Nm-2,沖刷臨界切應(yīng)力τe=0.45Nm-2,侵蝕常數(shù)E=5.5×10-4kgm-2s-1。2.1.3ω本文將泥沙沉降速度取為常值。海水中細(xì)顆粒泥沙的沉降速度決定于其絮凝當(dāng)量的大小。試驗(yàn)表明,這樣的當(dāng)量粒徑一般約為0.015～0.03mm,其相應(yīng)沉降速度為0.01～0.06cm/s,本文取ω=0.04cm/s。2.2垂向分布及含沙濃度過程線驗(yàn)證在所模擬的枯、洪季資料中,在模擬海區(qū)內(nèi)均有10個(gè)水位測站和22個(gè)流速、流向和含沙濃度的逐時(shí)測量資料(見圖1),域內(nèi)有6個(gè)水位站、14個(gè)流速、流向測站和12個(gè)含沙濃度測站可用于模擬驗(yàn)證(其余各站均在模型邊界處,用作邊界條件)。關(guān)于潮流的模擬驗(yàn)證情況(水位過程線,垂線平均流速、流向過程線,流速垂向分布等)詳見文獻(xiàn)。限于篇幅,這里給出6個(gè)測站的垂線平均含沙濃度過程線驗(yàn)證,見圖2(橫軸表示時(shí)間h,縱軸表示流速m/s);以F4測站為代表的含沙濃度垂向分布驗(yàn)證見圖3(橫軸單位為kg/m3,縱軸為σ坐標(biāo),驗(yàn)證了1991年12月19日10時(shí)至20日10時(shí)共25個(gè)小時(shí)的資料),圖2和圖3中,實(shí)線表示計(jì)算結(jié)果,黑方點(diǎn)表示實(shí)測。由圖可見,計(jì)算值與實(shí)測值吻合較好;圖4、5、6分別給出了模擬計(jì)算的1992年7月17日17時(shí)的不同水層的懸沙平面場,圖7給出了模擬計(jì)算的對應(yīng)平面場相應(yīng)時(shí)刻的沿航道的懸沙剖面分布,圖中反映的懸沙分布特征(如西灘存在高含沙濃度區(qū),西灘高含沙濃度水體不進(jìn)入伶仃航道等)與衛(wèi)星遙感分析結(jié)果是一致的。2.3中灘和東灘含沙濃度的日變化伶仃洋的泥沙來源除了河流徑流輸沙外,灘地泥沙的再懸浮也是一個(gè)主要因素。伶仃洋含沙濃度分布有如下特點(diǎn):(1)伶仃洋的懸沙運(yùn)動(dòng)隨潮流呈周期性變化;(2)灘地處的含沙濃度大于深槽處含沙濃度(這是灘地泥沙再懸浮的必然結(jié)果);(3)伶仃洋的高含沙濃度區(qū)在西灘,中灘和東灘在外動(dòng)力比較強(qiáng)的情況下,含沙濃度也較高,但與西灘含沙濃度相比,泥沙濃度量值要小;(4)伶仃洋的含沙濃度峰值主要是由大落潮引起的,但含沙濃度的峰值與水流速的峰值并不同步,而是有2～3小時(shí)的遲后(與文獻(xiàn)結(jié)論一致);(5)洪季含沙濃度高于枯季,西部含沙濃度高于東部;(6)伶仃洋西灘高含沙濃度區(qū)與伶仃洋航道深槽之間存在一含沙濃度突變線(渾水峰面線),隨著漲落潮,峰面線呈不規(guī)則的蛇形分布,但峰線位置基本上在西灘上擺布,一般情況下,不進(jìn)入航道(圖4、5、6、7),這與現(xiàn)場情況和衛(wèi)星遙感分析結(jié)果是一致的。(7)無論是枯季還