泥沙參數(shù)選擇方法及在懸沙濃度研究中的應(yīng)用

董 佳，張寧川

（大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連116024）

目前，在港口和海岸工程工可研階段，論證工程對(duì)泥沙運(yùn)動(dòng)及海床演變影響多采用數(shù)值模擬方法［1-3］。該方法中需要確定諸多影響泥沙運(yùn)移的基本參數(shù)，如泥沙沉降速度、泥沙擴(kuò)散系數(shù)、控制侵蝕速度的比例值、軟底床的侵蝕系數(shù)、硬底床的侵蝕指數(shù)、床面糙率、河床淤積物的干密度、床層間的泥沙轉(zhuǎn)換率等。上述參數(shù)選擇恰當(dāng)與否將直接影響數(shù)值計(jì)算的可靠性。通常情況，研究者通過實(shí)測(cè)資料驗(yàn)證或經(jīng)驗(yàn)獲取得，然而當(dāng)實(shí)測(cè)資料缺乏或經(jīng)驗(yàn)不足時(shí)，往往難以得到符合實(shí)際的計(jì)算結(jié)果。

對(duì)于泥沙沉降速度可通過物理模型試驗(yàn)或較為成熟的經(jīng)驗(yàn)公式計(jì)算得到較為可靠的數(shù)值；對(duì)于以潮流為主要環(huán)境動(dòng)力的淤泥質(zhì)海床而言，經(jīng)過數(shù)值計(jì)算表明，泥沙擴(kuò)散系數(shù)，控制侵蝕速度的比例值，軟底床的侵蝕系數(shù)，硬底床的侵蝕指數(shù)等對(duì)懸沙濃度的影響較小。因而以潮流為主要?jiǎng)恿Φ挠倌噘|(zhì)海域懸沙濃度計(jì)算中，床面糙率、河床淤積物的干密度、床層間的泥沙轉(zhuǎn)換率3個(gè)參數(shù)的最佳取值需要進(jìn)行詳細(xì)探討。

首先建立二維水動(dòng)力模型，在通過實(shí)測(cè)資料驗(yàn)證的基礎(chǔ)上，進(jìn)行懸沙濃度計(jì)算，定量地討論了每個(gè)參數(shù)的選擇對(duì)懸沙濃度的影響程度。在床面糙率、河床淤積物的干密度、床層間的泥沙轉(zhuǎn)換率可變范圍內(nèi)，分別調(diào)整各個(gè)參數(shù)的取值，以計(jì)算得到的懸沙濃度變化與實(shí)測(cè)資料吻合為目標(biāo)，搜尋得到上述3個(gè)參數(shù)的最優(yōu)取值。在此基礎(chǔ)上，對(duì)洋山港全水域中各特征點(diǎn)的懸沙濃度進(jìn)行了模擬并和實(shí)測(cè)值吻合較好，驗(yàn)證了所選擇參數(shù)的可靠性。

1 水流泥沙數(shù)學(xué)模型

1.1 控制方程

水流泥沙數(shù)學(xué)模型控制方程包括潮流控制方程和懸沙濃度控制方程。

對(duì)于潮流控制方程，采用包含二階紊動(dòng)項(xiàng)的平面二維非恒定流體運(yùn)動(dòng)方程

式中：S為懸沙濃度；Sm為泥沙的沉降和再懸浮濃度；Cb為底層泥沙濃度；τb為瞬時(shí)床底剪切應(yīng)力；τcd為淤積臨界剪切應(yīng)力；τce為侵蝕臨界剪切應(yīng)力；ε為泥沙水平擴(kuò)散系數(shù)；α為軟底床的侵蝕系數(shù)；n為硬底床的侵蝕指數(shù)；E為控制侵蝕速度的比例因子。

1.2 定解條件

（1）初始條件：在計(jì)算的初始時(shí)刻給出潮位、懸沙濃度等初始值，其選取范圍較為寬松，因?yàn)槌跏紬l件的誤差可以在正確的邊界條件控制下很快消失。

（2）潮流模型開邊界條件：在海域開邊界，海面水位的邊界條件由邊界處主要分潮的調(diào)和常數(shù)計(jì)算可得

式中：Em為該點(diǎn)相對(duì)于平均海平面的水位；ai、wi、φi分別是第i格分潮的振幅、頻率和遲角。

（3）潮流模型固體側(cè)邊界條件：在研究區(qū)域的固體側(cè)邊界處，一般假定垂直于固體海岸的法向速度為零，即。

（4）懸沙濃度計(jì)算模型開邊界條件：在入流時(shí)段S為給定值，在出流時(shí)段由下式計(jì)算得

（5）懸沙濃度計(jì)算模型固體邊界條件：一般假定固邊界處，法向泥沙通量為零，即。

1.3 數(shù)值計(jì)算方法與網(wǎng)格生成方法

數(shù)值計(jì)算采用中心有限體積法?？紤]到非結(jié)構(gòu)網(wǎng)格具有節(jié)點(diǎn)隨意編號(hào)、網(wǎng)格隨意加密以提高精度、可以對(duì)復(fù)雜的固岸邊界和水下地形進(jìn)行比較精確的模擬、少量的節(jié)點(diǎn)數(shù)可得到滿意的精度等優(yōu)點(diǎn)，故在數(shù)值計(jì)算中網(wǎng)格采用非結(jié)構(gòu)網(wǎng)格。

1.4 水動(dòng)力模型驗(yàn)證

以洋山港工程水域?yàn)樗畡?dòng)力模型驗(yàn)證對(duì)象。洋山深水港位于上海市南匯嘴東南海域的崎嶇列島海域，該工程近十幾年積累了豐富且可靠性較高水動(dòng)力、泥沙現(xiàn)場(chǎng)同步觀測(cè)資料。模型范圍為東經(jīng)121.8°～122.3°，北緯30.3°～30.8°，該范圍計(jì)算網(wǎng)格共剖分為35 394個(gè)。

水動(dòng)力模型驗(yàn)證包括計(jì)算水域潮位、流速及流向驗(yàn)證等內(nèi)容。圖1給出了潮位計(jì)算值與實(shí)測(cè)值的比較示例，圖1結(jié)果表明，不同位置潮位實(shí)測(cè)值和計(jì)算值吻合良好。圖2給出了流速、流向計(jì)算值與實(shí)測(cè)值的比較示例，圖2結(jié)果表明，不同位置流速及流向的實(shí)測(cè)值和計(jì)算值的吻合程度很高。模擬計(jì)算結(jié)果較好地再現(xiàn)了洋山港區(qū)非正規(guī)半日潮、水流流動(dòng)往復(fù)性較強(qiáng)的特點(diǎn)。

2 泥沙參數(shù)取值及對(duì)懸沙濃度的影響

2.1 懸沙濃度計(jì)算的主要影響參數(shù)篩選

如前所述，影響泥沙運(yùn)動(dòng)的參數(shù)包括泥沙沉降速度、泥沙擴(kuò)散系數(shù)、控制侵蝕速度的比例值、軟底床的侵蝕系數(shù)、硬底床的侵蝕指數(shù)、床面糙率、河床淤積物的干密度、床層間的泥沙轉(zhuǎn)換率。

對(duì)泥沙沉降速度，理論上與泥沙粒度有關(guān)。對(duì)于泥沙粒徑d<0.01 mm的泥沙，將發(fā)生絮凝現(xiàn)象，絮凝沉速在0.000 3～0.000 6 m/s，可直接取絮凝沉速為0.000 4 m/s；而泥沙粒徑在0.03

表1 泥沙沉降速度Tab.1 Velocity of sediment

對(duì)于二維模型中泥沙水平擴(kuò)散系數(shù)ε，一般是基于渦粘系數(shù)類比公式計(jì)算而得，即ε=λεt（εt為渦粘系數(shù)），渦粘系數(shù)比例因子在0.5～1.5。

在泥沙的沉降速度已確定的條件下，分別計(jì)算了λ為0.5、1.0、1.5時(shí)的懸沙濃度值，結(jié)果幾乎沒有什么變化，表明ε對(duì)懸沙濃度計(jì)算影響很小，數(shù)值計(jì)算時(shí)可取渦粘系數(shù)比例因子λ為1.0。

對(duì)于控制侵蝕速度的比例值E，硬泥層一般取0.000 1（kg·m-2）/s，軟泥層可0.000 005～0.000 02（kg·m-2）/s取值，分別取軟底床的為0.000 005（kg·m-2）/s、0.000 01（kg·m-2）/s、0.000 015（kg·m-2）/s、0.000 02（kg·m-2）/s 4 組值進(jìn)行懸沙濃度的計(jì)算，結(jié)果表明懸沙濃度變化很小，即表明E對(duì)懸沙濃度計(jì)算的影響不大，可直接取值0.000 02（kg·m-2）/s。

對(duì)于硬底床的侵蝕指數(shù)n，一般取10；對(duì)于軟底床的侵蝕系數(shù)α通常在4～20變化，對(duì)于淤積偏重的地區(qū)可選擇下限值4，侵蝕偏重的地區(qū)取20，而對(duì)于沖淤相當(dāng)?shù)目稍谥虚g挑選一個(gè)值即可。

綜上所述，泥沙沉降速度ω、泥沙擴(kuò)散系數(shù)ε、控制侵蝕速度的比例值E、軟底床的侵蝕系數(shù)α、硬底床的侵蝕指數(shù)n，在懸沙濃度的數(shù)值計(jì)算中分別按照上述方法取值，對(duì)結(jié)果的影響不大。由此篩選出對(duì)懸沙濃度計(jì)算影響最大的3個(gè)參數(shù)分別為：床面糙率、河床淤積物的干密度以及床層間的泥沙轉(zhuǎn)換率。

下面分別討論床面糙率、河床淤積物的干密度、床層間的泥沙轉(zhuǎn)換率對(duì)懸沙濃度計(jì)算結(jié)果的影響。

2.2 床面糙率k

床面糙率是描述海床對(duì)水流阻力的參數(shù)，可采用謝才系數(shù)、曼寧系數(shù)、尼古拉茲糙率等表示。當(dāng)采用尼古拉茲糙率時(shí)，床面糙率依賴于床面形狀和底床泥沙顆粒的粒徑。進(jìn)行含沙量計(jì)算時(shí)，可將計(jì)算范圍的局地床面形狀均值化，即床面糙率僅依賴底床泥沙顆粒的粒徑，其值可按n倍泥沙顆粒的粒徑計(jì)算，即k=nd。

在擴(kuò)散系數(shù)ε等其他5個(gè)參數(shù)固定不變、且按照2.1節(jié)所述取值的條件下，分別選取床面糙率k=1.0 d、2.5 d、4.0 d、5.5 d，分析床面糙率k的取值變化對(duì)懸沙濃度的影響。

根據(jù)經(jīng)驗(yàn)，當(dāng)d<0.5 mm時(shí)，有些學(xué)者直接取k=0.001 m。對(duì)于洋山港工程水域，底質(zhì)平均中值粒為0.037 mm。因此，在討論床面糙率對(duì)懸沙濃度的影響時(shí)，也將k=0.001 m作為一種處理方法，與k=nd的處理方法進(jìn)行比較。圖3給出了采用不同的床面糙率時(shí)懸沙濃度的計(jì)算結(jié)果。

由圖3可見，對(duì)于床面糙率，若直接取0.001 m，計(jì)算得到的懸沙濃度隨時(shí)間的變化過程與實(shí)測(cè)結(jié)果差別較大。若床面糙率按照粒徑的倍數(shù)來計(jì)算，當(dāng)糙率取不同的粒徑倍數(shù)時(shí)，計(jì)算的懸沙濃度和實(shí)測(cè)濃度的誤差也將隨之變化（表2）。

將表2的不同粒徑倍數(shù)情況下所計(jì)算得到的平均懸沙濃度和實(shí)測(cè)的平均懸沙濃度進(jìn)行對(duì)比，將粒徑倍數(shù)值與誤差關(guān)系繪出曲線（圖4）發(fā)現(xiàn)，當(dāng)選取k≈2 d時(shí)，計(jì)算的平均懸沙濃度和實(shí)測(cè)濃度的誤差最小，即床面糙率的最佳取值應(yīng)為k≈2 d。

表2 床面糙率對(duì)平均懸沙濃度值的影響匯總Tab.2 Influence of bed roughness on value of suspended sediment concentration

2.3 河床淤積物的干密度

河床淤積物的干密度是一個(gè)重要的泥沙參數(shù)，對(duì)于沖淤體積和重量的換算、泥沙起動(dòng)、河床沖刷及水流挾沙能力等均有相當(dāng)重要的影響［6］，干密度是不隨時(shí)間變化的。河床的最上層為最弱層，主要是浮泥或新淤積的泥沙，稱之為軟泥；其下面層的密度及強(qiáng)度都不斷遞增，稱之為硬泥。

河床淤積物的干容量與孔隙率成反比，孔隙率與泥沙粒徑成反比，河床淤積物的干容重隨著時(shí)間的變化很小，根據(jù)大量的實(shí)驗(yàn)得出：軟泥的干密度［7］會(huì)在100～400 kg/m3浮動(dòng)，具體的數(shù)值是取決于各自的環(huán)境，和新淤積泥沙的濕容重也有密不可分的關(guān)系，而硬泥則是由底質(zhì)粒徑按照公式γd=1 750×d0.183［8］計(jì)算，其中粒徑的單位為mm。

港區(qū)的海床分為兩層：上層為軟底床，即泥沙最新輸運(yùn)落淤的床層，分別計(jì)算了干密度為150 kg/m3、250 kg/m3、400 kg/m3這3種不同的密度量級(jí)（圖5）。區(qū)域中新落淤的泥沙干密度取250 kg/m3時(shí)計(jì)算的泥沙濃度和實(shí)測(cè)值擬合較好，其他的誤差在10%～60%不等；下層為硬底床，按照底層粒徑和公式γd=1 750×d0.183計(jì)算干密度。

2.4 床層間的泥沙轉(zhuǎn)換率

在模擬計(jì)算時(shí)，通常將海床分為可動(dòng)層和固結(jié)層兩層，層間的泥沙在動(dòng)力作用下的泥沙交換量的多少采用泥沙轉(zhuǎn)換率描述。該參數(shù)對(duì)于研究泥沙運(yùn)動(dòng)機(jī)理、考慮床沙和運(yùn)動(dòng)泥沙的關(guān)系、計(jì)算懸沙濃度及海床沖淤強(qiáng)度等均為最重要影響參數(shù)之一。其取值范圍一般在0.1～0.001（kg·m-2）/s，該范圍上、下限量值上相差100倍，各個(gè)取值對(duì)懸沙濃度計(jì)算結(jié)果的影響是顯而易見的。

在此，與分析床面糙率k對(duì)懸沙濃度影響一樣，在探討泥沙轉(zhuǎn)換率對(duì)懸沙濃度影響時(shí)，也將擴(kuò)散系數(shù)ε等其他 5 個(gè)參數(shù)固定不變（按照 2.1 節(jié)所述取值），分別選取 m 為 0.000 1（kg·m-2）/s、0.001（kg·m-2）/s、0.01（kg·m-2）/s、0.1（kg·m-2）/s，計(jì)算了 4 組不同的轉(zhuǎn)換率對(duì)應(yīng)的懸沙濃度隨時(shí)間的變化過程，與實(shí)測(cè)結(jié)果的比較匯總于圖6。

由圖 6 可知，當(dāng) m 取值為 0.01（kg·m-2）/s和 0.1（kg·m-2）/s時(shí)，即固結(jié)率較高時(shí)，由于大部分懸沙落於到床面后固結(jié)，使得懸沙濃度大幅度降低。換言之，懸沙濃度模擬計(jì)算時(shí)，如果過高估計(jì)了泥沙轉(zhuǎn)換率，將使得計(jì)算懸沙濃度小于實(shí)際懸沙濃度。反之，計(jì)算懸沙濃度將大于實(shí)際懸沙濃度。計(jì)算表明，當(dāng)泥沙轉(zhuǎn)換率為0.001（kg·m-2）/s時(shí)，模擬的泥沙濃度和實(shí)測(cè)的泥沙濃度值擬合最好。即最佳床層間的泥沙轉(zhuǎn)換率為0.001（kg·m-2）/s。

3 懸沙濃度計(jì)算在洋山港實(shí)際工程中的驗(yàn)證

對(duì)于影響懸沙濃度計(jì)算的各個(gè)泥沙參數(shù)，采用前節(jié)泥沙參數(shù)的取值方法，分別取泥沙沉降速度：絮凝沉速ω=0.000 4 m/s，非絮凝沉速按表1取值；泥沙擴(kuò)散系數(shù)：渦粘系數(shù)比例因子λ=1.0；控制侵蝕速度的比例值：硬泥層 E=0.000 1（kg·m-2）/s，軟泥層 E=0.000 02（kg·m-2）/s；硬底床的侵蝕指數(shù)：n=10；軟底床的侵蝕系數(shù)：對(duì)于港區(qū)南部（淤積）α=4，對(duì)于港區(qū)北部（沖刷）α=20；床面糙率：k=2 d；泥沙干密度：對(duì)軟底床γd=250 kg/m3；對(duì)硬底床γd=1 750×d0.183；床層間的泥沙轉(zhuǎn)換率：m=0.001（kg·m-2）/s。

選定上述泥沙參數(shù)后，對(duì)洋山港全水域各特征代表點(diǎn)的懸沙濃度分別進(jìn)行數(shù)值計(jì)算，并與實(shí)測(cè)結(jié)果進(jìn)行對(duì)比。選取的代表點(diǎn)位置分布在洋山港外南北附近、港區(qū)、汊道等典型水域，代表點(diǎn)的具體位置見圖7。圖8給出了各代表點(diǎn)垂線平均懸沙濃度計(jì)算和實(shí)測(cè)結(jié)果比較示例。

圖 8-a 為小洋山北部 N2 點(diǎn)（北緯°39′26.54″，東經(jīng) 2°04′04.85″）大潮期間垂線平均含沙量時(shí)間變化過程實(shí)測(cè)結(jié)果和計(jì)算結(jié)果的比較示例；圖 8-b～圖 8-f依次為大洋山南部 S1 點(diǎn)（30°33′19.14″，122°04′00.86″）、三期碼頭中部 M11 點(diǎn)（30°36′39.82″，122°04′30.31″）、西口門中部 W2 點(diǎn)（30°37′32.10″，122°00′13.05″）、大山塘—大貼餅汊道 SC3 點(diǎn)（30°35′30.58″，122°02′25.58″）、以及小洋山一顆珠山汊道 K2 點(diǎn)（30°38′41.28″，122°02′23.38″）的大潮期間垂線平均含沙量時(shí)間變化過程計(jì)算結(jié)果和實(shí)測(cè)結(jié)果的比較示例。表3給出了上述各代表點(diǎn)大潮和小潮期間的時(shí)均垂線平均含沙量計(jì)算與實(shí)測(cè)值的比較。

表3 洋山港工程水域典型泥沙測(cè)點(diǎn)時(shí)均垂線平均懸沙濃度值的計(jì)算值與實(shí)測(cè)值比較Tab.3 Comparison between computed vertical suspended sediment concentration in the average time and measured results at representative point in the Yangshan Port

由圖8和表3的比較可知，在洋山港周圍水域，按照前述推薦的各參數(shù)取值，計(jì)算得到的含沙量濃度隨時(shí)間變化過程及在一個(gè)潮周期內(nèi)的時(shí)均值與實(shí)測(cè)結(jié)果的吻合程度是令人滿意的。

4 結(jié)語

基于經(jīng)過實(shí)測(cè)資料驗(yàn)證的二維水動(dòng)力模型，在較準(zhǔn)確模擬潮流動(dòng)力的基礎(chǔ)上，對(duì)于影響懸沙濃度計(jì)算有關(guān)的8個(gè)參數(shù)通過數(shù)值計(jì)算的方法：首先篩選出影響最大的3個(gè)參數(shù)為床面糙率、河床淤積物的干密度、床層間的泥沙轉(zhuǎn)換率；然后在其各個(gè)參數(shù)取值的可變范圍內(nèi)，分別調(diào)整各個(gè)參數(shù)的取值，計(jì)算對(duì)比得到上述3個(gè)參數(shù)的最優(yōu)取值。最后，就洋山港實(shí)際工程進(jìn)行了懸沙濃度模擬計(jì)算，驗(yàn)證了所選擇參數(shù)的可靠性。

研究表明，對(duì)于以潮流為主要環(huán)境動(dòng)力的淤泥質(zhì)海床而言，泥沙沉降速度可通過物理模型試驗(yàn)或較為成熟的經(jīng)驗(yàn)公式計(jì)算得到較為可靠的數(shù)值；泥沙擴(kuò)散系數(shù)，控制侵蝕速度的比例值，軟底床的侵蝕系數(shù)，硬底床的侵蝕指數(shù)等對(duì)懸沙濃度的影響較小。對(duì)于床面糙率、軟底床淤積物的干密度、床層間的泥沙轉(zhuǎn)換率，最佳取值分別為 2 d、250 kg/m3、0.001（kg·m-2）/s。

［1］ZUO S H，ZHANG N C，LI B，et al.Numerical Simulations of Tidal Current&Sediment and Sea Bed Erosion and Deposition for Yangshan Deep-water Harbor of Shanghai［J］.International Journal of Sediment Research，2009，24（3）：287-298.

［2］左書華，張寧川，張征.島群海域環(huán)境下泥沙運(yùn)動(dòng)及地形沖淤變化數(shù)值模擬研究［J］.泥沙研究，2011（2）：1-8.ZUO S H，ZHANG N C，ZHANG Z.Numerical modeling of sediment transport and seabed erosion and siltation in archipelago sea area［J］.Journal of Sediment Research，2011（2）：1-8.

［3］李文丹，李孟國(guó)，龐啟秀.臺(tái)山核電站取排水工程潮流泥沙數(shù)值模擬研究［J］.水道港口，2011，32（2）：94-100.LI W D，LI M G，PANG Q X.The tide and sediment numerical simulation study of Taishan nuclear power plant project［J］.Journal of Waterway and Harbor，2011，32（2）：94-100.

［4］XIN W J.Computational Techniques of 2D Tidal Flow in Estuaries and Bays［J］.China Ocean Engineering，1995，9（4）：395-404.

［5］曹祖德，楊樹森，楊華.粉沙質(zhì)海岸的界定及其泥沙運(yùn)動(dòng)特點(diǎn)［J］.水運(yùn)工程，2003，352（5）：15-20.CAO Z D，YANG S S，YANG H.Definition of silt-sandy beach and its characteristics of sediment movement［J］.Port&Waterway Engineering，2003，352（5）：15-20.

［6］潘慶燊.三峽工程泥沙問題研究［M］.北京：中國(guó)水利水電出版社，1999.

［7］韓其為，王玉成，向熙瓏.淤積物的初期干容重［J］.泥沙研究，1981（1）：1-6.HAN Q W，WANG Y C，XIANG X L.Initial specific weight of deposits［J］.Journal of Sediment Research，1981（1）：1-6.

［8］羅肇森.潮汐通道口攔門沙航道的淤積計(jì)算［J］.海洋工程，1992（3）：10-18.LUO Z S.Prediction of sedimentation for the navigation channel of a tidal inlet with mouth bay［J］.Ocean Engineering，1992（3）：10-18.