朱明栓，張思慧，江中林

(1.福建船政交通職業(yè)學院 土木工程學院，福建省福州市倉山區(qū)首山路80號 350007；2.福建省興禹源水利工程設計有限公司 設計部，福建省福州市鼓樓區(qū)福飛南路104號 350001)

閩江是福建的母親河，也是我國重要的入海河流，其發(fā)源于武夷山脈南麓，為山溪性河流，水量豐沛，為福建省內(nèi)第一大江河，橫穿省會城市福州而過注入東海。閩江下游河段屬感潮河段，長期受徑流和潮流的雙重影響，致使河流流態(tài)和形態(tài)復雜，流向多變且順逆不定。閩江下游人類活動頻繁，社會經(jīng)濟發(fā)達，近十幾年來，一方面閩江上游河流電站的梯級開發(fā)，特別是水口電站的建成，使上游來沙量急劇減少，另一方面為滿足經(jīng)濟建設對“閩江沙”進行了大量的開采。上游來沙量的減少和對“閩江沙”的過量開采，極大地改變了天然徑流的水力條件和泥沙的運動規(guī)律，導致閩江下游河床急劇下切。為此，閩江下游河道問題已經(jīng)成為制約社會經(jīng)濟發(fā)展、危害生態(tài)環(huán)境、影響人民生活的突出問題。閩江南港航道整治工程的盡快實施對開發(fā)南港航運資源、減輕福州市區(qū)的防洪壓力、促進福州中心區(qū)域的防洪減壓、促進福州中心區(qū)域閩江旅游資源的開發(fā)與經(jīng)營有直接意義，并有利于閩江下游河沙的可持續(xù)開發(fā)利用和水環(huán)境的改善，對加快實施福州市的城市“東擴南進”戰(zhàn)略、福建海西發(fā)展戰(zhàn)略等都有十分重要的意義。因此迫切需要對該河段的水動力場和泥沙運動規(guī)律等問題開展深入的研究，以便于該段河道整治工程方案及相關工作的順利開展。

對于描述水平尺度遠大于垂直尺度的河流、海洋、湖泊等水體運動狀態(tài)的物理量，當沿水深方向的變化相對于水平方向的變化小得多時，可將三維問題通過沿水深方向積分取平均，得到沿水深積分平均的平面二維簡化形式[1]。文中應用平面二維水流數(shù)學模型研究閩江下游河段水流運動特征，為該段河道整治及感潮河流泥沙分布情況研究提供參考。

1 建模區(qū)域河道概況

閩江自淮安起，下游河道分為南北港，分別繞南臺島兩側流向下游[2]。閩江南港也稱烏龍江，長約34km，位于福州市內(nèi)南臺島南面，中段有大樟溪水流匯入，在馬尾處與北港匯合，繼而折向北進入通海河段，流經(jīng)青洲、閩安峽谷至亭江注入東海，長約12km。河道坡降較上游變小，河水流速減緩，并受潮水頂托作用，致使河道曲折而寬淺，沉積作用顯著，沙洲、邊灘發(fā)育，泥沙淤積嚴重。南港是天然的泄洪排沙水道，接受大部分中游來沙，江面寬淺，大部分江面寬為1000～3000m，最寬處達4000m以上。模型建模區(qū)域為南港入口淮安至閩江入?？谕そ拥谰唧w情況如圖1所示。

2 二維水流數(shù)學模型

2.1 基本方程

連續(xù)性方程:

運動方程[3]:

2.2 定解條件與求解

初始條件：

z(x，y，t)|t=0=zA(x，y)

u(x，y，t)|t=0=uA(x，y)

v(x，y，t)|t=0=vA(x，y)

邊界條件：

對水邊界Γ1

z(x，y，t)|(x，y)Γ1=zB(x，y，t)|(x，y)∈Γ1

對岸邊界Γ2

u(x，y)cos(n，x)+v(x，y)cos(n，y)=0

其中，Γ1為水邊界，Γ2為岸邊界；cos(n，x)，cos(n，y)為外法線的方向余弦。?Ω=Γ1+Γ2，Ω為水、岸邊界圍成的平面區(qū)域[4-5]。

文中采用有限元中的Galerkin法來求解平面二維河道水流方程[6]，對控制方程進行離散，設置變量并用Fortran語言編制計算程序。

3 模型中關鍵問題的處理

3.1 動邊界的處理

在天然河道及河口處通常存在較大面積的沙脊和潮灘，它們隨著水位的漲落時而裸露時而淹沒，相應的水域面積也隨著減小或增大，因此水邊界的變化幅度比較大。在對此區(qū)域進行數(shù)值模擬計算時，為了較為真實準確地模擬落潮歸槽、漲潮漫灘的潮流運動，使計算持續(xù)下去，有必要引進動邊界技術。

建模區(qū)域位于閩江南港，此處江面開闊，河道流速減緩，且受潮水頂托作用，灘涂、沙洲發(fā)育良好，漲潮、落潮時的情況不一，從淺灘流態(tài)和計算穩(wěn)定性來考慮需要做動邊界處理。當前對動邊界的處理方法有很多，常見的有水位判別法、凍結法、開挖法、切削法、窄縫法和線邊界法。因為水位判別法具有清晰的物理概念、簡單的實現(xiàn)過程和良好的計算效果，所以廣泛推廣應用[7-8]。同時，水位判別法還能對灘槽、露灘、沙脊等多種地貌并存的河流獲得較高精度的潮流、潮汐全貌，故文中采用此法。

在實際計算中一般把潮灘區(qū)的動邊界按二維問題來處理。對于平面上任一潮灘點(i，k)，瞬時水深hi，k=Di，k+ζi，k，其中Di，k，ζi，k分別為(i，k)點的灘地高程及潮位。

落潮過程中，ζi，k<0，當hi，k≤hmin時，認為該點干出，不參與計算，且令其流速為0；反之，則認為該點淹沒，參與計算。漲潮過程中，水邊界線隨潮位的上升向高灘推進，則計算網(wǎng)格點增多。因新增網(wǎng)格點原無潮位值，故取其周圍4點(i+1，k)，(i-1，k)，(i，k+1)，(i，k-1)中已淹水的諸點潮位的平均值。

ζi，k=

漲潮時，ζi，k>0，當hi，k≥hmin時，認為該點淹沒，參與計算，反之則不參與計算。

3.2 網(wǎng)格的劃分

河口地區(qū)的特點是水淺、岸線形狀和水下地形復雜、存在各種形狀的沙洲、灘涂和人工建筑物(導流堤等)，導致該地區(qū)的潮流場通常十分復雜。就網(wǎng)格劃分形式而言，不規(guī)則三角形網(wǎng)格非常實用，當前被廣泛采用。不規(guī)則三角形網(wǎng)格有很多優(yōu)點：1)可以很好地擬合水下地形和固邊界形狀，岸線和水下地形概化好；2)網(wǎng)格疏密自由控制，網(wǎng)格布設隨意；3)通用性好，矩形網(wǎng)格可以看作是由若干個直角三角形組成的三角形網(wǎng)格，是不規(guī)則三角形網(wǎng)格的一種特殊形式。

閩江下游河段屬于感潮河段，長期受徑流和潮流的雙重影響，致使河流流態(tài)和形態(tài)復雜，并且地形錯綜復雜、區(qū)域邊界曲折多變，因此對該區(qū)域流態(tài)和水深的模擬存在一定的困難。根據(jù)上述地形變化的不規(guī)則特點和實際需要，模型采用靈活多變的不規(guī)則三角形單元計算網(wǎng)格，對重點區(qū)域進行網(wǎng)格加密，采用整體和局部相結合的模型建模方法建立二維水流數(shù)學模型。通過在河道的交匯口和彎道的地方適當加密，模型能較好地反映實際河流的地形地貌。同時為了提高工作效率，便于修改，編制了計算機單元自動剖分程序，實現(xiàn)了計算區(qū)域網(wǎng)格單元剖分自動化。

4 模型的調(diào)試和驗證

4.1 資料選取

收集到的資料包括：

(1)2019年4月與9月實測的1∶1000和1∶5000閩江干流淮安至亭江段河道地形圖；

(2)2013年12月和2018年2月南港至亭江段水文實測資料，包含潮位、流速資料；

(3)南港近遠期規(guī)劃資料；

(4)南港現(xiàn)有堤防情況。

模型研究區(qū)域為閩江下游淮安至河口亭江處，上游邊界選取南港淮安斷面、南港中段大樟溪江口特大橋斷面和馬尾附近北港壁頭村斷面，下游邊界選取河口亭江斷面。

4.2 邊界條件及網(wǎng)格劃分

采用2017年12月11日～12日的實測潮位作為模型區(qū)域進、出口斷面的邊界控制條件。

模型區(qū)域采用三角形網(wǎng)格進行剖分，本模型共有6461個節(jié)點，12262個三角形網(wǎng)格單元，網(wǎng)格步長控制在200～700m之間，模型區(qū)域的網(wǎng)格圖如圖2～圖4所示。

圖2 模型網(wǎng)格剖分圖Fig.2 Model mesh subdivision

圖3 大樟溪匯入處的網(wǎng)格剖分圖Fig.3 Mesh subdivision map of the confluence of Da Zhang River

圖4 南北港匯流處的網(wǎng)格剖分圖Fig.4 Mesh subdivision of the confluence of the North and South ports

4.3 參數(shù)的選取

為了計算的穩(wěn)定性要求，水流程序計算過程中時間步長必須滿足:

式中：ΔLmin為計算區(qū)域內(nèi)最小的三角形邊長；Hmax為計算區(qū)域內(nèi)最大的水深。

4.4 糙率系數(shù)

糙率的精度直接影響到模型的計算精度[9]。由于影響河道糙率的因素較為復雜，文中采用試算法對其進行調(diào)試。若模型計算誤差過大，則不斷調(diào)整糙率值，直至將誤差控制在一定的范圍內(nèi)。最后確定合理糙率值的范圍為0.025～0.045之間。

4.5 富裕水深的確定

本河段內(nèi)灘地寬闊，且受潮汐和徑流的雙重影響，水流漲落起伏大，隨著水流的漲落，水邊界也會隨之移動，因此，如何合理的模擬水流邊界的移動，直接關系到模擬計算的成敗[10]。模型中采取水位判別法來處理這類動邊界問題，并經(jīng)過調(diào)試取富裕水深為0.05m。

4.6 模型驗證

模型采用2017年12月11日～12日的實測潮位、流速進行驗證，計算值和實測值比較如圖5～圖9所示。由圖可見計算的潮位、流速與實測值吻合情況良好，量值大小上雖有些波動，但基本相當。

圖5 洪塘大橋處潮位過程驗證Fig.5 Verification of tide level process at Hongtang Bridge

圖6 浦上大橋處潮位過程驗證Fig.6 Verification of tide level process at Pushang Bridge

圖7 灣邊斷面處潮位過程驗證Fig.7 Verification of tidal level process at bay side section

圖8 科貢斷面處測點流速大小驗證圖Fig.8 Verification diagram of flow velocity at measuring points at Kegong section

圖9 灣邊斷面處測點流速大小驗證圖Fig.9 Verification diagram of velocity at measuring points at bay side section

圖10～圖13為計算區(qū)域內(nèi)漲急、落急流場平面分布圖，從圖中可以看出流場平面分布和灘槽分布相一致，在大樟溪匯入口、馬尾匯流口附近，水流流態(tài)較為平順，能較好地復現(xiàn)水流形態(tài)。以上分析表明該模型的參數(shù)選取是合理的。

圖10 河道落潮時的大樟溪與閩江交匯處流場圖Fig.10 Flow field diagram of the confluence of Dazhang River and Minjiang River at ebb tide

圖11 河道落潮時的馬尾匯流處流場圖Fig.11 Flow field diagram of Mawei confluence at ebb tide

圖12 河道漲潮時的大樟溪與閩江交匯處流場圖Fig.12 Flow field diagram of the confluence of Dazhang River and Minjiang River at flood tide

圖13 河道漲潮時的馬尾匯流處流場圖Fig.13 Flow field diagram of Mawei Dazhang confluence at flood tide

5 誤差分析

表1是河道中幾個主要斷面實測水位平均值與計算水位平均值的比較及相對誤差分析，通過該表可以看出模型精度較高，誤差均在4.55%以內(nèi)。引起誤差的原因可能是河道匯流、河床形態(tài)復雜等。

表1 主要斷面實測水位與計算水位的比較Tab.1 Comparison between measured water level and calculated water level of main section

6 結語

閩江下游長期受徑流和潮流的雙重影響，特別是文中研究的南港河段，灘涂、沙洲發(fā)育良好，河流流態(tài)復雜多變。文中針對天然河道岸線曲折多變、地形復雜的特點，采用靈活多變的不規(guī)則三角形單元計算網(wǎng)格，應用有限元法建立起二維水流數(shù)學模型，并在模型中加入了動邊界處理技術，使結果更符合實際情況。根據(jù)實測資料對二維水流模型進行調(diào)試、驗證，結果表明計算值與實測值吻合良好，符合該區(qū)域潮流的運動特性。該研究成果為后期進行泥沙數(shù)值模擬和航道整治工程數(shù)值模型研究提供了良好的基礎。