1、- 1 - 黄河小浪底水库三维水沙数学模型初步研究建模和率定摘 要 本文首先通过转换得到贴体坐标系下三维动量方程和 方程；然后，采用k基于交错网格的有限体积法对其离散，再采用 SIMPLE-C 方法对速度场和压力场进行耦合求解，所有代数方程都采用 TDMA 方法求解；通过求解基于二维水深方程的Poisson 方程来确定自由面，固体边界则采用壁面函数法处理；最后，以三维方腔流为算例，对程序的可靠性进行验证。关 键 词 贴体坐标；有限体积法；壁面函数；三维数学模型；方腔流1 前言目前小浪底水库调水调沙调度运用亟待解决的工程泥沙是不同泄水建筑物组合运用条件下坝区输沙流态和出库水流、含沙量、级配过程预
2、测分析等。该问题具有强烈的三维性，一维、平面二维和立面二维模型能反映断面平均、垂线平均或横向平均的水流泥沙运动情况，但不能充分反映具有强烈三维性的小浪底水库坝区水沙运动变化规律。因此，要全面准确预测分析不同泄水建筑物运用条件下出库水沙过程，满足小浪底水库调水调沙实时调度运用，必须尽快建立小浪底库区三维水沙模型。2 控制方程离散与求解2.1 贴体坐标下的控制方程通过贴体坐标系统（Body-Fitted Coordinates） 1可以将复杂的物理域的流动问题转换到规则、简单的计算域中进行计算。在计算域中，实现复杂边界条件下的水沙流数值模拟，相应贴体坐标下控制方程可以表示为下列通用的张量形式： (
3、)(,)kmkjUStx通用控制方程中主要符号的具体形式见表 1。表 1 通用控制方程中各符号的具体形式- 2 - 式中： 为时间； 为物理坐标； 为计算坐标； 为水流密度； 为流速； 分别tjxkju代表 。 ,kwvu2.2 控制方程离散采用有限体积法离散控制方程，其中变量布置采用交错网格。取图 1 所示的控制体积，通用控制方程在节点 P 的离散方程为：PEWNSTBAA采用混合格式 2离散对流项，该格式虽然只有一阶精度，但是一种绝对稳定的离散格式，精度可以满足工程需要。扩散项和源项采用中心差分离散，最终得到各系数表达式如下：图 1 控制体积示意图， ；|,2eeEFAD|,2wWFA，
4、；nnNssS， ；|,ttT|,bbB0()PEWNSTPewnstbFF其中， 为对流强度，表达式如下：， ， ；()eeU()w()nV， ， ；ssVttb- 3 - 是扩导系数，其表达式如下：D1121233()0.5().().2eenstbgg1121233()0.5().().wnstbDgg205()()ss wt2205()()ssewt331312.().t ensg3131322.().().()bbewnsDgg其中， （m=1,2,3；k=1,2,3；j=1,2,3） 。mkkjjx2.4 初边界条件2.4.1 自由面边界自由表面的运动对非恒定流的泥沙输移起着不可忽
5、略的作用，泥沙的自由表面条件应根据自由表面的上浮通量与沉降通量相平衡给出的，仅当自由表面不变化，或变化十分缓慢的情况下，自由表面法向才与垂向重合。事实上，若床面形态复杂、地形变化较大时，自由表面均会产生运动，甚至出现剧烈的波动，导致自由表面发生一定的弯曲和变形；再者，自由表面的运动，会直接引起水流的静水压力的变化，进而对水体的流场，特别是二次流的精细结构产生影响。进行水流的数值计算时，由于自由表面未知，不能直接将自由表面的大气压力作为一个主要的自由水面边界条件。所以，得到的压力场类似于一个封闭管中的压力场，在计算过程中无法直接得到自由水面的位置。近十年来，处理自由表面问题主要有标记结点法、空隙
6、比法和标高函数法等方法，本模型采用自由面位置的 Poisson 方程求解自由面位置 2。2.4.2 壁面边界采用壁面函数法 3，即在粘性底层内不布置任何节点，把靠近壁面的第一个节点布置在粘性底层之外的完全湍流区，要求第一个计算节点与壁面间的无因次距离（ 为第一个节点距离壁面的距离）在 30100 之间。设计算壁面相邻的*PzuPz第一个节点到壁面的无量纲距离 ，定义摩阻速度 ，则计算边界*Pzu1/42*uck上平行于壁面的流速 满足对数关系式：u*2.5ln.uz- 4 - 第一个内点 处湍动能和耗散率分别为：P12*kuc13*Ppuz2.4.3 床面附近含沙量已知床面近邻某一节点 的含沙
7、量 ，则床面附近的含沙量：LLlS()*1elsbtzblLlblS式中： 为床面附近的挟沙力； 为床面泥沙交换层厚度，水库中取值*bl b0.005m，河道中取沙波厚度的 。232.4.4 床面附近挟沙力由张瑞瑾从扩散理论出发提出了含沙量沿垂线分布公式可以得出床面附近挟沙力：（4-6） * *()expsbsbl lkSSffdu式中： 为相对水深，其表达式为 ； 、 分别为距离水面和床面的1zHsb相对位置； 为含沙量沿垂线分布函数； 为相对水深 处第 粒径组水流挟沙力，f*lSl其计算公式为： 3*0mml SllsLsUSkPgh式中： 、 为系数、指数； 为浑水容重； 为泥沙容重；
8、为垂线平均流sU速； 垂线平均水深； 为悬移质级配； 为浑水沉速。hSl sl2.3 控制方程求解求解 u,v,w 的关键是如何求解压力场，本文采用压力修正法 1来得到压力场，该方法的基本思想为：对于给定的压力场，按次序求解 u,v,w 的代数方程，由此得到的速度场未必满足质量守恒的要求，因而必须对压力场进行修正。为了避免由于速度场和压力场分开求解可能带来的压力震荡问题，采用如图 2 所示的交错网格对速度场- 5 - 和压力场进行求解。图 2 交错网格布置示意图（a）主控制体 （b） 控制体 （c） 控制体uv设原来的压力为 ，与之相对应的速度为 ， 分别表示压力与*p*,wu,vup速度的修
9、正量，将其带入离散后的动量方程，采用 SIMPLE-C 方法可以得到速度修正值的计算公式为： PuPxxipApvpyyipwpzziA将速度修正值带入连续方程可以得到压力修正方程为：式中：WNSTBPE PApC, , ,ueABuwvNn, ,vSsTtwBbPEWNTA*ewnstbCUV uvwpppyxziiiBAyvxziiiAuvwpppywxziii- 6 - 上式中的 、 分别为速度 u、v、w 的主节点系数。upAvpw模型采用三对角方程组的解法 TDMA（Tri-Diagonal Matrix Algorithm）进行迭代求解。3 模型率定3.1 丁坝绕流水槽的长度为 8
10、m，宽度为 30cm。丁坝位于 x=4m 的位置，丁坝尺寸为15cm3cm5cm。入流流量为 3600m3/s，水深为 9cm。计算模拟时，三维模型的网格点在 x、y、z 为 2513119，丁坝附近网格逐渐加密。下图为 t=24s 时，z=2、4cm截面上丁坝附近流场分布。图 3 不同截面上丁坝附近流场分布图3.2 三维方腔流三维方腔流计算区域大小为 1m1m1m，雷诺数 Re=100，划分为 414141 个均匀网格。计算边界条件：z=0 时，u=1,v=0,w=0；x=0,x=1,y=0,y=1,z=1 时，均为无滑移边界条件，即 u=0,v=0,w=0.图 4 为在 x=0.45，y=
11、0.45 处速度 u 的计算结果与 Chen4的计算结果比较图；图5 为截面 x=0.45 的速度矢量图。分析表明水流计算稳定性较好，计算精度可以满足要求。- 7 - 图 4 x=0.45，y=0.5 处速度 u 沿 z 方向分布图 5 x=0.45 处的流速矢量图4 结论与建议（1）本文所建立的贴体坐标系下的小浪底水库水库模型对复杂边界有着良好的适应性，且针对水库模型的特殊性，采用绝对稳定的混合格式离散对流项，从而保证了数值计算的稳定。（2） 通过三维方腔流的计算结果与 Chen 的结果进行对比，证明模型有着良好的稳定性，且精度可以满足要求。主要参考文献1 王福军.著.计算流体动力学分析CF
12、D 软件原理与应用.清华大学出版社，20042 陶文铨，著.数值传热学(第 2 版).西安交通大学出版社. 2002.9.3 Chen C J, Ho K S, Cheng W C. The Finite Analytic Method. IIHR Report. 1982, 232-V, The University of Iowa, Iowa City, USA, 88-100.4 唐学林.水流泥沙运动的三维数值模拟研究: 博士后研究报告.清华大学，2005基金项目：“十一五”国家科技支撑计划重点项目（2006BAB06B02）资助作者简介：余欣（1971-） ，男，河南南阳人，高级工程师
13、，黄河数学模拟系统首席专家，主要从事河流泥沙动力学研究- 8 - - 9 - Primary Research of Three-dimensional Mathematical Model of Flow and Sediment of Xiaolangdi Reservoir on Yellow River-Modeling and CalibrationYu Xin1,2 Wang Ming1 Tang Xuelin2 Wang Min1(1. Yellow River Institute of Hydraulic Research, Zhengzhou , 450003 2. Chin
14、a Agriculture University, Beijing , 100083)Abstract: This paper first gets the 3-D momentum equation and k equation under body-fitted coordinates through conversion; then, use finite volume method based on staggered grid to discretize it; and then use SIMPLE-C method to seek coupled solution of velo
15、city and pressure field. All the algebraic equations are solved through TDMA. Determine the free surface through the solution of Poisson Equation based on 2-D water depth equation; apply wall function to processing solid boundary; finally, use 3-D lid-driven cavity flow as the example to verify the
16、reliability of the program. Key Words: body-fitted coordinates; finite volume method; wall function; 3-D mathematical model; lid-driven cavity flow 10 1 ForewordThe engineering sediment that needs urgent handling in the regulation of flow and sediment at Xiaolangdi Reservoir is the comprehensive exe
17、rcise of forecast and analysis of discharged sediment flow, discharged water flow and grading process and so on with the combined use of the different regulating buildings in dam area. This issue has strong 3-D features. 1-D, plane 2-D and vertical 2-D models could reflect the moving of the cross-se
18、ction average, depth average and horizontal average water flow and sediment, but they cannot fully reflect the rule of the water and sediment movement change in the dam area of Xiaolangdi Reservoir, which has strong 3-D features. Therefore, in order to fully and accurately forecast and analyze the d
19、ischarge process of flow and sediment through different discharge buildings and meet the need of the actual scheduling and exercise of the regulation of flow and sediment at Xiaolangdi Reservoir, it is urgent to build up the 3-D model of flow and sediment at the dam area of Xiaolangdi Reservoir. 2 T
20、he Discretization and Solution of Governing Equation2.1 Governing Equation under Body-fitted CoordinatesThe flow question in the complicated physical domain could be conversed into regular and simple computational domain for computation through body-fitted coordinates 1. The simulation of water and
21、sediment flow value under complicated boundary conditions could be realized in computational domain; the corresponding governing equation under body-fitted coordinates could be presented by universal tensor form as below: ()(,)kmkjUStxThe specific forms of the major symbols in the universal governing equation are listed in the table-1 Table-1 the specific forms of the symbols in universal governing equation In the equation: t stands for time, jx for physical coordinate, k for computational coordinate, for water flow density, ju for flow velocity; for ,kwvu respectively.
