沙運(yùn)動(dòng)規(guī)律數(shù)學(xué)模型研究及清淤方案優(yōu)選

1、1 基本理論和方法1.1 模型的基本方程 水流模型方程水流連續(xù)方程：(1)水流運(yùn)動(dòng)方程：(2)(3)其中：x、y為空間坐標(biāo)，t為時(shí)間,h為水深，z為水位z=z0+h,z0為床面高程，u,v為流速在x、y方向上的分量，f為柯氏力系數(shù)f=2sin,為地球自轉(zhuǎn)速度，為緯度，Cs為謝才系數(shù)，曼寧公式Cs=h1/6/n,n為糙率，g為重力加速度，為渦動(dòng)粘滯系數(shù)，拉普拉斯算子。 泥沙模型方程    懸沙連續(xù)方程：(4)    河床沖淤演變過(guò)程：(5)其中：C為懸移質(zhì)含沙量，Kx,Ky分別為x、y方向的泥沙擴(kuò)散系數(shù)，P為水流作用下泥沙的起懸量，D為

2、水流作用下泥沙的沉降量，Cbx、Cby分別為x、y方向的河口浮泥層推移輸沙率，z0為床面高程，Cm為床面極限含沙量。1.2 泥沙模型方程中物理量的確定 泥沙沉降量D的確定    泥沙沉降量D與懸移含量C和泥沙沉速*有關(guān)，可表示為：D=1C*(6)其中：1為校正值。 泥沙起懸量P的確定    泥沙起懸量P與飽和挾沙狀態(tài)下的含沙量C*和泥沙沉速*有關(guān)，可表示為：P=2C*(7)    C*通過(guò)LaneKalinske公式確定為：C*=5.55I00.5(u*T/)exp(-(*/u*T)2)1.61(8)式中

3、u*T為潮流與波浪合成摩阻流速，u2*T=(u2+v2)/C2s+fwu20/4(9)u0=H/Tsh(2h/L)，波浪底部最大軌跡流速H波高，T波周期，L波長(zhǎng)，fw波浪阻力系數(shù)。    ab=H/2sh(2h/L)-1，近底波浪質(zhì)點(diǎn)振幅，r*=120d50床面粗度，d50床面泥沙中徑；    ,泥沙的平均沉速，1，2分別為泥沙最大(d1)與最小(d2)顆粒的沉速；I0為床面活動(dòng)泥沙含量。    這樣，將床面層附近泥沙起懸量和沉降量區(qū)分開(kāi)來(lái)處理的方法，具有明確的物理意義，可以清楚地反映出懸移泥沙與床面層內(nèi)

4、運(yùn)動(dòng)泥沙及床面活動(dòng)泥沙發(fā)生相互交換過(guò)程中矛盾的兩個(gè)方面。在波浪掀沙和潮流輸沙過(guò)程中，平衡只是相對(duì)的、暫時(shí)的，而不平衡則是絕對(duì)的、長(zhǎng)期的。泥沙起懸量作為一個(gè)起懸函數(shù)(Pickup Function)的形式出現(xiàn)，有助于反映泥沙在床面層附近的交換機(jī)理，同時(shí)采用Lane-Kalinske公式，其優(yōu)越之處在該公式在建立時(shí)，充分考慮了底部水流湍流脈動(dòng)特性，反映了一定的力學(xué)機(jī)理，從目前流體力學(xué)研究來(lái)看，所不完善的是尚未區(qū)分大尺度擬序結(jié)構(gòu)和小尺度隨機(jī)脈動(dòng)的作用，這方面我們正進(jìn)行深入研究，但作為工程應(yīng)用，該公式已經(jīng)過(guò)了許多實(shí)驗(yàn)室和現(xiàn)場(chǎng)天然實(shí)測(cè)資料的驗(yàn)證，證明是可靠的1,2?？傊@樣選取，有助于逐步弄清泥沙在床

5、面層附近的交換機(jī)理，也有助于從力學(xué)機(jī)理上反映掀沙和輸沙過(guò)程的不平衡特性。1.3 模型方程的數(shù)值求解方法 水流方程的數(shù)值解    (1)雜交有限元方法    雜交有限元方法是三角形單元和等參四邊形單元相交錯(cuò)結(jié)合的有限元方法，這種方法對(duì)于具有復(fù)雜邊界和長(zhǎng)寬尺寸懸殊的水域是非常合適的3,4。    設(shè)單元中的近似函數(shù)?。?10)其中i(x,y)為基函數(shù)，Z、U、V取相同的基函數(shù)；如果單元為三角形，i(x,y)為線性三角元基函數(shù)；如果單元是等參四邊形，則i(x,y)為四結(jié)點(diǎn)等參基函數(shù)。  &#

6、160; 將單元插值近似函數(shù)(10)代入方程(1)、(2)、(3)的Galerkin積分式中，可得雜交有限元方程：(11)(12)(13)    式中A、B1、C1、D1、G、B2、C2、D2、B4、C4分別為系數(shù)矩陣，其元素分別為：Aij=(ij),，B3的元素與B1相同；C3的元素與C2相同。()代表在三角形單元或四邊形單元內(nèi)積分。對(duì)于不同的區(qū)域三角形單元或等參四邊形單元；i,j,k分別為1，2，3或1，2，3，4.    (2)質(zhì)量集中方法    在過(guò)行長(zhǎng)系列時(shí)間過(guò)程量模擬時(shí)，為節(jié)省計(jì)算時(shí)間和計(jì)算機(jī)

7、內(nèi)存，提高計(jì)算效率，按以下三種單元組合情況(圖1)，分別采用集中質(zhì)量法對(duì)單元系數(shù)矩陣進(jìn)行化簡(jiǎn)，具體方法可參見(jiàn)文獻(xiàn)3。(a)三角形單元(b)等參四邊形單元(c)三角形-等參四邊形單元圖1    (3)時(shí)間方向的差分離散和半隱格式    為計(jì)算穩(wěn)定，采用半隱差分格式，先利用Un、Vn和Zn通過(guò)離散式(12)和(13)，計(jì)算出Un+1和Vn+1值：AUn+1=(A-tB1-tG)Un-(C1·t-tfA)Vn-tD1Zn(14)AVn+1=(-tfA-tB2)Un-(C2·t-A+tG)Vn-tD2Zn(15)然后通過(guò)離

8、散式(11)，利用Un+1、Vn+1和Zn，求得Zn+1值：AZn+1=AZn-(tB3+tC3)Zn-(tB4Un+1+tC4Vn+1)(16) 泥沙方程的數(shù)值解 泥沙擴(kuò)散方程和河床演變方程采用集中質(zhì)量迎風(fēng)有限元方法，具體如下：   (1)泥沙濃度場(chǎng)的數(shù)值模擬    為保證含沙度場(chǎng)的計(jì)算精度，同時(shí)考慮漲落潮流的實(shí)際背景，采用集中質(zhì)量迎風(fēng)有限元素方法。下面主要以三角形單元為例加以說(shuō)明，等參四邊形方法與之類似。    取插值函數(shù)為：(17)其中i(x,y)為基函數(shù)；如果單元為三角形，i(x,y)為線性三角元基函數(shù)；如

9、果單元是等參四邊形，則i(x,y)為四結(jié)點(diǎn)等參基函數(shù)。    采用Garlerkin方法，方程(4)可寫(xiě)為：(18)    然后在計(jì)算節(jié)點(diǎn)Pi的相關(guān)單元中尋找其迎風(fēng)單元，并將dcdt在迎風(fēng)單元上離散，設(shè)(x1,y1),(xi,yi),(xi,xi)為節(jié)點(diǎn)Pi迎風(fēng)單元ei的三個(gè)頂點(diǎn)坐標(biāo)，i為迎風(fēng)單元ei的面積。    其中：，這樣，由式(18)可得：(19)    Pik(k=1,2,ni)為Pi的相關(guān)節(jié)點(diǎn)，其中：aiik=aii·iki+aii·iki，這

10、里，E代表Pi節(jié)點(diǎn)的相關(guān)單元點(diǎn)數(shù)。    (2)海床演變的數(shù)值模擬    利用向前差分格式離散床面隨時(shí)間的變化項(xiàng)，然后代入方程(6)，整理得：(20)2 獨(dú)流減河口潮流泥沙數(shù)學(xué)模型的研究2.1模型域的確定及剖分 河口海域模型域：以工農(nóng)兵閘及以遠(yuǎn)10m等深線為東西相界，距離15km；南北以工農(nóng)兵閘為起點(diǎn)，向北3km，向南4km，共7km，105km2的區(qū)域，河口區(qū)域網(wǎng)格剖分考慮清淤及導(dǎo)堤方案的需求，采用嵌套模型，即：在河口5km以內(nèi)，工農(nóng)兵閘附近，采取加密模型網(wǎng)格；該小模型網(wǎng)格為200m×100m,200m×50m等

11、網(wǎng)格，而在河口距工農(nóng)兵閘5km以遠(yuǎn)，采用稀疏的大網(wǎng)格，為1000m×1000m的大模型網(wǎng)格，在大網(wǎng)格與加密小網(wǎng)格之間，采去變網(wǎng)格，使之大小相互匹配，從而使小模型形成統(tǒng)一的整體(參見(jiàn)圖2).圖2 獨(dú)流減河河口區(qū)網(wǎng)格剖分2.2 河口水流數(shù)學(xué)模型的調(diào)試驗(yàn)證 河口模型采用1997年6月13日14日，17日18日，21日22日全潮觀測(cè)資料，以6月13日14日中潮資料作為調(diào)試資料，以大潮17日18日，小潮21日22日作為驗(yàn)證資料，河口地形取1996年10月地形，各測(cè)點(diǎn)平面分布如圖2，共8個(gè)全潮觀測(cè)點(diǎn)。以圖3 水位過(guò)程驗(yàn)證曲線(3號(hào)點(diǎn))圖4 流速模過(guò)程驗(yàn)證曲線(3號(hào)點(diǎn))  &#

12、160; 從以上水位、流速過(guò)程線計(jì)算及實(shí)測(cè)圖可以看出，水位計(jì)算值與實(shí)測(cè)值非常吻合；流速模吻合也較好；同時(shí)漲落潮流速橢圓圖5也基本上符合實(shí)際情況，因此計(jì)算結(jié)果正確，潮流模型是可靠的，可進(jìn)一步用于規(guī)劃工程方案的預(yù)報(bào)。2.3 泥沙模型的調(diào)試與驗(yàn)證 采用非平衡泥沙數(shù)學(xué)模型，結(jié)合97年6月全潮觀測(cè)的泥沙過(guò)程資料，對(duì)模型進(jìn)行調(diào)試與驗(yàn)證。模型驗(yàn)證主要通過(guò)測(cè)點(diǎn)含沙量過(guò)程和床面沖淤變化情況來(lái)確定。    (1)含沙量驗(yàn)證：由于篇幅所限，僅給出一組點(diǎn)的驗(yàn)證結(jié)果(圖5).    (2)床面沖淤變化情況：在含沙量驗(yàn)證正確的基礎(chǔ)上，對(duì)床面演變情況進(jìn)行了驗(yàn)證計(jì)算

13、，采用6月1314日中潮作為典型潮型，并考慮波浪作用，取多年平均波高，以兩個(gè)月為一計(jì)時(shí)段，即每計(jì)算兩個(gè)潮段，一天的沖淤量，使其乘以60，與原來(lái)地形疊加，再進(jìn)行計(jì)算，如此累積，從96年10月計(jì)算到98年3月，其中考慮97年3月挖泥情況，2km內(nèi)泥沙回淤計(jì)算結(jié)果見(jiàn)表1，床面地形變化實(shí)測(cè)與計(jì)算結(jié)果見(jiàn)圖6.2.4 河口泥沙回淤計(jì)算結(jié)果分析 從表1可以看出，96年10月挖泥后，頭兩個(gè)月回淤量17萬(wàn)m3；之后逐月減少；而再次挖泥后，在97年5月至6月之間，回淤量又有所增加，以后又呈逐減趨勢(shì)，直至最后達(dá)到平衡，但泥沙總的回淤量隨時(shí)間的增長(zhǎng)亦仍呈增長(zhǎng)趨勢(shì)，只是增長(zhǎng)幅度有所減緩。表1 回淤量及總回淤量(萬(wàn)m3)

14、時(shí)間回淤量/兩個(gè)月累積回淤量1996.1096.1117.317.31996.1297.116.033.31997.297.314.547.8挖泥清淤-407.81997.597.616.6924.491997.797.812.7937.281997.997.108.3245.501997.1197.125.5851.081998.198.24.4255.42圖5 含沙量過(guò)程驗(yàn)證曲線(3號(hào)點(diǎn))    利用1998年3月與1996年10月地形對(duì)照，2km范圍內(nèi)，泥沙回淤總量為63.9萬(wàn)m3，計(jì)算值(55.42萬(wàn)m3)與實(shí)測(cè)值(63.9萬(wàn)m3)較為接近，泥沙沖淤分布計(jì)

15、算與實(shí)測(cè)分布基本一致，含沙量和床面演變計(jì)算表明泥沙數(shù)學(xué)模型計(jì)算結(jié)果是正確的，所建模型是可靠的，可用于其它問(wèn)題的研究和工程方案的預(yù)報(bào)計(jì)算。圖6(a) 獨(dú)流減河口98年3月實(shí)測(cè)地形(單位：m)圖6(b) 獨(dú)流減河口98年3月計(jì)算地形(單位：m)3 河口清淤方案泄流能力回淤量的計(jì)算及清淤方案的優(yōu)選3.1 清淤方案泄流能力及回淤量計(jì)算 根據(jù)歷年清淤情況，經(jīng)海委規(guī)劃設(shè)計(jì)科研處篩選，提出共12個(gè)清淤方案(見(jiàn)表2)，然后針對(duì)不同的清淤方案，分別進(jìn)行河口泄流能力和泥沙回淤能力計(jì)算，從中篩選推薦清淤方案。利用河口大模型嵌套加密小模型進(jìn)行回淤計(jì)算，以2km內(nèi)各清淤方案的回淤總量作比較，優(yōu)選最佳清淤方案。 

16、;   參照海河口大、中、小潮出現(xiàn)的機(jī)率為19.5%，48.5%，32%，通過(guò)各種潮型的搭配，計(jì)算泥沙回淤，最終得出各清淤方案下河口的泥沙年回淤量。為明顯起見(jiàn)，將各清淤方案下河口的清淤量、年回淤量、清淤的工程效益(年回淤量與清淤量之比值)、河口日均最大泄流量列于表2.表2 清淤方案總表方案部位(m)清寬(m)清深(m)挖方量(萬(wàn)m3)年回淤量(萬(wàn)m3)日均泄流量(m3/s)年回淤量與清淤量之比10+000-0+260200-2.620+260-0+520200-150-2.624038.3365795.8%0+520-0+840150-2.6220+000-0+260200

17、-2.00+260-0+520200-150-2.035.430.9338088.2%0+520-0+840150-2.020+000-0+2602000+260-0+520200-150-3.044.548.23784108%0+520-0+84015040+000-1+1002000+100+0+200200-100-2.531.140.373280130%0+200-1+00010050+000-0+1002000+100-0+200200-100-2.536.243.83450120.9%0+200-1+20010060+000-0+1002000+100-0+200200-100-2

18、.540.3745.93560113.6%0+200-1+40010070+00-0+1002000+100-0+200200-80-2.525.3638.83225152.9%0+200-1+0008080+00-0+1002000+100-0+200200-60-2.519.635.03190178.6%0+200-1+0006090+00-0+1001500+100-0+200150-100-2.416.030.03100187.5%0+200-0+700100100+000-0+260200-2.00+260-0+520200-150-2.041.435.0361584.5%0+520

19、-1+000150-2.0110+000-0+1002000+100-0+200200-100-2.542.3640.5376395.6%0+200-1+500100120+00-0+1002000+100-0+200200-80-2.534.8938.03698109.5%0+200-1+500803.2 清淤方案優(yōu)選    (1)比較相同清淤長(zhǎng)度、相同清淤寬度、不同清淤深度的清淤方案(1、2、3)，可以看出，清淤深度越深(方案3)，河口日均泄量也越大，回淤量也越大。而清淤深度越淺的方案(方案2)，回淤量也越小。在這種清淤方式下，清淤深度越小，清淤量較小，回淤率也較小，相應(yīng)地在滿足泄流要求時(shí)