论文 黄河小浪底水库三维水沙数学模型初步研究.docVIP

黄河小浪底水库三维水沙数学模型初步研究 －建模和率定 摘 要 本文首先通过转换得到贴体坐标系下三维动量方程和方程；然后，采用基于交错网格的有限体积法对其离散，再采用SIMPLE-C方法对速度场和压力场进行耦合求解，所有代数方程都采用TDMA方法求解；通过求解基于二维水深方程的Poisson方程来确定自由面，固体边界则采用壁面函数法处理；最后，以三维方腔流为算例，对程序的可靠性进行验证。 关 键 词 贴体坐标；有限体积法；壁面函数；三维数学模型；方腔流 1 前言 目前小浪底水库调水调沙调度运用亟待解决的工程泥沙是不同泄水建筑物组合运用条件下坝区输沙流态和出库水流、含沙量、级配过程预测分析等。该问题具有强烈的三维性，一维、平面二维和立面二维模型能反映断面平均、垂线平均或横向平均的水流泥沙运动情况，但不能充分反映具有强烈三维性的小浪底水库坝区水沙运动变化规律。因此，要全面准确预测分析不同泄水建筑物运用条件下出库水沙过程，满足小浪底水库调水调沙实时调度运用，必须尽快建立小浪底库区三维水沙模型。 2 控制方程离散与求解 2.1 贴体坐标下的控制方程 通过贴体坐标系统（Body-Fitted Coordinates）[1]可以将复杂的物理域的流动问题转换到规则、简单的计算域中进行计算。在计算域中，实现复杂边界条件下的水沙流数值模拟，相应贴体坐标下控制方程可以表示为下列通用的张量形式： 通用控制方程中主要符号的具体形式见表1。 表1 通用控制方程中各符号的具体形式 式中：为时间；为物理坐标；为计算坐标；为水流密度；为流速；分别代表。 2.2 控制方程离散 采用有限体积法离散控制方程，其中变量布置采用交错网格。取图1所示的控制体积，通用控制方程在节点P的离散方程为： 采用混合格式[2]离散对流项，该格式虽然只有一阶精度，但是一种绝对稳定的离散格式，精度可以满足工程需要。扩散项和源项采用中心差分离散，最终得到各系数表达式如下： 图1 控制体积示意图 ，； ，； ，； 其中，为对流强度，表达式如下： ，，； ，，； 是扩导系数，其表达式如下： 其中，（m=1,2,3；k=1,2,3；j=1,2,3）。 2.4 初边界条件 2.4.1 自由面边界 自由表面的运动对非恒定流的泥沙输移起着不可忽略的作用，泥沙的自由表面条件应根据自由表面的上浮通量与沉降通量相平衡给出的，仅当自由表面不变化，或变化十分缓慢的情况下，自由表面法向才与垂向重合。事实上，若床面形态复杂、地形变化较大时，自由表面均会产生运动，甚至出现剧烈的波动，导致自由表面发生一定的弯曲和变形；再者，自由表面的运动，会直接引起水流的静水压力的变化，进而对水体的流场，特别是二次流的精细结构产生影响。进行水流的数值计算时，由于自由表面未知，不能直接将自由表面的大气压力作为一个主要的自由水面边界条件。所以，得到的压力场类似于一个封闭管中的压力场，在计算过程中无法直接得到自由水面的位置。近十年来，处理自由表面问题主要有标记结点法、空隙比法和标高函数法等方法，本模型采用自由面位置的Poisson方程求解自由面位置[2]。 2.4.2 壁面边界 采用壁面函数法[3]，即在粘性底层内不布置任何节点，把靠近壁面的第一个节点布置在粘性底层之外的完全湍流区，要求第一个计算节点与壁面间的无因次距离（为第一个节点距离壁面的距离）在30～100之间。设计算壁面相邻的第一个节点到壁面的无量纲距离，定义摩阻速度，则计算边界上平行于壁面的流速满足对数关系式： 第一个内点处湍动能和耗散率分别为： 2.4.3 床面附近含沙量 已知床面近邻某一节点的含沙量，则床面附近的含沙量： 式中：为床面附近的挟沙力；为床面泥沙交换层厚度，水库中取值0.005m，河道中取沙波厚度的。 2.4.4 床面附近挟沙力 由张瑞瑾从扩散理论出发提出了含沙量沿垂线分布公式可以得出床面附近挟沙力： （4-6） 式中：；分别为距离水面和床面的相对位置；为含沙量沿垂线分布函数；为相对水深处第粒径组水流挟沙力，其计算公式为： 式中：为系数、指数；为浑水容重；为泥沙容重；为垂线平均流速；垂线平均水深；为悬移质级配；为浑水沉速。 2.3 控制方程求解 求解u,v,w的关键是如何求解压力场，本文采用压力修正法[1]来得到压力场，该方法的基本思想为：对于给定的压力场，按次序求解u,v,w的代数方程，由此得到的速度场未必满足质量守恒的要求，因而必须对压力场进行修正。为了避免由于速度场和压力场分开求解可能带来的压力震荡问题，采用如图2所示的交错网格对速度场和压力场进行求解。 图2 交错网格布置示意图 （a）主控制体 （b）控制体 （c）控制