何 昆，楊元平，陳 剛，錢盛杰，張芝永

（1.浙江省水利河口研究院（浙江省海洋規(guī)劃設(shè)計(jì)研究院）， 浙江 杭州 310017；2.海寧市鹽倉(cāng)江堤管理所， 浙江 海寧 314400）

1 問題的提出

丁壩作為一種護(hù)岸建筑物與河道整治建筑物，在水利工程中應(yīng)用廣泛。涌潮是一種發(fā)生于喇叭形河口或海灣地區(qū)的特殊潮汐現(xiàn)象，具有潮差大、水動(dòng)力強(qiáng)的特點(diǎn)，易對(duì)丁壩造成沖刷損壞。因此，研究涌潮作用下丁壩的水力特性，對(duì)丁壩的布置優(yōu)化及安全運(yùn)行具有重要意義。

對(duì)丁壩水力特性的研究一直是河流動(dòng)力學(xué)科的熱門研究方向之一，研究方法主要有物理模型試驗(yàn)和數(shù)值模擬等。1953 年Ahmed[1]首次對(duì)丁壩繞流的水面線變化情況進(jìn)行測(cè)量；彭靜[2]等基于顏料示蹤和油膜技術(shù)對(duì)丁壩群附近的流場(chǎng)分布進(jìn)行可視化記錄研究；Jennifer[3]等采用激光多普勒測(cè)速儀（LDV）對(duì)丁壩附近的水位和流速分布進(jìn)行觀測(cè)，并在此基礎(chǔ)上研究丁壩水流的三維紊動(dòng)特性；凌建明[4]等運(yùn)用二維水深平均k-ε 紊流模型，采用有限元的方法對(duì)控制方程進(jìn)行離散，通過計(jì)算得到不同工況下繞流丁壩附近的流場(chǎng)分布；劉超[5]等采用有限體積法離散求解二維淺水方程，輔以RNG k-ε湍流模型，模擬丁壩壩頭附近的湍流流場(chǎng)分布；胡子俊[6]等采用RANS 方法對(duì)丁壩繞流進(jìn)行模擬，將不同湍流模型與試驗(yàn)觀測(cè)結(jié)果進(jìn)行對(duì)比后得出RNG k-ε 湍流模型的精度較高，通過數(shù)值計(jì)算得到的流場(chǎng)變化形態(tài)與試驗(yàn)觀測(cè)和客觀流場(chǎng)運(yùn)動(dòng)規(guī)律吻合的結(jié)論。

但目前對(duì)丁壩的三維數(shù)值模擬或物理模型模擬，大部分針對(duì)恒定水流形式，對(duì)涌潮作用下丁壩水力特性的研究較少。本文通過建立三維數(shù)值水槽模型，對(duì)涌潮作用下丁壩的水力特性規(guī)律進(jìn)行探索，結(jié)果可為涌潮作用地區(qū)的丁壩設(shè)計(jì)提供參考。

2 理論模型簡(jiǎn)介

2.1 控制方程

模擬的水體運(yùn)動(dòng)為不可壓縮黏性流體在笛卡爾坐標(biāo)系中的運(yùn)動(dòng)，其運(yùn)動(dòng)的控制方程為連續(xù)性方程及動(dòng)量守恒方程（Navier-Stokes 方程）。其中，由于采用VOF 算法進(jìn)行水體自由表面追蹤，在連續(xù)性方程和動(dòng)量方程中加入體積和面積分?jǐn)?shù)參數(shù)，具體表達(dá)式如下：

連續(xù)性方程：

動(dòng)量方程：

式中：ρ為流體密度，kg/m3；p為壓強(qiáng)，Pa；VF為流體的體積分?jǐn)?shù)；Ax，Ay，Az是x，y，z3 個(gè)方向流體的面積分?jǐn)?shù)；u，v，w為x，y，z3 個(gè)方向的速度分量，m/s；Gx，Gy，Gz為x，y，z3 個(gè)方向的重力加速度，m/s2；fx，fy，fz為在x，y，z3 個(gè)方向的流體黏滯力加速度，m/s2。

2.2 紊流模型

流體與丁壩作用后通常會(huì)產(chǎn)生變形破碎和湍流漩渦，為更好地模擬水體運(yùn)動(dòng)狀態(tài)，本文采用RNG k-ε 模型封閉水流運(yùn)動(dòng)控制方程。其方程式如下所示：

式中：kT為紊動(dòng)動(dòng)能，J；PT為紊動(dòng)動(dòng)能的表征項(xiàng)；GT是浮力作用下的紊動(dòng)動(dòng)能產(chǎn)生項(xiàng)，在不可壓縮流體中默認(rèn)為0；εT為紊流能量消散率；CDIS為無量綱參數(shù)，DiffT和DiffS為能量消散項(xiàng)。

2.3 方程的離散與求解

采用有限差分法進(jìn)行數(shù)值離散，將空間計(jì)算域離散成若干個(gè)結(jié)構(gòu)化矩形網(wǎng)格單元求解方程。計(jì)算網(wǎng)格為交錯(cuò)網(wǎng)格，壓強(qiáng)、流體體積分?jǐn)?shù)、渦黏滯系數(shù)、湍流動(dòng)能、耗散率等標(biāo)量均定義在控制體的中心，速度等矢量則定義在垂直于矢量方向的控制體表面中心。

方程離散后，選用GMRES 迭代算法進(jìn)行壓力迭代求解。動(dòng)量方程中的求解，本文選擇顯式解法中二階精度格式進(jìn)行離散，該格式能夠更好地模擬水流的變形破碎。

2.4 模擬范圍與網(wǎng)格劃分

計(jì)算模型范圍在X方向（水流方向）以丁壩為中心的長(zhǎng)度取20L，Y方向取5L，Z方向取H1。X、Y方向采用非均勻矩形網(wǎng)格劃分，并對(duì)丁壩附近網(wǎng)格進(jìn)行局部加密計(jì)算。其中非加密部分X、Y方向網(wǎng)格尺度最大取0.2L， 加密部分X、Y方向網(wǎng)格尺度最大取0.1L（見圖1）。Z方向網(wǎng)格為均勻網(wǎng)格，網(wǎng)格尺度取0.1H1。（L為壩長(zhǎng)，H1為壩高）。

圖1 XY 平面網(wǎng)格劃分示意圖

2.5 邊界條件

數(shù)值模擬水池左側(cè)即Xmin使用速度邊界條件，通過給定邊界流體的流速和水位來模擬生成恒定流以及涌潮；右側(cè)Xmax使用出流邊界；頂端Zmax使用壓力邊界條件，底部Zmin為墻邊界條件；前側(cè)Ymin為墻邊界條件，后側(cè)Ymax為對(duì)稱邊界條件。

3 模型驗(yàn)證

涌潮作用下的丁壩水力特性物理模型試驗(yàn)存在參數(shù)測(cè)量難度大且精度不高的問題，因此，本文采用恒定流作用下的丁壩水流結(jié)構(gòu)進(jìn)行數(shù)學(xué)模型驗(yàn)證。驗(yàn)證試驗(yàn)在浙江省河口海岸重點(diǎn)實(shí)驗(yàn)室的玻璃水槽中進(jìn)行。水槽長(zhǎng)50.0 m，寬1.2 m，高0.6 m。試驗(yàn)所用丁壩長(zhǎng)度L=40 cm，于水槽中間段中心位置對(duì)稱布置，試驗(yàn)水深d=20 cm。丁壩由厚度為10 mm的有機(jī)玻璃制成，其高度取22 cm，略高于水位以模擬非淹沒式剛性壩。試驗(yàn)立面示意見圖2。

圖2 試驗(yàn)立面示意圖

試驗(yàn)水流條件見表1。流速測(cè)量采用“小威龍”剖面流速儀，縱向流速測(cè)點(diǎn)斷面布置依據(jù)距離結(jié)構(gòu)物近則密、遠(yuǎn)則疏的原則，并考慮流場(chǎng)對(duì)稱性。為提高效率，橫向取1/2 結(jié)構(gòu)布置測(cè)點(diǎn)，從中心點(diǎn)起沿水流方向（縱向）布置0～4 共5 個(gè)測(cè)流斷面，斷面間距均為10 cm。每個(gè)縱向斷面的上游布置S1～S3 共3 個(gè)斷面，下游布置0～9 共10 個(gè)斷面，斷面間距及具體布置見圖3。

表 1 試驗(yàn)水流條件表

圖3 流速測(cè)點(diǎn)平面布置圖 單位：cm

各工況下，表層斷面（z=18 cm）處試驗(yàn)測(cè)量結(jié)果與數(shù)值模擬結(jié)果的速度云圖對(duì)比見圖4。從圖4 中可以看出，數(shù)值模擬與試驗(yàn)的流速分布情況及規(guī)律基本相同，由于丁壩的挑流作用及過水?dāng)嗝婵s窄，壩頭附近流速均顯著增大。同時(shí)由于丁壩的掩護(hù)作用，下游丁壩正后方區(qū)域流速顯著減小。由于試驗(yàn)測(cè)量斷面數(shù)據(jù)插值與數(shù)值模擬數(shù)據(jù)插值相比較為稀疏，因此壩頭附近的流速與實(shí)測(cè)流速有一定差異，壩下游側(cè)的回流區(qū)也有一定誤差?？傮w而言，本文建立的數(shù)學(xué)模型能夠相對(duì)準(zhǔn)確地模擬丁壩影響下的水力特性。

圖4 表層斷面速度云圖對(duì)比圖

4 水力特性分析

采用參考文獻(xiàn)7 中的方法及數(shù)值生成涌潮，并在此基礎(chǔ)上模擬涌潮作用下非淹沒丁壩附近的水力特性情況。以潮前水深d=10 cm、潮頭高度H=5 cm 及傳播速度c=0.45 m/s 的涌潮為例分析說明。

4.1 水位分析

涌潮作用下的水位分布情況見圖5。由圖5 可知，涌潮在潮頭剛剛到達(dá)丁壩時(shí)，由于丁壩的阻流效應(yīng)，迎潮側(cè)水位液面快速抬高，最大壅高可達(dá)2H。背潮側(cè)水位未發(fā)生明顯變化，僅接近丁壩的潮頭線明顯晚于其他位置處的潮頭線，且越接近壩根處，潮頭線的滯后時(shí)間越長(zhǎng)。

圖5 涌潮作用下水位分布圖

涌潮潮頭經(jīng)過丁壩后，迎潮側(cè)水位壅高明顯減小，最大壅高僅0.8H。涌潮作用下，背潮側(cè)的回流區(qū)范圍內(nèi)水位下降明顯，最大下降幅度達(dá)H。

4.2 流速分析

涌潮作用下的平均流速分布情況見圖6。由圖6 可知，當(dāng)涌潮潮頭剛剛到達(dá)丁壩時(shí)，丁壩的束流效果較弱，僅在沿丁壩長(zhǎng)度方向L范圍內(nèi)的水流產(chǎn)生繞流，流向發(fā)生變化。流速大小發(fā)生變化的范圍基本在丁壩附近，壩頭附近局部范圍內(nèi)流速增大，最大流速達(dá)1.7c。丁壩迎潮側(cè)水流受丁壩壅阻，流速減緩。背潮側(cè)流速分布基本以潮頭線為界，潮頭線至丁壩范圍內(nèi)流速為涌潮傳播速度，潮頭線至丁壩范圍外流速接近于0。

圖6 涌潮作用下平均流速分布圖

在涌潮潮頭經(jīng)過丁壩后，丁壩的束流效果增強(qiáng)，沿丁壩長(zhǎng)度方向4L范圍內(nèi)的水流流向發(fā)生變化，且在丁壩背潮側(cè)產(chǎn)生回流漩渦。壩頭背潮側(cè)附近流速增大的范圍相比剛到達(dá)時(shí)大幅增加，沿涌潮傳播方向長(zhǎng)度達(dá)L，垂直涌潮傳播方向長(zhǎng)度達(dá)1.2L，該范圍內(nèi)的最大流速略減小至1.5c，且丁壩壩根附近迎背側(cè)均存在流速接近于0 的區(qū)域。

4.3 近底切應(yīng)力分析

涌潮作用下的近底切應(yīng)力分布情況見圖7。由圖7 可知，在涌潮潮頭剛剛到達(dá)丁壩時(shí)，切應(yīng)力與流速分布相似，背潮側(cè)及靠近丁壩1.5L范圍內(nèi)的迎潮側(cè)，切應(yīng)力接近于0，其余位置處的切應(yīng)力達(dá)0.4 Pa，僅壩頭位置局部范圍內(nèi)的切應(yīng)力較大，近底最大切應(yīng)力達(dá)1.7 Pa。

圖7 涌潮作用下近底切應(yīng)力分布圖

涌潮潮頭經(jīng)過丁壩后，在沿丁壩長(zhǎng)度方向L范圍內(nèi)的切應(yīng)力基本都接近于0。除壩頭背潮側(cè)位置沿涌潮傳播方向長(zhǎng)度達(dá)L，垂直涌潮傳播方向長(zhǎng)度達(dá)1.5L的范圍內(nèi)切應(yīng)力相對(duì)較大，平均切應(yīng)力達(dá)1.1 Pa，最大紊動(dòng)切應(yīng)力達(dá)1.9 Pa，其余位置處的切應(yīng)力基本在0.4 Pa 以內(nèi)。由此可見，壩頭背潮側(cè)位置切應(yīng)力相對(duì)較大，泥沙易在此起動(dòng)形成沖刷作用，而迎潮側(cè)則由于切應(yīng)力始終較小，水流攜帶的泥沙易在此沉降形成淤積作用。

5 結(jié) 語

在利用物理模型試驗(yàn)成果對(duì)數(shù)學(xué)模型進(jìn)行驗(yàn)證的基礎(chǔ)上，模擬研究了涌潮作用下非淹沒丁壩的水位、流速、切應(yīng)力分布情況，得到以下主要結(jié)論：

涌潮潮頭剛到達(dá)丁壩時(shí)，多表現(xiàn)為波浪的特性，水位在丁壩迎潮側(cè)有明顯壅高，流速及切應(yīng)力分布以潮頭線為界，僅壩頭局部范圍內(nèi)擾流數(shù)值較大；涌潮潮頭經(jīng)過丁壩后，多表現(xiàn)為水流的特性，丁壩迎潮側(cè)水位壅高，背潮側(cè)水位下降，流速及切應(yīng)力在壩頭背潮側(cè)沿涌潮傳播方向L，垂直涌潮傳播方向長(zhǎng)度1.5L的范圍內(nèi)較大，這也是涌潮作用下丁壩壩頭沖刷需要重點(diǎn)防護(hù)的區(qū)域。

由于實(shí)際工程中丁壩的結(jié)構(gòu)形式及布置方式各有不同，本次模擬的水力特性與實(shí)際有一定差異，對(duì)丁壩群等與工程實(shí)際更為符合的情況有待進(jìn)一步研究。