彎道水流二維數(shù)值模擬及岸線整治的水動(dòng)力分析方法研究

張 麗

（交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究所工程泥沙交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，天津300456）

彎曲河道的邊界極不規(guī)則，岸線蜿蜒曲折，在數(shù)值計(jì)算時(shí)邊界擬合困難，很難保證計(jì)算精度。為克服這一點(diǎn)，很多學(xué)者做了研究，尋求解決處理復(fù)雜邊界的辦法。邊界擬合坐標(biāo)變換法［1］是目前應(yīng)用較多的處理復(fù)雜邊界的方法。本文結(jié)合前人的計(jì)算經(jīng)驗(yàn)［2-3］，采用邊界擬合坐標(biāo)變換技術(shù)，建立邊界擬合坐標(biāo)系下平面二維水流數(shù)學(xué)模型，并根據(jù)二維恒定流動(dòng)量方程，提出彎道岸線沖擊力、剪切力的計(jì)算方法，以此為理論基礎(chǔ)，計(jì)算黎河西梁子河段岸線應(yīng)力沿程分布，為彎曲河道的岸線防護(hù)提供理論依據(jù)。

1 邊界擬合坐標(biāo)系下平面二維水流運(yùn)動(dòng)的基本方程

本文利用邊界擬合坐標(biāo)技術(shù)，將直角坐標(biāo)系下的控制方程轉(zhuǎn)化為邊界擬合坐標(biāo)系下的控制方程。

邊界擬合坐標(biāo)系下的水流連續(xù)方程為

邊界擬合坐標(biāo)系下的水流運(yùn)動(dòng)方程為

式中：u，v分別為流速在ξ，η方向的分量；z為水位，z=h+zb；h為水深；zb為床面高程；n為糙率系數(shù)；g為重力加速度；ε為紊動(dòng)粘滯系數(shù)。

2 邊界擬合坐標(biāo)系下水流控制方程的離散與求解

本文采用有限差分法對(duì)各方程進(jìn)行離散，考慮到對(duì)流項(xiàng)在流動(dòng)方向上對(duì)導(dǎo)數(shù)差分計(jì)算的影響，因此對(duì)對(duì)流項(xiàng)的離散采用迎風(fēng)格式。對(duì)于每個(gè)方程，在ξ和η方向上均采用交替方向隱式算法（ADI）［4-5］。為避免出現(xiàn)鋸齒形或棋盤形壓力場(chǎng)的不穩(wěn)定結(jié)果，物理變量采用交錯(cuò)網(wǎng)格排列，即把速度u、v及水位z分別布置于3套不同的網(wǎng)格系統(tǒng)。這樣，相鄰的壓力點(diǎn)構(gòu)成了對(duì)其間速度的壓力梯度，就避免了由不均勻壓力場(chǎng)而形成的均勻壓力梯度的結(jié)果（圖1）。

圖1 變量節(jié)點(diǎn)布置示意圖Fig.1 Sketch of variable node

3 模型驗(yàn)證計(jì)算

3.1 計(jì)算區(qū)域及網(wǎng)格剖分

本文選取河北省遵化縣境內(nèi)的黎河西梁子河段40+910斷面至41+860斷面，全長(zhǎng)950 m（包含3個(gè)彎道）為計(jì)算區(qū)域（圖2）。

本文將計(jì)算區(qū)域劃分為128×12的網(wǎng)格圖，最小空間步長(zhǎng)7 m左右，曲線網(wǎng)格見圖3。

圖2 計(jì)算區(qū)域地形圖和水文觀測(cè)斷面位置Fig.2 Topography of calculation area and location of hydrological observation section

圖3 計(jì)算區(qū)域曲線網(wǎng)格圖Fig.3 Curvilinear grid of calculation area

3.2 邊界的處理

邊界擬合坐標(biāo)系下水流控制方程的邊界條件亦需要進(jìn)行推導(dǎo)處理。各種邊界條件下邊界節(jié)點(diǎn)的位置及節(jié)點(diǎn)上物理量的表示情況見圖4。以nΔt→（n+1/2）Δt前半個(gè)時(shí)間段內(nèi)ξ方向?yàn)槔渲蠭S表示左邊界右面的最小整數(shù)內(nèi)點(diǎn)，IE表示右邊界左面的最大整數(shù)內(nèi)點(diǎn)。

（1）左端閉邊界處理。當(dāng)左端為閉邊界時(shí)，此時(shí)左邊界以外節(jié)點(diǎn)的水位為零，且邊界上法線方向的流速為零，即有。

（3）右端閉邊界處理。此時(shí)閉邊界上法線方向的速度為零，即，在右端閉邊界上不需要進(jìn)行其他特別的數(shù)學(xué)處理。

對(duì)于在（n+1/2）Δt→（n+1）Δt時(shí)間段內(nèi)，上、下邊界的處理推導(dǎo)完全類似，此處不再贅述。

圖4 邊界節(jié)點(diǎn)位置示意圖Fig.4 Sketch of boundary nodes location

3.3 變換網(wǎng)格的地形生成

當(dāng)變換網(wǎng)格生成以后，還需要求出變換網(wǎng)格各節(jié)點(diǎn)處所對(duì)應(yīng)的地形高程才能進(jìn)行數(shù)值計(jì)算。由于進(jìn)行了坐標(biāo)變換，網(wǎng)格節(jié)點(diǎn)位置會(huì)發(fā)生較大的偏移，這時(shí)的網(wǎng)格節(jié)點(diǎn)高程就必須通過初始網(wǎng)格各節(jié)點(diǎn)的高程來插值求解。本文選取最短距離——斷面梯度插值法進(jìn)行插值。

對(duì)于河道來說，引入兩岸的邊界條件，初始化區(qū)域內(nèi)網(wǎng)格節(jié)點(diǎn)。此時(shí)初始網(wǎng)格節(jié)點(diǎn)按斷面布置，經(jīng)坐標(biāo)變換后，變換網(wǎng)格節(jié)點(diǎn)的分布也是基本趨于河道斷面形式布置，這樣就主要在河道斷面上，沿?cái)嗝娉跏季W(wǎng)格節(jié)點(diǎn)梯度方向插值求解變換后網(wǎng)格節(jié)點(diǎn)高程。

3.4 數(shù)學(xué)模型驗(yàn)證

針對(duì)河北省境內(nèi)的黎河西梁子河段水域，在2001年測(cè)繪的地形圖（圖2）和同年9月份輸水期60 m/s流量條件下，采用2001年6月24日9個(gè)特征過流斷面實(shí)測(cè)水文資料（觀測(cè)斷面位置見圖2中的加粗?jǐn)嗝妫﹣眚?yàn)證。水流模型的初始條件為u=0，v=0，z=21.85 m；邊界條件為上游給定水位21.98 m，下游給定水位21.81 m，岸邊界處采用不可滲透邊界條件v*=-uyξ+vxξ=0，即法向速度為0。各斷面流速驗(yàn)證結(jié)果［6］見圖5，河段沿程水位驗(yàn)證結(jié)果［6］見圖 6。

從圖5和圖6可以看出模型計(jì)算結(jié)果和原型實(shí)際觀測(cè)資料（流速和水位）吻合很好，斷面平均流速分布基本相同，沿程水位基本相似。可見，模型很好地反映了原型黎河西梁子河段的水流運(yùn)動(dòng)特性，也說明了模型理論和計(jì)算方法是正確的，模型具有一定工程適應(yīng)性。

圖5 各斷面垂線平均流速橫向分布驗(yàn)證圖Fig.5 Verification of average flow velocity

圖6 沿程水位驗(yàn)證圖Fig.6 Verification of water level

4 彎曲河道岸線應(yīng)力分析

河道岸線限制了水流流速的大小和方向，水流受到岸線的約束力，同時(shí)水流也對(duì)岸線施加了沖擊力和剪切力，沿程應(yīng)力的劇烈變化可能會(huì)引起河道岸線的破壞，在彎道上這種劇烈變化出現(xiàn)的可能性尤其大。考慮到流速被河岸突然遏制，必然會(huì)引起動(dòng)量變化，即垂直于河岸的動(dòng)量分量形成對(duì)河岸的沖擊力，平行于河岸的動(dòng)量分量形成對(duì)河岸的剪切力［7］。取河道兩斷面間的彎道，進(jìn)行平面動(dòng)力分析（圖7）。

圖 7 中 A1，A2為相鄰上、下游斷面；θ1，θ2分別為相鄰上、下游過流斷面的法線方向角；Fx，F(xiàn)y分別為相鄰過流斷面間岸線對(duì)水體的約束力；θ為相鄰過流斷面間岸線切線的角度；定義當(dāng)θ1，θ2分別為彎道入口和出口過流斷面時(shí)，θ2-θ1為彎道角度，即河道轉(zhuǎn)彎角度；L為中泓線弧長(zhǎng)。

二維恒定流動(dòng)量方程

圖7 彎道動(dòng)力分析示意圖Fig.7 Sketch of dynamic analysis on channel bends

式中：β為斷面動(dòng)量修正系數(shù)；ρ為流體密度；Q為河道流量；v1x，v1y，v2x，v2y為上、下游過流斷面平均流速v1，v2在x，y兩個(gè)坐標(biāo)方向的投影；∑Fx，∑Fy為作用在上、下游過流斷面間流體上的所有外力在x，y兩個(gè)坐標(biāo)方向投影的代數(shù)和。

在以上理論分析的基礎(chǔ)上，結(jié)合水流數(shù)學(xué)模型計(jì)算結(jié)果，研究黎河河道主流對(duì)岸線的沖擊力和剪切力沿程變化。在60 m3/s流量無護(hù)岸工程條件下，岸線所受水流沖擊力、剪切力沿程分布如圖7所示。結(jié)合文中的圖2，從圖8可以看出，22斷面、108斷面和120斷面附近河道顯著收縮或擴(kuò)張，45斷面、72斷面和90斷面附近河道轉(zhuǎn)彎，這些斷面附近的岸線上沖擊應(yīng)力和剪切應(yīng)力都出現(xiàn)了極值，并發(fā)生了劇烈變化，對(duì)這些區(qū)域的河岸應(yīng)重點(diǎn)加強(qiáng)岸線防護(hù)。

圖8 黎河西梁子河段兩岸的沖擊應(yīng)力和剪切應(yīng)力沿程變化Fig.8 Variation of impact force and shearing force along Xiliangzi channel of Lihe river

5 結(jié)論

本文采用邊界擬合坐標(biāo)變換技術(shù)，建立邊界擬合坐標(biāo)系下平面二維水流數(shù)學(xué)模型，并對(duì)河北省境內(nèi)的黎河西梁子河段水流情況進(jìn)行模擬，計(jì)算結(jié)果與實(shí)測(cè)資料吻合較好，說明模型很好地反映了原型黎河西梁子河段的水流運(yùn)動(dòng)特性，也說明模型理論和計(jì)算方法是正確的。同時(shí)根據(jù)二維恒定流動(dòng)量方程，提出彎道岸線沖擊力、剪切力的計(jì)算方法，以此為理論基礎(chǔ)，計(jì)算黎河西梁子河段岸線應(yīng)力沿程分布，為彎曲河道的岸線防護(hù)提供理論依據(jù)。

［1］Thompson J F.Boundary-fitted coordinate system for numerical solution of partial differential equation［J］.J.Computer Phys.，1982（47）：1-108.

［2］李永強(qiáng)，陳紅中，石守亮.天然河道水沙平面二維數(shù)學(xué)模型研究［J］.華北水利水電學(xué)院學(xué)報(bào)，1999，20（1）：10-13.LI Y Q，CHEN H Z，SHI S L.The study and the application of two-dimension water flow and sediment mathematical model in the river［J］.Journal of North China Institute of Water Conservancy and Hydroelectric Power，1999，20（1）：10-13.

［3］魏文禮，金忠青.復(fù)雜邊界河道流速場(chǎng)的數(shù)值模擬［J］.水利學(xué)報(bào)，1994（11）：26-30.WEN W L，JIN Z Q.Numerical simulation of velocity field for the river of complex boundary［J］.Journal of Hydraulic Engineering，1994（11）：26-30.

［4］曹祖德，王洪運(yùn).水動(dòng)力泥沙數(shù)值模擬［M］.天津：天津大學(xué)出版社，1994.

［5］李義天，趙明登，曹志芳.河道平面二維水沙數(shù)學(xué)模型［M］.北京：中國(guó)水利水電出版社，2001.

［6］陳育權(quán).彎道水流二維數(shù)值模擬及岸線整治的水動(dòng)力學(xué)方法研究［D］.天津：天津大學(xué)，2005.

［7］李大鳴，姜明保，張彤宇.彎曲河道護(hù)岸受力分析的水動(dòng)力學(xué)方法［J］.水運(yùn)工程，2002（7）：23-26.LI D M，JIANG M B，ZHANG T Y.Hydrodynamic Method of Force Analysis of Revetment in Channel Bends［J］.Port&Waterway Engineering，2002（7）：23-26.