朱明栓，張思慧，江中林

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

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

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

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

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

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

2.1 基本方程

連續(xù)性方程:

運(yùn)動(dòng)方程[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)

邊界條件：

對(duì)水邊界Γ1

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

對(duì)岸邊界Γ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法來(lái)求解平面二維河道水流方程[6]，對(duì)控制方程進(jìn)行離散，設(shè)置變量并用Fortran語(yǔ)言編制計(jì)算程序。

3 模型中關(guān)鍵問(wèn)題的處理

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

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

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

在實(shí)際計(jì)算中一般把潮灘區(qū)的動(dòng)邊界按二維問(wèn)題來(lái)處理。對(duì)于平面上任一潮灘點(diǎn)(i，k)，瞬時(shí)水深hi，k=Di，k+ζi，k，其中Di，k，ζi，k分別為(i，k)點(diǎn)的灘地高程及潮位。

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

ζi，k=

漲潮時(shí)，ζi，k>0，當(dāng)hi，k≥hmin時(shí)，認(rèn)為該點(diǎn)淹沒，參與計(jì)算，反之則不參與計(jì)算。

3.2 網(wǎng)格的劃分

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

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

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

4.1 資料選取

收集到的資料包括：

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

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

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

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

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

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

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

模型區(qū)域采用三角形網(wǎng)格進(jìn)行剖分，本模型共有6461個(gè)節(jié)點(diǎn)，12262個(gè)三角形網(wǎng)格單元，網(wǎng)格步長(zhǎ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ù)的選取

為了計(jì)算的穩(wěn)定性要求，水流程序計(jì)算過(guò)程中時(shí)間步長(zhǎng)必須滿足:

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

4.4 糙率系數(shù)

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

4.5 富裕水深的確定

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

4.6 模型驗(yàn)證

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

圖5 洪塘大橋處潮位過(guò)程驗(yàn)證Fig.5 Verification of tide level process at Hongtang Bridge

圖6 浦上大橋處潮位過(guò)程驗(yàn)證Fig.6 Verification of tide level process at Pushang Bridge

圖7 灣邊斷面處潮位過(guò)程驗(yàn)證Fig.7 Verification of tidal level process at bay side section

圖8 科貢斷面處測(cè)點(diǎn)流速大小驗(yàn)證圖Fig.8 Verification diagram of flow velocity at measuring points at Kegong section

圖9 灣邊斷面處測(cè)點(diǎn)流速大小驗(yàn)證圖Fig.9 Verification diagram of velocity at measuring points at bay side section

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

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

圖11 河道落潮時(shí)的馬尾匯流處流場(chǎng)圖Fig.11 Flow field diagram of Mawei confluence at ebb tide

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

圖13 河道漲潮時(shí)的馬尾匯流處流場(chǎng)圖Fig.13 Flow field diagram of Mawei Dazhang confluence at flood tide

5 誤差分析

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

表1 主要斷面實(shí)測(cè)水位與計(jì)算水位的比較Tab.1 Comparison between measured water level and calculated water level of main section

6 結(jié)語(yǔ)

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