基于σ坐標變換和水位函數法的立面二維水動力學模型

孫 先，楊文俊，王國棟

(長江科學院 a.研究生部;b.科研計劃處，武漢 430010)

立面二維水流數學模型能給出水流運動沿水深方向變化的特征，相對于平面二維模型只能反映出河床平面上的水流地形變化［1］，對于研究寬深比較小的河道和水庫中的水流物質輸運有著不可比擬的優勢。因其求解涉及到3種介質2個交界面，網格劃分、自由水面的追蹤、動水壓力的求解3個問題成為了模型建立的難點。

σ坐標變換由 Phillips［2］在1957年提出，它能夠很好地擬合河底地形，并隨著地形和水面的變化而實時調整。通過σ坐標變換后，可以將原本不規則的計算區域變為一個固定的矩形，給網格剖分和編程都帶來了方便。將笛卡爾坐標系下的立面二維水動力學控制方程轉化到σ坐標系下，然后在交錯網格上，采用水位函數法求解自由水面，運用有限差分法分步離散求解控制方程。最后使用模型模擬微幅波運動，通過模擬值與理論值的比較，驗證模型的可靠性。

1 σ坐標系下的立面二維水動力學控制方程

基于Boussinesq假設的不可壓縮牛頓流體的Navier-stokes方程，對流體運動基本方程沿橫向積分平均后，可得笛卡爾坐標系下立面二維水動力運動基本方程組:

式中:u，w分別為x*方向和z*方向的流速;Pd為動水壓力;ε為紊動擴散系數;ξ為自由水面;ρ為水體密度;g為重力加速度。

立面二維計算區域的σ坐標變換如圖1所示。

圖1 立面二維計算區域的σ坐標變換Fig.1 Transformation of coordinate σ of vertical 2-D computational domain

由以上微分關系可得，σ坐標系下法向流速為

對于笛卡爾坐標系下連續性方程通過以上坐標變換得

x方向水流運動方程變為

z方向水流運動方程變為

2 控制方程的離散以及求解

本節中，采用分步解法，在交錯網格上對σ坐標系下的控制方程進行離散。

2.1 分步解法求解x方向動量方程

在本步驟中忽略動水壓力Pd以及水平擴散作用項的影響，原x方向動量方程式(7)可以按照分別考慮對流作用、擴散作用以及靜水壓力作用分步求解，具體過程如下:

(1)僅考慮對流作用下的方程離散，

(2)考慮垂向擴散作用下的方程離散，

(3)考慮靜水壓力作用下的方程離散，

由式(13)易得

2.2 水位函數法求解沿程水深

對連續性方程(6)，沿垂向積分可得水位方程

2.3 求解動水壓力

整理得

對比式(8)可得

整理得

將式(11)代入笛卡爾坐標系下連續性方程

可得σ坐標系下連續性方程的另一種形式，即

在標量控制體上離散該方程可得

式中:

定義

CDX反映了笛卡爾坐標系下計算網格的不均勻程度。計算動水壓力時，不考慮河床地形和水深隨時間變化，因此可忽略式(18)中含叉項，將式(18)、式(20)代入式(23)，整理后可寫成如下形式:

其中:

3 模型驗證

經典的微幅波具有解析解，被廣泛應用于檢測模型模擬非靜水壓力水流運動的能力及模擬非恒定水流運動時的質量、能量守恒特性。是否能成功模擬微幅波運動，是區別靜水模型和非靜水模型的重要標志之一［4］。

微幅波運動屬于勢流范疇，其勢函數和彌散方程分別為［5］:

式中:?為勢函數;L為波高，即波峰到波谷的距離;ω為角頻率;g為重力加速度;k為波數;H為水深。

?對橫坐標x*、垂向坐標z*求導可得速度場:

式中，T為微幅波的振動周期。同時，動水壓強分布為

如圖2所示，在封閉的方槽中有一微幅uninodal波，初始時刻水體靜止并處于最大勢能狀態，以z=0m為平衡水面。在開始運動后，當水面運動至水平時，水體的所有勢能均轉化為動能，接著它又在另一邊形成最大勢能狀態然后再返回，如此周期運動。圖中方槽長L=10m，水深H=10m，初始水面形態為

式中:A為振幅;k=2π/nL，對于uni-nodal波 n=2。取A=0.1m，為水深的1%，滿足微幅波理論的適用條件。當不計動水壓力，按照淺水波計算時，波速，微幅波的振動周期為2 s;當考慮動水壓力，按照深水波計算，由式(35)可得，角頻率ω=1.769 rad/s，進而可得微幅波運動的波速、周期分別為c=ω/k=5.63 m/s，T=2π/ω=3.55 s。

計算網格取Δx*=Δz*=1m，以式(39)為初始水面，根據微幅波理論的假定，在計算中不計對流項、擴散項及底面阻力。時間步長Δt*=0.01s，經過4s計算，得到x*=0處的水位變化過程如圖3所示。

圖2 封閉方槽中的微幅波Fig.2 Micro-amplitude wave in the closed square groove

圖3 x*=0處水位變化模型解與解析解對比圖Fig.3 Comparison of water level fluctuation between model solution and theoretical solution when x*=0

在使用非靜水模型計算時，T/4，5T/8時刻的流場、動水壓強如下圖4所示。

通過圖3解析解與模型計算值之間的對比可以得到，x*=0水位變化計算結果基本與解析解之間相對誤差控制在10%以內。由圖4中特定時刻速度場和動水壓強場分布模型結果和解析解兩者的值大致相同，變化趨勢基本吻合?？芍疚乃⒌牧⒚娑S水動力學具有良好的精度。

4 結論

本文建立了動水壓力條件下的立面二維數學模型，并運用該模型對微幅波運動進行模擬，模擬的結果與理論計算值吻合良好，證明此模型可以模擬較復雜的存在動水壓力的水流情況。不過結果仍存在一定誤差。誤差的主要來源是:①離散過程中差分格式帶來的誤差;②計算過程中對CDX等項的忽略。這些誤差均可以通過采用高階差值函數和改進算法減小。另外，天然水體均是固液兩相流，本文僅探討了水流模型，在結合泥沙模型后才能準確地模擬天然水體。

圖4 特定時刻速度、動水壓強解析解與模型計算值的比較Fig.4 Comparison of speed and hydrodynamic pressure between theoretical solution and model solution in specific moments

［1］任坤杰，金 峰，徐勤勤.滑坡涌浪垂面二維數值模擬［J］.長江科學院院報，2006，23(2):1－4.(REN Kun-jie，JIN Feng，XU Qin-qin.Vertical Two-Dimensional Numerical Simulation for Landslide-Generated Waves［J］.Journal of Yangtze River Scientific Research Institute，2006，23(2):1－4.(in Chinese))

［2］PHILLIPSN A.A Coordinate System Having Some Special Advantages for Numerical Forecasting［J］.Journal of Meteorology，1957，14:184－185.

［3］馮小香.水庫壩前沖刷漏斗形態數值模擬研究［D］.武漢:武漢大學，2006.(FENG Xiao-xiang.Numerical Simulation of Scoured Funnel in Front of Dam［D］.Wuhan:Wuhan University，2006.(in Chinese))

［4］胡德超.三維水沙運動及河床變形數學模型研究［D］.北京:清華大學，2009.(HU De-chao.Study on the Three-dimensional Numerical Model for Free-surface Flow and Suspended Sediment Transport［D］.Beijing:Tsinghua University，2009.(in Chinese))

［5］吳宋仁.海岸動力學［M］.北京:人民交通出版社，2000.(WU Song-ren.Coastal Dynamics［M］.Beijing:China Communications Press，2000.(in Chinese))