(浙江省水利水電勘測設(shè)計院，杭州 310002)

1 研究背景

浙江沿海地區(qū)在灘涂保護與利用開發(fā)中，海堤龍口段的地基大多數(shù)為深厚淤泥軟土地基，受進排水要求限制，海堤龍口段拋填石料厚度一般較薄，其預(yù)壓固結(jié)程度不如一般海堤段的預(yù)壓固結(jié)程度，其抗沖刷能力相對較差。龍口寬度太窄則龍口度汛期漲落潮時水流相對較急，容易造成龍口嚴(yán)重沖刷，影響工程安全穩(wěn)定性；龍口寬度太寬則造成堵口難度大，影響工程總體工期，且不易合龍。因此，選擇合適的龍口寬度在深厚淤泥地基海堤工程中顯得尤為重要。

吳繼偉[1]采用水量平衡法編制了海堤龍口水力計算軟件，但其未考慮龍口在堵口過程中底坎、口門寬度的動態(tài)變化過程，且對堵口過程中拋石底坎的透水性以及水閘的影響等要素未深入考慮，軟件計算模型仍有待深入細化。

趙庚潤等[2]利用流體有限元軟件Fluent中的k-ε湍流模型對海堤龍口二維流動進行了數(shù)值模擬計算，認為落潮時龍口最大流速可達8.70 m/s，漲潮時龍口最大流速可達10.60 m/s；其中大流速區(qū)一般位于堤頂或背流側(cè)的區(qū)域；甚至龍口附近會出現(xiàn)水躍或者卷浪型流態(tài)。

胡志根等[3]采用物理模型試驗對大比降、非對稱基床中的合龍進行研究，觀測分析堵口合龍進占方式、龍口位置、龍口水流流態(tài)及其龍口水力參數(shù)之間的關(guān)系。

楊磊等[4]通過相關(guān)模型試驗對寬戧堤龍口堵口范圍內(nèi)的水流特性進行試驗研究，認為當(dāng)寬戧堤堤頂寬度B與總水頭H0的比值B/H0在311～517之間時，龍口內(nèi)水面形態(tài)將逐漸由單降型流態(tài)漸變?yōu)殡p降型流態(tài)；但當(dāng)B/H0>517時,龍口戧堤堤頂?shù)膶挾刃?yīng)逐漸減弱。

工程實際海堤龍口計算中外海潮位等水力邊界條件不斷變化，導(dǎo)致龍口水力計算較復(fù)雜，一般需采用簡化小程序計算(如文獻[1])?；诖耍疚囊罁?jù)水力學(xué)理論，結(jié)合MATLAB軟件，采用黃朝煊等[5]研究中給出Runge-Kutta法對海堤龍口水力計算進行分析研究，以便工程應(yīng)用推廣。

2 龍口水力計算基本方程

根據(jù)文獻[6]和文獻[7]可知，由圍區(qū)內(nèi)進(出)水量平衡原理得龍口水力計算基本方程為

(1)

q外=q閘+q泄+q滲，

(2)

(3)

(4)

式中：q內(nèi)為內(nèi)陸流域來水平均流量(m3/s)；q外為受外海潮位影響流進(出)圍區(qū)的流量(m3/s)；q閘為水閘過流流量(m3/s)；q泄為口門堰上溢流流量(m3/s)；q滲為計算時段內(nèi)海堤龍口拋石體滲流平均流量(m3/s)；V[Z(t)]為t時刻圍區(qū)內(nèi)庫容量(m3)；Z(t)為t時刻圍區(qū)內(nèi)水位(m)；h(t)為t時刻外海側(cè)潮位(m)；F[Z(t)]為圍區(qū)水位為Z(t)時的水面面積；mμ為龍口水力計算的流量系數(shù)；Bω為實際等效過流計算寬度；k為參數(shù)，堰流時k=1.5。

由于外海潮位、圍區(qū)庫容等邊界條件為非線性性關(guān)系，因此，以上龍口水力計算方程為非線性常微分方程，目前還難以給出其解析解。

3 基于四階Runge-Kutta法的龍口水力計算

根據(jù)以上方程式(1)—式(4)，代換后可得水力計算控制方程為

(5)

記

則式(5)可簡化為

(6)

可采用經(jīng)典的四階Runge-Kutta法求解控制方程得

(7)

其中：

K1i=fti,Zti；

K4i=fti+Δt,Zti+ΔtK3i。

式中： Δt為計算時間段的步長；Zti為在ti=iΔt時刻的圍區(qū)水位；K1i,K2i,K3i,K4i為i時段的積分參數(shù)，由Z(ti)，ti等計算確定。

根據(jù)以上數(shù)值計算法，可由外海側(cè)起調(diào)潮位、圍區(qū)側(cè)起調(diào)水位Zt0等邊界條件計算求解出圍區(qū)水位隨時間的過程線Z(t)，其中Runge-Kutta法的求解流程如圖1所示，可利用MATLAB[8-9]軟件數(shù)值編程計算實現(xiàn)，其中計算時間步長Δt一般可取0.1～0.5 h。

圖1 Runge-Kutta法求解過程示意圖Fig.1 Sketch of the solution process of Runge-Kuttamethod

其中深厚淤泥地基上的海堤龍口典型斷面如圖2所示。

圖2 深厚淤泥地基上的海堤龍口典型斷面示意圖Fig.2 Schematic diagram of the closure gap section of sea dyke on deep silt foundation

4 龍口最大水位差的解析分析

由于外海潮位的近似周期性變化特性，假設(shè)外海側(cè)潮位過程線為

(8)

其中，

式中：h為潮位；T0為潮位半周期，即潮位周期為2T0；H為最大潮差;χ為平均潮位。

龍口、水閘等總過流能力的計算公式為

(9)

在忽略內(nèi)陸來水情況下，依據(jù)式(1)流量平衡方程得：

(10)

(11)

其中圍區(qū)側(cè)起調(diào)水位取為Z0，參數(shù)Z0為海堤龍口底檻高程；值得說明的是，在落潮時需對式(11)取絕對值。

將式(11)代入式(10)可得

(12)

按照工程計算經(jīng)驗，在最不利水力條件即圍區(qū)側(cè)水位與外海側(cè)潮位之間最大水位差情況(如圖3所示)[11]一般出現(xiàn)在高潮位之后，即滿足dH0/dt=0時，對式(11)兩邊求導(dǎo)，變化可得

(13)

圖3 潮位與圍區(qū)內(nèi)水位關(guān)系曲線示意圖Fig.3 Relation between tide level and water level in the closure area

在龍口水力計算中，主要水力參數(shù)是圍區(qū)水位與外海潮位間水位差最大值H0max，其中根據(jù)浙江省水利水電勘測設(shè)計院龍口水力計算經(jīng)驗，最大水位差出現(xiàn)的時間一般為tm=t0+(1.0～1.5)T0，即低潮位之前(0～0.5)T0。記龍口內(nèi)外兩側(cè)水位差最大時刻為tm，此時最大水位差為H0max，相應(yīng)圍區(qū)水位為Z(tm)。

將式(13)代入式(12)，則得

(14)

由于式(14)具有2個待求變量tm和Ztm，還需要1個方程才能求解，根據(jù)水量平衡建立方程，即

(15)

式中參數(shù)K為圍區(qū)庫容修正系數(shù)，一般取值為1.0～1.2。

由式(14)和式(15)，通過數(shù)學(xué)變換，可得以下非線性代數(shù)方程組，即

根據(jù)方程組(16)，求解非線性方程組可解出待求參數(shù)tm和Z(tm)。

(17)

依據(jù)浙江省水利水電勘測設(shè)計院龍口水力計算經(jīng)驗，一般μ=1.0～1.50，借助MATLAB軟件計算求解出μ和tm，并根據(jù)式(8)和式(11)給出最大水位差H0max計算式。

設(shè)圍區(qū)側(cè)的庫容曲線滿足二次拋物線，即

V[Zt]=a2Z2t+a1Zt+a0。(18)

對式(18)求導(dǎo)可知

F[Zt]=2a2Zt+a1。

(19)

式中常參數(shù)a0,a1,a2根據(jù)實測庫容關(guān)系曲線擬合得出，將V[Z(t)]以及F[Z(t)]代入以上非線性代數(shù)方程(17)，得到最大水位差時龍口圍區(qū)側(cè)水位與龍口寬龍口底坎高程之間的關(guān)系函數(shù)。

外海潮位與圍區(qū)水位之間的最大水位差為

(20)

進而依據(jù)《灘涂治理工程規(guī)范》[6]和《海堤工程設(shè)計規(guī)范》[7]可知最大流速為

(21)

根據(jù)式(17)、式(20)及式(21)，可對龍口水力特性進行計算分析，并能繪制漲潮、落潮時的圍區(qū)潮位最大水位差與龍口寬度關(guān)系曲線，進而能給出龍口段漲潮、落潮時的最大水流流速等水力特性，為工程實際中深厚淤泥地基海堤龍口抗沖刷設(shè)計提供計算依據(jù)。

5 算例分析

5.1 項目概況

某圍墾工程圍墾面積2 296.6 hm2，工程主要由海堤、水閘、施工道路等組成。根據(jù)防護對象的規(guī)模和重要性，確定工程防潮(洪)標(biāo)準(zhǔn)為50 a一遇。該工程海堤級別為3級，水閘級別同為3級，施工道路、河道等次要建筑物為4級，圍堰等臨時建筑物為5級。海堤地基土層主要為Ⅲ0層淤泥、Ⅲ1層淤泥質(zhì)粉質(zhì)黏土、ⅢSis粉砂、粉土混淤泥、Ⅲ2層淤泥、Ⅳ1淤泥質(zhì)粉質(zhì)黏土、Ⅳ2層粉質(zhì)黏土、黏土組成。

其中龍口布置原則：①龍口位置宜考慮圍區(qū)涂面較低部位，以利于圍區(qū)潮流進出；②龍口位置應(yīng)離水閘有一定距離，以免影響水閘施工和水閘泄流影響堵口；③龍口型式采用“寬高式龍口”龍口。根據(jù)龍口水力條件，為減小龍口流速防止沖刷和簡化保護措施，目前浙江省大中型圍墾工程和堵港工程的度汛龍口型式均采用“寬高式龍口”，即龍口寬度盡可能寬一些，底檻高程適當(dāng)抬高。

5.2 龍口水力計算

該工程圍堤龍口度汛設(shè)計潮型，以非汛期5 a一遇設(shè)計高潮位3.94 m為控制水位，將典型潮位過程作適當(dāng)修正后采用；圍堤度汛期設(shè)計潮位見表1。

表1 龍口度汛期設(shè)計潮位過程線Table 1 Process of design tide level at closure gap in flood season

其中通過MATLAB軟件對外海設(shè)計潮位進行數(shù)值擬合分析得

h(t)=0.76+2.67sin0.507t+1.407 。 (22)

其中潮位過程線與擬合函數(shù)關(guān)系如圖4所示。

圖4 龍口度汛期設(shè)計潮位過程線及數(shù)值擬合Fig.4 Process and numerical fitting of design tide level at closure gap in flood season

該工程圍堤內(nèi)側(cè)庫容與圍區(qū)水位關(guān)系通過實測地形曲線量算，見表2。

表2 圍區(qū)庫容-水位關(guān)系

通過MATLAB軟件數(shù)值擬合，給出庫容-水位關(guān)系曲線(見圖5)，擬合函數(shù)為

V[Z(t)]=39.2Z2(t)+460Z(t)+453 。

(23)

圖5 圍區(qū)庫容-水位關(guān)系曲線及數(shù)值擬合Fig.5 Measured and fitted relation between storage capacity of closure area and water level

根據(jù)本文方法，對該圍墾工程龍口水力條件進行計算研究，其中漲潮、落潮時龍口最大流速Vmax與龍口底高程Z0、龍口寬度B之間關(guān)系曲線族見圖6。

圖6 龍口最大流速Vmax與龍口底高程、 龍口寬度B之間關(guān)系曲線族Fig.6 Set of curves of maximum flow velocity Vmax at closure gap versus bottom height and width of closure gap

通過以上數(shù)值計算分析可知：當(dāng)?shù)讬懜叱淘?.0以上時，龍口處漲落潮流速隨龍口寬度增大而減小，且流速變化較明顯；當(dāng)?shù)讬懜叱倘?.0 m以上時，龍口處漲落潮流速隨龍口寬度變化相對較小。

根據(jù)龍口布置及控制漲落潮流速等一般要求，同時結(jié)合本工程涂面高程實際情況，以最大控制流速不大于3 m/s，通過以上數(shù)值計算，最終設(shè)計選定龍口底檻高程0.5 m，龍口設(shè)計寬為300 m。

6 結(jié)論與建議

針對深厚淤泥軟土灘涂的海堤龍口抗水流沖刷能力較差的問題，本文對海堤龍口水力學(xué)特性進行了研究分析，主要結(jié)論如下：

(1) 依據(jù)水力學(xué)理論及水量平衡原理建立了海堤龍口水力計算基本方程，該方程為非線性常微分方程，一般情況下難以給出精確解析解，可依據(jù)數(shù)值分析理論利用經(jīng)典四階Runge-Kutta法求解，結(jié)合MATLAB計算軟件數(shù)值編程給出數(shù)值解，其中計算時間步長Δt一般可取0.1～0.5 h。

(2) 利用正弦函數(shù)擬合外海潮位周期變化特性，利用二次函數(shù)擬合圍區(qū)庫容-水位關(guān)系曲線，結(jié)合相關(guān)數(shù)值等效分析法及海堤龍口計算經(jīng)驗，給出了海堤龍口最大水位差與龍口寬度之間的近似解析計算方法，為工程實際應(yīng)用提供方便。

(3)最后通過工程案例，利用MATLAB軟件計算對比分析，表明該案例中底檻高程在0.0以上時，龍口處漲落潮流速隨龍口寬度增大而減小，且流速變化較明顯；當(dāng)?shù)讬懜叱倘?.0 m以上時，龍口處漲落潮流速隨龍口寬度變化相對較小，與工程實際經(jīng)驗一致。

值得說明的是，由于漲落潮時口門上水流流態(tài)較復(fù)雜，涉及高淹沒出流情況，因此本文簡化計算方法難免存在一定不足，筆者將進行后續(xù)深入研究。