王亞南，郭維東，李東斌，王 穎，錢 彤

(1.沈陽農(nóng)業(yè)大學(xué) 水利學(xué)院，沈陽 110866；2.河北省水利水電勘測設(shè)計院，天津 300250；3.阜新市凌河保護區(qū)管理局，遼寧 阜新 123000；4.朝陽市凌河保護區(qū)管理局，遼寧 朝陽 122000)

我國海岸線長達1.8萬多km，河口及海岸交通便利、資源豐富[1]。河口港的擴建、河口三角洲的整治、河口建閘工程以及漁業(yè)和石油的開發(fā)等，均會導(dǎo)致嚴重的河口泥沙問題，因此對河口地區(qū)泥沙運動規(guī)律的研究具有重要的作用。國內(nèi)外對河口泥沙數(shù)學(xué)模型的研究經(jīng)歷了一維、二維到三維的過程。鄭國棟[2]建立了橫向、縱向、垂向3個維向，分別分析了珠江三角洲的地貌變化, 張心鳳[3]建立二維動邊界非恒定潮流泥沙數(shù)學(xué)模型，對潮流泥沙沖淤情況進行了研究。三維泥沙數(shù)學(xué)模型由于模型結(jié)構(gòu)復(fù)雜、計算工作量較大的因素發(fā)展較緩慢[4]。隨著計算機技術(shù)的發(fā)展以及實際工程和研究的需要，簡單的三維泥沙數(shù)學(xué)模型逐漸應(yīng)用到計算中[5-8]。

根據(jù)河口平面形態(tài)將河口分為喇叭形、彎曲形、順直形、三角洲形[9]。三角洲形河口平面形態(tài)較喇叭口型相似，是由順直形和喇叭口形轉(zhuǎn)化過而來，其最大特點就是在河口出口處形成河口三角洲。本文擬采用數(shù)值模擬方法，對三角洲型感潮河口在不同初始條件和邊界條件下進行數(shù)值模擬。通過地形資料實現(xiàn)所研究河口地區(qū)的建模；選擇數(shù)學(xué)模型，設(shè)定初始條件和邊界條件，并采用實測資料驗證模型的準確性，對三角洲型河口泥沙輸移的規(guī)律進行了計算和分析。

1 基本控制方程

控制方程主要包括：連續(xù)方程、動量方程、輸運方程、泥沙控制方程、地質(zhì)演化方程。

(1)連續(xù)性方程：

(1)

(2)動量方程：

x方向上：

(2)

y方向上：

(3)

(3)輸運方程如下：

(4)

(4)泥沙控制方程。

輸沙公式如下：

qt=qb+qs

(5)

式中：qt為總輸沙率；qs為懸移質(zhì)輸沙率；qb為推移質(zhì)輸沙率。

(6)

式中：T為時間；u為瞬時速度；c為瞬時懸沙濃度；D為瞬時水深；d為懸沙粒徑。

(10)

式中：θc為泥沙啟動的臨界剪切應(yīng)力；θ′為與底表面摩阻有關(guān)的無量綱剪切應(yīng)力；p為一層內(nèi)所有砂礫都啟動的概率；s為泥沙的相對密度；d為泥沙粒徑；ψ為橫向長度尺度；Uf為一個波浪周期內(nèi)的摩擦系數(shù)；μs為靜摩擦系數(shù)(μs=tanφs，φs為休止角)；θc,0為臨界Shield系數(shù)；β為動摩阻參數(shù)。

(5)地形演化控制方程。河床變化率的計算采用泥沙連續(xù)方程：

(11)

式中：n為河床孔隙率；t為時間；Sx為x方向上推移質(zhì)的總輸沙率；Sy為y方向上推移質(zhì)的輸沙率；z為河床高程；x,y為水平笛卡爾坐標；ΔS為泥沙源匯項，對于平衡輸沙，源匯項為0。

為研究地形的變化引入了地貌演化的計算公式：

(12)

式中：ΔtHD為時間步長；Znew為變化后的水位；Zold為變化前的水位。

2 研究方案及模型建立

2.1 三角洲型河口簡介

三角洲型河口是在河口出口處形成河口三角洲。例如長江三角洲和珠江三角洲。隨著河口外延及口門的淤堵，上游處水位不斷被抬高，當水位高出兩岸地面一定高度時，水流就會改道，從自然堤的薄弱處決口，開辟一條新通道，并在新的口門重演上述過程，在這個不斷改道的過程中，形成了河口三角洲，三角洲型河口示意圖如圖1所示。

圖1 三角洲型河口示意圖 Fig.1 Sketch map of delta type Estuary

本次模型計算區(qū)域為河北省秦皇島市洋河河口，采用2009年版的電子海圖及地質(zhì)勘查四隊提供的2011年6月實測地形圖，為方便計算，本文模擬典型三角洲型河口，河床縱比降為0。模型中水域面積合計14.39 km2，其中海域面積13.00 km2。感潮河口段長度2.2 km、河口平均寬度850 m。計算區(qū)域網(wǎng)格劃分結(jié)果如圖2所示。

圖2 三角洲網(wǎng)格劃分結(jié)果圖Fig.2 The Result of Mesh Genetation in Yang River Estuary

2.2 研究方案

本次研究方案主要對不同流量、不同來沙中值粒徑條件下的河口水位、流速及河底高程變化進行了模擬。方案如表1。

2.3 系數(shù)選取

水平渦黏系數(shù)估計公式如下：

表1 三角洲型河口模擬方案Tab.1 The Calculation scheme table of delta style

(13)

式中：E為渦黏系數(shù)；A為單元面積；u和v為流速在x和y方向的流速；C為系數(shù)，范圍在0.094～0.3之間。

糙率系數(shù)的選取根據(jù)現(xiàn)有的實測資料以及大量實驗總結(jié)出的經(jīng)驗值。本文選取的糙率值為0.032(m·s-1)。

2.4 模型驗證

模型率定的計算區(qū)域上邊界定在洋河大橋以下的橡膠壩下游斷面，距洋河大橋60 m，河口寬度540 m，河口長度1.2 km。海域邊界取在距河口1.9 km邊界，水域總覆蓋面積約6.25 km2，洋河河口水系圖如圖3所示。驗證指標選為水位和流速，潮位數(shù)據(jù)及泥沙數(shù)據(jù)均來自秦皇島礦產(chǎn)水文工程地質(zhì)大隊勘察報告，在模擬時段內(nèi)，實測數(shù)據(jù)和模擬數(shù)據(jù)計算結(jié)果如圖4、圖5。

圖3 洋河河口水系圖Fig.3 The Yanghe River River Estuary map

圖4 數(shù)值模擬與實測數(shù)據(jù)比較(水位)Fig.4 The contrast between numerical simulation and measured data(Water level)

圖5 數(shù)值模擬與實測數(shù)據(jù)比較(流速)Fig.5 The Contrast Between Numerical Simulation and Measured Data(Current Speed)

為分析數(shù)值模擬的精度，分別計算其平均絕對誤差(MAE)和均方根(RMSE)。計算結(jié)果見表2。

表2 水位與流速數(shù)值模擬評價效果Tab 2 Evaluation effect of numerical simulation：water level and current speed

以上計算精度均滿足《海岸與河口潮流泥沙模擬技術(shù)規(guī)程》(JTS/T231-2-2010)中的精度要求。因此該模型在水動力模擬方面的結(jié)果是可信的。

3 計算結(jié)果分析

3.1 方案一結(jié)果分析

在相同控制流量Q1=500 m3/s、不同來沙中值粒徑d50=0.1 mm、d50=0.2 mm及d50=0.3 mm條件下，計算A、B、C點的水位、流速及河床高程的變化趨勢。

3.1.1 水位流速數(shù)值模擬

由圖6、圖7、圖8可知，d50=0.1 mm工況的水位變化在模擬時段的前段時間波動較劇烈，在A、B點上，5 d后水位逐漸趨于平穩(wěn)，A點水位范圍為2.2～2.5 m，B點水位范圍為1.6～2.4 m，在C點上，10 d后逐漸趨于平穩(wěn)，水位范圍為1.2～2.2 m；流速的變化趨勢為：在A、B上波動較平穩(wěn)，A點變化范圍為0.6～0.7 m/s，B點變化范圍為0.4～0.55 m/s；在C點，模擬時候的后半段流速明顯增大，變化范圍為0.2～0.6 m/s。流速波動趨勢與水位波動趨勢呈反向關(guān)系。d50=0.2 mm工況和d50=0.3 mm與d50=0.1 mm工況的水位變化趨勢不一致，在同一時刻，d50=0.1 mm工況的瞬時流速大于其余兩種工況。

圖6 方案一A點水位、流速曲線Fig.6 The curve of A point water level and current speed in plan 1

圖7 方案一B點水位、流速曲線Fig.7 The curve of B point water level and current speed in plan 1

圖8 方案一C點下水位、流速曲線Fig.8 The curve of C point water level and current speed in plan 1

3.1.2 河床高程變化數(shù)值模擬

在A點，3種中值粒徑工況下的河床高程呈現(xiàn)不同的特點，d50=0.1 mm工況曲線中可以看出在開始的1天內(nèi)，河床高程增長迅速，6月16日之后保持平穩(wěn)，在6月26日之后河床高程再次增長；d50=0.2 mm工況則以同一個河床增長率增長，波動不大。d50=0.3 mm工況下，河床高程變化較小，略有淤積，在6月26日之后略有下降。B點位于三角洲洲頭，在三種中值粒徑下保持河床高程增長，d50=0.1 mm工況河床變化最為劇烈，d50=0.3 mm工況則最為平緩。C點位于河口出口三角洲之后，其沖淤規(guī)律與A、B點不同，中值粒徑為d50=0.2 mm及d50=0.3 mm工況存在沖刷，且d50=0.3 mm沖刷較嚴重。見圖9-圖11。

圖9 方案一A點河床高程曲線Fig.9 The curve of A point bed level in plan 1

圖10 方案一B點河床高程曲線Fig. 10 The curve of B point bed level in plan 1

圖11 方案一C點河床高程曲線Fig.11 The curve of C point bed level in plan 1 注：圖9-11中，橫軸為模擬時間，左側(cè)縱軸為Q=500 m3/s，d50=0.1mm工況下的河床高程；右側(cè)縱軸則為Q=500 m3/s，d50=0.2 mm及d50=0.3 mm工況下的河床高程。

3.1.3 淤積部位數(shù)值模擬

在研究三角洲形河口的淤積部位的過程中，除了提取典型斷面單點處的水力及泥沙要素外，還提取了平面二維的河床高程圖，進而可以從整體上觀察不同工況下的沖淤部位的規(guī)律。

圖12-圖14表示是Q=500 m3/s，d50=0.1 mm、d50=0.2 mm、d50=0.3 mm工況下，模擬進行了15 d以后的結(jié)果。從圖12-圖14可知，在三角洲型河口存在比較明顯的沖淤區(qū)域，當上游來水經(jīng)過河口三角洲時，三角洲洲頭形成低流速區(qū)，這個區(qū)域中的水流掀沙能力很低，形成淤積區(qū)，Q=500 m3/s，d50=0.1 mm工況下淤積較嚴重，最大淤積厚度為1.2 m。河道斷面分汊處造成沖刷，Q=500 m3/s，d50=0.2 mm工況下沖刷較劇烈，最大沖刷深度為0.06 m。水流從分汊河道流出后在三角洲尾部形成滯留區(qū)，流速降低產(chǎn)生淤積，Q=500 m3/s,d50=0.1 mm工況下，最大淤積厚度為1.2 m。

圖12 Q=500 m3/s, d50=0.1 mm計算15 d后河底高程等值線圖Fig.12 The curve of bed level isolines in the time 15 days with the model in Q=500 m3/s, d50=0.1 mm

圖13 Q=500 m3/s, d50=0.2 mm計算15 d后河底高程等值線圖Fig.13 The curve of bed level isolines in the time 15 days with the model in Q=500 m3/s, d50=0.2 mm

圖14 Q=500 m3/s, d50=0.3 mm計算15 d后河底高程等值線圖Fig.14 The curve of bed level isolines in the time 15 days with the model in Q=500 m3/s, d50=0.3 mm

3.2 方案二結(jié)果分析

在相同來沙中值粒徑d50=0.16 mm，不同控制流量Q2=200 m3/s、Q3=10 m3/s條件下，模擬計算A、B、C點的水位、流速及河床高程的變化趨勢。

3.2.1 水位流速數(shù)值模擬

圖15表明，兩種工況在流量不同、中值粒徑不同情況下，同一時刻水位相差較小，整體變化的趨勢與潮水位的變化趨勢相同。

圖15 方案二A、B、C點水位曲線Fig.15 The curve of ABC point water level in plan 2

從圖16與圖17可知，在三角洲型河口中流速變化較為有規(guī)律，接近正弦曲線。在同一時刻A點波動范圍為0.14～0.34 m/s，B點波動范圍為0.04～0.22 m/s，C點波動范圍為0～0.14 m/s，A點流速最大，C點流速最小。由于流量的差異，Q=200 m3/s工況整體波動范圍比Q=10 m3/s工況流速波動范圍大。在Q=10 m3/s工況下，A、B、C流速曲線較為接近，與Q=200 m3/s工況相比較，波動沒有明顯規(guī)律。在同一時刻，B點流速最大，C點流速最小。

圖16 Q2=200 m3/s, d50=0.16 mm ABC點流速曲線Fig.16 The curve of A、B、C point Current Speed in Q2=200 m3/s

圖17 Q2=10 m3/s, d50=0.16 mm ABC點流速Fig.17 The curve of A、B、C point Current Speed in Q2=10 m3/s

3.2.2 河床高程變化數(shù)值模擬

為研究單點處的河床高程變化，提取了A、B、C三點在整個模擬時段內(nèi)的河床高程值，圖18表示的是在以上兩種工況下三點的河床高程與時間的變化曲線。

圖18 方案二ABC點河床高程Fig.18 The curve of ABC point bed level in plan 2

由圖18可知，在中值粒徑相同情況下，流量、位置對沖淤規(guī)律有一定的影響。在Q=200 m3/s工況下A、B兩點為淤積區(qū)域，且B點淤積程度較大，C點高程為負值，為沖刷區(qū)域。在Q=10 m3/s工況下，流量較前一種工況小，在三角洲挾沙量小，淤積程度較小。分汊流帶來的泥沙在三角洲后滯留區(qū)的補給量較小，導(dǎo)致在枯水期C點的沖刷程度加大。

3.2.3 淤積部位數(shù)值模擬

為研究河口沖刷范圍及規(guī)律，提取了在模擬開始15 d后對應(yīng)的河口河床高程等值線圖。圖19與圖20表示Q=200 m3/s及10 m3/s，d50=0.16 mm工況下的兩個時間上的平面二維河床高程等值線圖。

圖19 Q=200 m3/s, d50=0.16 mm計算15 d后河底高程等值線圖Fig.19 The curve of bed level isolines in the time 15 days in Q=200 m3/s

圖20 Q=10 m3/s,d50=0.16mm計算15 d后河底高程等值線圖Fig.20 The curve of bed level isolines in the time 15 days in Q=10 m3/s

從圖19、圖20可知，Q=200 m3/s與Q=10 m3/s工況下的沖淤結(jié)果是不同的，在中水期(Q=200 m3/s)河道是在徑流和潮流共同作用下，在分汊河道和三角洲尾部正中形成沖刷區(qū)域。在洲頭和三角洲尾部兩端形成淤積區(qū)域。在枯水期潮流作用大于徑流占主導(dǎo)地位，進水口至洲頭區(qū)域河床高程存在很小的沖刷，沖刷深度在三角洲至外海區(qū)域受海流作用明顯。

4 結(jié) 語

(1)豐水期不同泥沙中值粒徑對河口水位的影響較小，潮水位對河口水位影響較大，起主導(dǎo)作用，但是河口三角洲的存在使得不同點的水位呈現(xiàn)差異。隨著時間的變化，波動趨于穩(wěn)定。在對應(yīng)的流速曲線和水位曲線成反向關(guān)系。三角洲型河口也存在比較明顯的沖淤區(qū)域，三角洲兩側(cè)形成沖刷區(qū)，三角洲洲頭以及三角洲尾部形成淤積區(qū)。

(2)中水期和枯水期，河口水位的變化規(guī)律和豐水期是一致的。在中水期，流速變化較為有規(guī)律，接近正弦曲線。在同一時刻進口的瞬時流速最大，三角洲后點處流速最小。在枯水期，流速曲線波動范圍較大，同一時刻三角洲頭處瞬時流速最大，進水口瞬時流速最小。中值粒徑相同時，流量及點位對沖淤規(guī)律有一定的影響。在中水期，三角洲兩側(cè)和三角洲尾部正中形成沖刷區(qū)域，在洲頭和三角洲尾部兩端形成淤積區(qū)域；在枯水期，進水口至洲頭區(qū)域河床高程存在很小的沖刷，在三角洲至外海區(qū)域受海流作用明顯形成沖刷區(qū)域，出口兩端是淤積區(qū)。

□

[1] 李文丹. 復(fù)雜河口海岸地區(qū)水動力數(shù)值模擬研究[D]. 哈爾濱：哈爾濱工程大學(xué),2008.

[2] 鄭國棟,顧立忠,李虎成,等. 珠江三角洲河道地貌變化對網(wǎng)河水情影響研究[J]. 中國農(nóng)村水利水電,2010,(7):33-36.

[3] 張心鳳,李 越. 銀湖灣A區(qū)圍墾規(guī)劃潮流泥沙數(shù)值模擬研究[J]. 中國農(nóng)村水利水電,2014,(8):62-66,72.

[4] 邵宇陽. 近岸潮流和泥沙運動三維數(shù)值模擬[D]. 南京：河海大學(xué),2005.

[5] 郭維東,田茹妍,李東斌,劉冰. 彎曲型感潮河口泥沙輸移規(guī)律的數(shù)值模擬[J]. 人民黃河,2016,(4):1-6.

[6] 趙張益. 河口海岸三維水沙運動的間斷有限元模型研究[D]. 天津：天津大學(xué),2014.

[7] 王效遠. 考慮波浪破碎影響的近岸三維泥沙數(shù)學(xué)模型[D]. 天津：天津大學(xué),2009.

[8] 王崇浩,韋永康. 三維水動力泥沙輸移模型及其在珠江口的應(yīng)用[J]. 中國水利水電科學(xué)研究院學(xué)報,2006,(4):246-252.

[9] 熊紹隆,曾 劍. 潮汐河口分類指標與河床演變特征研究[J]. 水利學(xué)報,2008,(12):1 286-1 295.