三角洲型感潮河口泥沙輸移規(guī)律研究

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

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

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

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

1 基本控制方程

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

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

(1)

(2)動(dòng)量方程：

x方向上：

(2)

y方向上：

(3)

(3)輸運(yùn)方程如下：

(4)

(4)泥沙控制方程。

輸沙公式如下：

qt=qb+qs

(5)

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

(6)

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

(10)

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

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

(11)

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

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

(12)

式中：ΔtHD為時(shí)間步長(zhǎng)；Znew為變化后的水位；Zold為變化前的水位。

2 研究方案及模型建立

2.1 三角洲型河口簡(jiǎn)介

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

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

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

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

2.2 研究方案

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

2.3 系數(shù)選取

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

表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)有的實(shí)測(cè)資料以及大量實(shí)驗(yàn)總結(jié)出的經(jīng)驗(yàn)值。本文選取的糙率值為0.032(m·s-1)。

2.4 模型驗(yàn)證

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

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

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

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

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

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

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

3 計(jì)算結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

圖11 方案一C點(diǎn)河床高程曲線Fig.11 The curve of C point bed level in plan 1 注：圖9-11中，橫軸為模擬時(shí)間，左側(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ò)程中，除了提取典型斷面單點(diǎn)處的水力及泥沙要素外，還提取了平面二維的河床高程圖，進(jìn)而可以從整體上觀察不同工況下的沖淤部位的規(guī)律。

圖12-圖14表示是Q=500 m3/s，d50=0.1 mm、d50=0.2 mm、d50=0.3 mm工況下，模擬進(jìn)行了15 d以后的結(jié)果。從圖12-圖14可知，在三角洲型河口存在比較明顯的沖淤區(qū)域，當(dāng)上游來(lái)水經(jīng)過(guò)河口三角洲時(shí)，三角洲洲頭形成低流速區(qū)，這個(gè)區(qū)域中的水流掀沙能力很低，形成淤積區(qū)，Q=500 m3/s，d50=0.1 mm工況下淤積較嚴(yán)重，最大淤積厚度為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計(jì)算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計(jì)算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計(jì)算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é)果分析

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

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

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

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

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

圖16 Q2=200 m3/s, d50=0.16 mm ABC點(diǎn)流速曲線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點(diǎn)流速Fig.17 The curve of A、B、C point Current Speed in Q2=10 m3/s

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

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

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

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

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

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

圖19 Q=200 m3/s, d50=0.16 mm計(jì)算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計(jì)算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ū)域。在洲頭和三角洲尾部?jī)啥诵纬捎俜e區(qū)域。在枯水期潮流作用大于徑流占主導(dǎo)地位，進(jìn)水口至洲頭區(qū)域河床高程存在很小的沖刷，沖刷深度在三角洲至外海區(qū)域受海流作用明顯。

4 結(jié) 語(yǔ)

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

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

□

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

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

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

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

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

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

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

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

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