二維潮流泥沙數(shù)學(xué)模型 在挖入式港池中的應(yīng)用研究

談寶林 王俊虎

摘 要：本文根據(jù)長(zhǎng)江口南港水域歷年水文及地形資料建立的長(zhǎng)江口二維潮流泥沙數(shù)學(xué)模型，模擬分析計(jì)算中船長(zhǎng)興造船基地二期工程挖入式港池方案實(shí)施后，周邊水域的水動(dòng)力變化以及港池泥沙淤積強(qiáng)度，港池內(nèi)水體的交換能力，為工程項(xiàng)目設(shè)計(jì)及后續(xù)運(yùn)營(yíng)提供技術(shù)支持，為類似項(xiàng)目提供參考。

關(guān)鍵詞：長(zhǎng)興二期;挖入式港池;泥沙回淤;水體交換

中圖分類號(hào)：U653.3? ? ? ? ? 文獻(xiàn)標(biāo)識(shí)碼：A? ? ? ? ? ? 文章編號(hào)：1006—7973（2021）12-0120-04

1項(xiàng)目背景

中船長(zhǎng)興造船基地二期工程位于上海市崇明區(qū)長(zhǎng)興鎮(zhèn)，西起躍進(jìn)港以東、東至橫沙小港以西，南臨長(zhǎng)江。岸線自上游至下游依次為挖入式港池、2#舾裝碼頭、2#船塢、3#舾裝碼頭、4#舾裝碼頭。

本次擬建挖入式港池：在用地上游最西側(cè)利用海慶河入海口現(xiàn)狀，港池大致成EN-WS走向，港池寬度240m。挖入式港池內(nèi)布置兩座舾裝碼頭，其中港池西側(cè)為5#舾裝碼頭，長(zhǎng)度720m、寬度為20m;港池東側(cè)為6#舾裝碼頭，長(zhǎng)度為760m、寬度為20m。港池北側(cè)布置1座室內(nèi)船塢，其余為工作船碼頭。

2建立長(zhǎng)江口潮流泥沙數(shù)學(xué)模型

2.1計(jì)算方法

2.1.1二維淺水控制方程

在笛卡爾直角坐標(biāo)系下，根據(jù)靜壓和勢(shì)流假定，沿垂向平均的二維潮流基本方程可表述為如下形式：

（1）

（2）

（3）

（4）

（5）

2.1.2基本方程的離散及求解

方程（1～4）垂向積分變成二維形式，寫成向量形式：

（2-6）

2.1.3邊界處理

邊界分為開邊界和閉邊界。由于本文采用的是有限體積法，水位、流速布置在網(wǎng)格中心點(diǎn)，網(wǎng)格邊界上沒有布置變量，因此不能夠通過網(wǎng)格邊界處理邊界條件，需用到特殊的邊界處理方法。

2.1.3.1開邊界

對(duì)于邊界處的網(wǎng)格，可求，關(guān)鍵是求解，開邊界又分為急流開邊界和緩流開邊界，因本文模型為緩流開邊界，故下文只討論緩流開邊界，根據(jù)相容關(guān)系處理。

2.1.3.2閉邊界

采用鏡像法處理。在閉邊界外側(cè)虛擬一個(gè)單元，邊界上的兩側(cè)的法向流速相反，切向流速相同，即表示單元法向和切向流速。

2.1.4關(guān)鍵技術(shù)問題

2.1.4.1動(dòng)邊界的處理

本次研究采用凍結(jié)法，根據(jù)計(jì)算單元水深判斷是否露灘，當(dāng)水深小于某一控制水深時(shí)，單元潮位“凍結(jié)”不變，要進(jìn)行下一時(shí)刻計(jì)算前，被凍結(jié)的單元水深由周邊有效水深進(jìn)行修正，如果水深大于控制水深則重新參與計(jì)算，為避免水量不平衡，動(dòng)邊界控制水深采用0.01m。

2.1.4.2糙率處理

糙率是潮流計(jì)算的主要計(jì)算參數(shù)之一，反映了潮流運(yùn)動(dòng)過程中的阻力特性，糙率選取正確與否對(duì)計(jì)算結(jié)果有直接影響。糙率是一個(gè)綜合參數(shù)，與床面泥沙特性、水深及地形形態(tài)都有一定關(guān)系，本項(xiàng)研究中根據(jù)經(jīng)驗(yàn)選用了附加糙率公式，考慮水深變化后的糙率響應(yīng)。

2.1.4.3樁基模擬

由于本次挖入式港池兩側(cè)建有舾裝碼頭，結(jié)構(gòu)為高樁承臺(tái)，因此本報(bào)告樁基的模擬采用概化方式。一般數(shù)學(xué)模型中對(duì)于樁基模擬常用的方法有加密網(wǎng)格法、附加糙率法、附加阻力法和碼頭區(qū)不過水法等，本次模擬中為更好地模擬樁基引起的阻力，采用碼頭區(qū)地形概化和加糙相結(jié)合的方法，首先將樁群看作是阻水建筑物，將碼頭平臺(tái)阻水部分當(dāng)作斷面突然縮水的建筑物考慮， 進(jìn)而，根據(jù)樁徑計(jì)算樁基的阻水體積，將其換算至樁所在的網(wǎng)格單元，計(jì)算網(wǎng)格地形凸起高度，修改網(wǎng)格地形。

2.2數(shù)學(xué)模型建立

本次數(shù)學(xué)模型采用嵌套模型計(jì)算方式。大模型為中國(guó)近海模型，嵌套模型為長(zhǎng)江口近海模型，長(zhǎng)江口近海模型的邊界條件由中國(guó)近海模型提供。

2.2.1 中國(guó)近海模型

中國(guó)近海模型，南邊界至越南的Qui Nhon—馬來西亞的Subic Bay連線，臺(tái)灣島，東南側(cè)邊界沿巴布延群島—八示戈島—石桓島—沖繩群島，并與日本的南九州地區(qū)相連，東側(cè)邊界由日本的福岡與韓國(guó)的釜山相連。

模型采用三角形網(wǎng)格作為計(jì)算單元，網(wǎng)格總數(shù)為271156個(gè)，最小網(wǎng)格邊長(zhǎng)約為244m，最大網(wǎng)格邊長(zhǎng)約為28598m。

2.2.2長(zhǎng)江口近海模型

局部小模型包括整個(gè)長(zhǎng)江口和杭州灣，數(shù)學(xué)模型上邊界至南京，外海邊界至長(zhǎng)江口口外-40m等深線左右。模型長(zhǎng)江口南京至徐六涇采用2016年實(shí)測(cè)地形資料，北支采用2010年實(shí)測(cè)地形資料，南支、南港、北港、南槽、北槽采用2019年最新實(shí)測(cè)地形資料，外海地形采用最新海圖拼接。

模型采用三角形網(wǎng)格，共計(jì)劃分單元122224個(gè)，節(jié)點(diǎn)總數(shù)62877個(gè)，最小網(wǎng)格邊長(zhǎng)約15m，最大網(wǎng)格邊長(zhǎng)約為14935m。

2.3模型率定

采用2017年7月10日～7月19日的長(zhǎng)江口洪季實(shí)測(cè)同步水文資料進(jìn)行驗(yàn)證，分別驗(yàn)證潮位、潮流和含沙量過程，均得到了較好的模擬，與實(shí)測(cè)值吻合較好。

2.4水流泥沙驗(yàn)證

采用2019年1月11日～1月30日長(zhǎng)江口枯季實(shí)測(cè)同步水文資料進(jìn)行驗(yàn)證。潮位各站模型計(jì)算潮位過程與實(shí)測(cè)值吻合較好。漲落潮流速偏差在10%以內(nèi)。含沙量偏差在20%以內(nèi)，部分測(cè)點(diǎn)含沙量偏差較大，主要因?yàn)楝F(xiàn)場(chǎng)影響含沙量因素多且復(fù)雜，數(shù)學(xué)模型無法全部考慮。

2.5航道回淤驗(yàn)證

以2016年北槽航道淤積資料為地形驗(yàn)證資料，回淤驗(yàn)證外海潮型采用大-中-小潮組合潮型，上游流量采用多年洪季平均流量和多年枯季平均流量組合的方式。具體計(jì)算方式為：“多年洪季平均流量+大-中-小潮潮型”計(jì)算30天，“多年枯季平均流量++大-中-小潮潮型”計(jì)算30天，然后采用時(shí)間加速計(jì)算航道年回淤量。

2016年北槽航道淤積量為5400萬方，數(shù)學(xué)模型計(jì)算量為6321萬方，O～V單元數(shù)學(xué)模型計(jì)算值較大，總回淤量偏差約為17%?？傮w來說，數(shù)學(xué)模型計(jì)算的航道淤積分布能夠反映北槽的航道淤積情況。

建立的長(zhǎng)江口大范圍二維潮流泥沙數(shù)學(xué)模型，以2017年、2019年實(shí)測(cè)水文、航道回淤資料作為水動(dòng)力、泥沙率定驗(yàn)證資料，結(jié)果滿足規(guī)程要求，說明該數(shù)學(xué)模型，能夠反映長(zhǎng)江口的水沙運(yùn)動(dòng)特征。

3水動(dòng)力變化分析

3.1計(jì)算條件

外海潮汐條件采用2017年洪季驗(yàn)證邊界條件，采用2018年長(zhǎng)興站的實(shí)測(cè)潮位觀測(cè)資料統(tǒng)計(jì)，累積頻率為10%的典型大潮潮差約為3.4m，累積頻率為90%的典型小潮潮差約為1.5m，數(shù)學(xué)模型計(jì)算的長(zhǎng)興站大潮潮差約為3.3m，略小于典型大潮，長(zhǎng)興站大潮潮差約為1.45m，基本接近典型小潮。

水動(dòng)力模型計(jì)算外海邊界采用大潮潮型，徑流采用長(zhǎng)江多年洪季平均徑流;泥沙模型計(jì)算外海邊界采用大中小組合潮型，徑流分別采用長(zhǎng)江多年洪季、枯季平均徑流;水體交換模型計(jì)算外海邊界分別采用大潮和小潮潮型，徑流采用長(zhǎng)江多年洪季平均徑流。

3.2潮位變化

沿著岸線，在港池口門處取一點(diǎn)，港池西北側(cè)取一點(diǎn)，港池的東南側(cè)位依次取三點(diǎn)，作為潮位變化分析采樣點(diǎn)。經(jīng)模擬分析，從潮位過程變化看，工程前后潮位過程基本沒有變化，從高低潮位變化看，工程方案實(shí)施后，長(zhǎng)興水道的高低潮位基本不變，因此方案的實(shí)施對(duì)周邊潮位基本沒有影響。

3.3流態(tài)變化

經(jīng)模擬分析本底條件下，漲落潮流向受南港岸線的約束，潮流為往復(fù)流特征，漲潮方向?yàn)闁|南向西北，落潮方向相反，為西北向東南，由于工程區(qū)附近的長(zhǎng)興水道等水域水深條件較好，因此，漲落潮流較平順，但在擬建挖入式港池附近形成局部的回流。

經(jīng)模擬分析方案實(shí)施后，南港及長(zhǎng)興水道大范圍的潮流流態(tài)變化相對(duì)較小，與本底基本一致。挖入式港池由于水域的開挖作用，漲落潮流態(tài)呈現(xiàn)出明顯的回流特征，具體表現(xiàn)在：漲潮時(shí)，在東南向西北的漲潮流作用下，港池口門處出現(xiàn)順時(shí)針的回流區(qū)，同時(shí)港池內(nèi)部在口門回流的作用下，形成逆時(shí)針方向的二次回流;落潮時(shí)，在西北向東南向的落潮流作用下，港池內(nèi)形成大范圍的逆時(shí)針回流。與漲潮回流不同的是，漲潮回流為二次回流，而落潮回流未出現(xiàn)二次回流，主要原因?yàn)槁涑眲?dòng)力較強(qiáng)，造成港池口門區(qū)的水流流速相對(duì)較大，因而回流的范圍增大。

3.4流場(chǎng)變化

經(jīng)模擬分析挖入式港式方案實(shí)施后，對(duì)南港主槽及長(zhǎng)興水道水域潮流基本不影響。港池內(nèi)部的流速變化是由于港池開挖后，其內(nèi)部自身的水動(dòng)力引起，而非與本底的流速對(duì)比變化。

3.5港池及近岸水流條件分析

經(jīng)模擬分析，5#、6#舾裝碼頭前沿流速基本在0.15～0.20m/s之間，流速相對(duì)較弱，口門處流速相對(duì)較大，漲落急流速可達(dá)到0.60～0.70m/s左右。對(duì)于挖入式港池內(nèi)部，自口門向港池內(nèi)部，漲落急流速均呈逐漸減小趨勢(shì)，對(duì)于口門附近，漲落急流速基本在0.20m/s左右，而港池內(nèi)部漲落急流速在0.05～0.10m/s左右，總體較小;港池口門處，漲落急流速可達(dá)到0.60～0.70m/s左右;港池上游，漲急流速約為0.90m/s左右，落急流速在1.30m/s左右，港池下游漲急流速基本在1.00m/s左右，落急流速在1.20～1.50m/s左右，落急流速總體大于漲急。碼頭前沿最大橫流均在0.10m/s以內(nèi)，幅度較小。

4挖入式港池回淤計(jì)算

4.1含沙量分布

經(jīng)模擬統(tǒng)計(jì)分析，港池自口門向底部，含沙量呈逐漸減小趨勢(shì)，但變化幅度相對(duì)較小，口門區(qū)平均含沙量約為0.23kg/m3，港池末端含沙量約為0.18kg/m3。

4.2泥沙回淤分析

經(jīng)模擬統(tǒng)計(jì)分析，港池口門處回淤強(qiáng)度較大，往港池內(nèi)部，回淤強(qiáng)度逐漸減小，口門處最大回淤強(qiáng)度約為2.01m/a，港池末端回淤強(qiáng)度約為0.76m/a，挖入式港池內(nèi)的平均回淤強(qiáng)度約為1.29m/a，回淤量約為22.9萬方/a。5#舾裝碼頭前沿淤積強(qiáng)度約為1.07m/a，6#舾裝碼頭前沿淤積強(qiáng)度約為1.27m/a，港池南側(cè)淤積強(qiáng)度總體大于北側(cè)。

4.3大風(fēng)天驟淤計(jì)算

采用麥德遜（Madsen）和格蘭特（Grant）的波浪推移質(zhì)輸沙率公式計(jì)算底沙驟淤，當(dāng)海向大風(fēng)（強(qiáng)臺(tái)風(fēng)）出現(xiàn)時(shí)，含沙量可能呈數(shù)倍增大，采用劉家駒公式作近似推算。經(jīng)計(jì)算，港池內(nèi)一次大風(fēng)過程（50年一遇S向波浪作用12h），推移質(zhì)驟淤強(qiáng)度約為0.29m/次，懸移質(zhì)驟淤約為0.05m/次，因此港池內(nèi)一次大風(fēng)過程驟淤強(qiáng)度約為0.34m/次。

5挖入式港池水體交換分析

本文采用數(shù)學(xué)模型示蹤劑方法，初始計(jì)算時(shí)，在模擬水域內(nèi)充滿濃度為1的示蹤劑，連續(xù)計(jì)算50h，分析港池內(nèi)示蹤劑濃度的變化。水體交換能力評(píng)價(jià)指標(biāo)主要有水體交換率、半交換周期、更替周期等，本文采用半交換周期分析水體交換能力。

經(jīng)模擬計(jì)算分析，大潮條件下港池內(nèi)的水體交換速度明顯大于小潮，大潮港池內(nèi)的水體交換半周期約為4.5h，2/3水體交換周期約為8h，小潮港池內(nèi)的水體交換半周期約為7.5h，2/3水體交換周期約為12.5h。

6結(jié)語

在分析南港歷年水文及地形資料的基礎(chǔ)上，建立了長(zhǎng)江口二維潮流泥沙水質(zhì)數(shù)學(xué)模型，以2017年洪季、2019年枯季同步水文測(cè)驗(yàn)資料為率定和驗(yàn)證資料，模型相似性良好，能夠復(fù)演南港水域的水沙運(yùn)動(dòng)。在此基礎(chǔ)上，計(jì)算了挖入式港池方案實(shí)施后，周邊水域的水動(dòng)力變化以及港池泥沙淤積強(qiáng)度，分析了港池內(nèi)水體的交換能力，主要結(jié)論如下：

（1）長(zhǎng)江口南港水域潮汐類型為正規(guī)半日潮，M2分潮占絕對(duì)主導(dǎo)地位，漲潮流向一般為295o～313o，落潮流向?yàn)?15o～122o。據(jù)1999～2019年南港實(shí)測(cè)水文測(cè)驗(yàn)統(tǒng)計(jì)，南港漲、落潮平均流速未見明顯趨勢(shì)性變化。采用2018年長(zhǎng)興站水位資料計(jì)算了臺(tái)風(fēng)“安比”引起的增水值，經(jīng)計(jì)算，安比臺(tái)風(fēng)在長(zhǎng)興站引起的最大增水約為0.5m。

（2）南港河段漲潮平均含沙量大于落潮平均含沙量，三峽水庫蓄水后，漲潮和落潮平均含沙量降低較明顯。南港懸沙中值粒徑大都在0.008～0.014mm之間，南港底質(zhì)以細(xì)砂、極細(xì)沙為主，中值粒徑大都大于0.062mm。

（3）南港河床地形變化與上游沙體運(yùn)動(dòng)密切相關(guān)，河槽斷面形態(tài)呈“U”“W”型交替變化;1998～2019年，工程區(qū)附近長(zhǎng)興水道水深條件總體較好，深槽水深基本大于10m。

（4）工程方案實(shí)施后，工程附近水域高低潮位變化均在0.002m以內(nèi)，幾乎不變，因此方案的實(shí)施對(duì)周邊潮位基本沒有影響;漲潮時(shí)，港池口門處出現(xiàn)順時(shí)針的回流，港池內(nèi)部形成逆時(shí)針方向的二次回流，落潮時(shí)，港池內(nèi)形成大范圍的逆時(shí)針回流，南港主槽及長(zhǎng)興水道處等水域潮流基本不變;挖入式港池碼頭前沿最大橫流均在0.10m/s以內(nèi)，幅度較小。

（5）港池自口門向底部，含沙量呈逐漸減小趨勢(shì)，口門區(qū)平均含沙量約為0.23kg/m3，港池末端含沙量約為0.18kg/m3。港池口門處回淤強(qiáng)度較大，往港池內(nèi)部，回淤強(qiáng)度逐漸減小，口門處最大回淤強(qiáng)度約為2.01m/a，港池末端回淤強(qiáng)度約為0.76m/a，挖入式港池內(nèi)的平均回淤強(qiáng)度約為1.29m/a，回淤量約為22.9萬方/a。

（6）港池內(nèi)一次大風(fēng)過程（50年一遇S向波浪作用12h），推移質(zhì)驟淤強(qiáng)度約為0.29m/次，懸移質(zhì)驟淤約為0.05m/次，港池內(nèi)一次大風(fēng)過程驟淤強(qiáng)度約為0.34m/次。

（7）大潮期，港池內(nèi)的水體交換半周期約為4.5h，2/3水體交換周期約為8h;小潮期，港池內(nèi)的水體交換半周期約為7.5h，2/3水體交換周期約為12.5h。

參考文獻(xiàn)：

[1] JTS/T231-2010，海岸與河口潮流泥沙模擬技術(shù)規(guī)程[S].

[2]羅小峰.中船長(zhǎng)興造船基地一期工程水工建筑物物理模型試驗(yàn)研究[Z].南京水利科學(xué)研究院，2005.

[3]陸永軍，左利欽，王紅川，李浩麟.波浪與潮流共同作用下二維泥沙數(shù)學(xué)模型[J].泥沙研究，2005（06）：3-14.

[4] 韓玉芳.中船長(zhǎng)興造船基地二期工程物理模型試驗(yàn)研究[Z].南京水利科學(xué)研究院，2007.