基于σ坐标变换和水位函数法的立面二维水动力学模型

孙 先，杨文俊，王国栋

(长江科学院 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中特定时刻速度场和动水压强场分布模型结果和解析解两者的值大致相同，变化趋势基本吻合。可知本文所建立的立面二维水动力学具有良好的精度。

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))