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

1、1 基本理論和方法1.1 模型的基本方程 水流模型方程水流連續(xù)方程：(1)水流運(yùn)動方程：(2)(3)其中：x、y為空間坐標(biāo)，t為時間,h為水深，z為水位z=z0+h,z0為床面高程，u,v為流速在x、y方向上的分量，f為柯氏力系數(shù)f=2sin,為地球自轉(zhuǎn)速度，為緯度，Cs為謝才系數(shù)，曼寧公式Cs=h1/6/n,n為糙率，g為重力加速度，為渦動粘滯系數(shù)，拉普拉斯算子。 泥沙模型方程    懸沙連續(xù)方程：(4)    河床沖淤演變過程：(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*通過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波長，fw波浪阻力系數(shù)。    ab=H/2sh(2h/L)-1，近底波浪質(zhì)點振幅，r*=120d50床面粗度，d50床面泥沙中徑；    ,泥沙的平均沉速，1，2分別為泥沙最大(d1)與最小(d2)顆粒的沉速；I0為床面活動泥沙含量。    這樣，將床面層附近泥沙起懸量和沉降量區(qū)分開來處理的方法，具有明確的物理意義，可以清楚地反映出懸移泥沙與床面層內(nèi)

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

5、面層附近的交換機(jī)理，也有助于從力學(xué)機(jī)理上反映掀沙和輸沙過程的不平衡特性。1.3 模型方程的數(shù)值求解方法 水流方程的數(shù)值解    (1)雜交有限元方法    雜交有限元方法是三角形單元和等參四邊形單元相交錯結(jié)合的有限元方法，這種方法對于具有復(fù)雜邊界和長寬尺寸懸殊的水域是非常合適的3,4。    設(shè)單元中的近似函數(shù)?。?10)其中i(x,y)為基函數(shù)，Z、U、V取相同的基函數(shù)；如果單元為三角形，i(x,y)為線性三角元基函數(shù)；如果單元是等參四邊形，則i(x,y)為四結(jié)點等參基函數(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)積分。對于不同的區(qū)域三角形單元或等參四邊形單元；i,j,k分別為1，2，3或1，2，3，4.    (2)質(zhì)量集中方法    在過行長系列時間過程量模擬時，為節(jié)省計算時間和計算機(jī)

7、內(nèi)存，提高計算效率，按以下三種單元組合情況(圖1)，分別采用集中質(zhì)量法對單元系數(shù)矩陣進(jìn)行化簡，具體方法可參見文獻(xiàn)3。(a)三角形單元(b)等參四邊形單元(c)三角形-等參四邊形單元圖1    (3)時間方向的差分離散和半隱格式    為計算穩(wěn)定，采用半隱差分格式，先利用Un、Vn和Zn通過離散式(12)和(13)，計算出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)然后通過離

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)泥沙濃度場的數(shù)值模擬    為保證含沙度場的計算精度，同時考慮漲落潮流的實際背景，采用集中質(zhì)量迎風(fēng)有限元素方法。下面主要以三角形單元為例加以說明，等參四邊形方法與之類似。    取插值函數(shù)為：(17)其中i(x,y)為基函數(shù)；如果單元為三角形，i(x,y)為線性三角元基函數(shù)；如

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

10、里，E代表Pi節(jié)點的相關(guān)單元點數(shù)。    (2)海床演變的數(shù)值模擬    利用向前差分格式離散床面隨時間的變化項，然后代入方程(6)，整理得：(20)2 獨流減河口潮流泥沙數(shù)學(xué)模型的研究2.1模型域的確定及剖分 河口海域模型域：以工農(nóng)兵閘及以遠(yuǎn)10m等深線為東西相界，距離15km；南北以工農(nóng)兵閘為起點，向北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)一的整體(參見圖2).圖2 獨流減河河口區(qū)網(wǎng)格剖分2.2 河口水流數(shù)學(xué)模型的調(diào)試驗證 河口模型采用1997年6月13日14日，17日18日，21日22日全潮觀測資料，以6月13日14日中潮資料作為調(diào)試資料，以大潮17日18日，小潮21日22日作為驗證資料，河口地形取1996年10月地形，各測點平面分布如圖2，共8個全潮觀測點。以圖3 水位過程驗證曲線(3號點)圖4 流速模過程驗證曲線(3號點)  &#

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

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

14、時間回淤量/兩個月累積回淤量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 含沙量過程驗證曲線(3號點)    利用1998年3月與1996年10月地形對照，2km范圍內(nèi)，泥沙回淤總量為63.9萬m3，計算值(55.42萬m3)與實測值(63.9萬m3)較為接近，泥沙沖淤分布計

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

16、;   參照海河口大、中、小潮出現(xiàn)的機(jī)率為19.5%，48.5%，32%，通過各種潮型的搭配，計算泥沙回淤，最終得出各清淤方案下河口的泥沙年回淤量。為明顯起見，將各清淤方案下河口的清淤量、年回淤量、清淤的工程效益(年回淤量與清淤量之比值)、河口日均最大泄流量列于表2.表2 清淤方案總表方案部位(m)清寬(m)清深(m)挖方量(萬m3)年回淤量(萬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)比較相同清淤長度、相同清淤寬度、不同清淤深度的清淤方案(1、2、3)，可以看出，清淤深度越深(方案3)，河口日均泄量也越大，回淤量也越大。而清淤深度越淺的方案(方案2)，回淤量也越小。在這種清淤方式下，清淤深度越小，清淤量較小，回淤率也較小，相應(yīng)地在滿足泄流要求時