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

張 麗

（交通運輸部天津水運工程科學研究所工程泥沙交通行業(yè)重點實驗室，天津300456）

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

1 邊界擬合坐標系下平面二維水流運動的基本方程

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

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

邊界擬合坐標系下的水流運動方程為

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

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

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

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

3 模型驗證計算

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

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

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

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

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

3.2 邊界的處理

邊界擬合坐標系下水流控制方程的邊界條件亦需要進行推導處理。各種邊界條件下邊界節(jié)點的位置及節(jié)點上物理量的表示情況見圖4。以nΔt→（n+1/2）Δt前半個時間段內(nèi)ξ方向為例，其中IS表示左邊界右面的最小整數(shù)內(nèi)點，IE表示右邊界左面的最大整數(shù)內(nèi)點。

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

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

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

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

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

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

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

3.4 數(shù)學模型驗證

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

從圖5和圖6可以看出模型計算結果和原型實際觀測資料（流速和水位）吻合很好，斷面平均流速分布基本相同，沿程水位基本相似?？梢姡Ｐ秃芎玫胤从沉嗽屠韬游髁鹤雍佣蔚乃鬟\動特性，也說明了模型理論和計算方法是正確的，模型具有一定工程適應性。

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

圖6 沿程水位驗證圖Fig.6 Verification of water level

4 彎曲河道岸線應力分析

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

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

二維恒定流動量方程

圖7 彎道動力分析示意圖Fig.7 Sketch of dynamic analysis on channel bends

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

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

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

5 結論

本文采用邊界擬合坐標變換技術，建立邊界擬合坐標系下平面二維水流數(shù)學模型，并對河北省境內(nèi)的黎河西梁子河段水流情況進行模擬，計算結果與實測資料吻合較好，說明模型很好地反映了原型黎河西梁子河段的水流運動特性，也說明模型理論和計算方法是正確的。同時根據(jù)二維恒定流動量方程，提出彎道岸線沖擊力、剪切力的計算方法，以此為理論基礎，計算黎河西梁子河段岸線應力沿程分布，為彎曲河道的岸線防護提供理論依據(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］李永強，陳紅中，石守亮.天然河道水沙平面二維數(shù)學模型研究［J］.華北水利水電學院學報，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］魏文禮，金忠青.復雜邊界河道流速場的數(shù)值模擬［J］.水利學報，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］曹祖德，王洪運.水動力泥沙數(shù)值模擬［M］.天津：天津大學出版社，1994.

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

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

［7］李大鳴，姜明保，張彤宇.彎曲河道護岸受力分析的水動力學方法［J］.水運工程，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.