張 麗

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

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

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

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

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

邊界擬合坐標(biāo)系下的水流運動方程為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖 7 中 A1，A2為相鄰上、下游斷面；θ1，θ2分別為相鄰上、下游過流斷面的法線方向角；Fx，F(xiàn)y分別為相鄰過流斷面間岸線對水體的約束力；θ為相鄰過流斷面間岸線切線的角度；定義當(dāng)θ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兩個坐標(biāo)方向的投影；∑Fx，∑Fy為作用在上、下游過流斷面間流體上的所有外力在x，y兩個坐標(biāo)方向投影的代數(shù)和。

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

圖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é)模型，并對河北省境內(nèi)的黎河西梁子河段水流情況進行模擬，計算結(jié)果與實測資料吻合較好，說明模型很好地反映了原型黎河西梁子河段的水流運動特性，也說明模型理論和計算方法是正確的。同時根據(jù)二維恒定流動量方程，提出彎道岸線沖擊力、剪切力的計算方法，以此為理論基礎(chǔ)，計算黎河西梁子河段岸線應(yīng)力沿程分布，為彎曲河道的岸線防護提供理論依據(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ù)學(xué)模型研究［J］.華北水利水電學(xué)院學(xué)報，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ù)雜邊界河道流速場的數(shù)值模擬［J］.水利學(xué)報，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］.天津：天津大學(xué)出版社，1994.

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

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

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