基于浸没边界法的水库变动水面模拟及验证

2020-03-10程永光吴家阳

水利水运工程学报 2020年1期

周睿，程永光，吴家阳

(武汉大学水资源与水电工程科学国家重点实验室，湖北武汉 430072)

水库水温时空变化预测是水电工程生态环境评价和改善研究领域中重要的研究方向，借助数学工具和数学方法来解决类似问题最经济可行[1]。预测水库水温，应解决水面大幅变动模拟及其计算的问题，可以通过改进浸没边界法加以实现。对于水库水温类型的判别，我国现行水库环境影响评价中普遍采用经验公式法[2]：径流-库容比数法和密度弗劳德数法。也有研究者提出其他一些方法，如蔡为武[3]考虑到水库调节性能、年内泄水状况、泄水孔口相对位置而提出的水库水温分层法；郄志红等[4]提出的一种判别水库水温分层模式的人工神经网络法；刘金禄[5]根据水库水温垂直结构分布的模糊性特点，应用多目标模糊模式识别和回归分析理论，研究了水库水温垂直分布类型的模糊识别回归预测法。国内研究者还提出了许多经验性水温估算方法[6-8]。

目前国内水库水温水质预测的难点在于对速度场和温度场耦合作用、水库形态等因素的处理，模型通用性不足。而新近发展起来的格子Boltzmann方法(Lattice Boltzmann Method，LBM)[9]是一种先进的流体力学数值模拟方法，具有算法简单、并行性好、边界条件易处理和压力直接计算等优点[10]。但是将LBM法用于模拟水库水温的研究比较有限，且模型预测也存在较大误差。因此为了完善水温预测方法，本文利用简便易行、计算效率高的浸没边界法，开展水库水温的时空变化规律研究，模拟真实水库自由液面大幅度变动[11]，并通过简单算例来证实其可靠性。

1 自由液面处理算法介绍

1.1 控制方程

首先引入阶梯函数即Heaviside函数H，规定在一种流体(比如液体)中H为1，在气体中H为0，自由液面则标记H为0.5。为将Heaviside函数与Dirac函数联系起来，将它表示为一维Dirac函数乘积的积分，

例如在二维条件下积分区域是由周线S包围的A(t)面积。

很显然，当(x，y)落在面积A(t)内时，Heaviside函数取值为1而在其他地方取值为0。为了计算Heaviside函数的梯度，设定两类变量，以有上标的变量为第二类变量。首先考虑梯度是针对第一类变量计算，而积分是针对第二类变量计算，这样梯度符号可以和积分号交换顺序。同时Dirac对于两类变量是反互易的，所以梯度符号也可以针对第二类变量计算，即

本模型考虑的是同一种不可压缩流体，密度为常数，则全场密度可用Heaviside函数表示

式中：ρ1和ρ0分别为H取1和0处的密度。而密度的梯度则可以表示为

在此▽ρ=ρ1-ρ0。

流体的运动假定由N-S方程描述，在考虑表面张力的情况下，有

上述方程适用于全流场(含密度场)。β为所涉及问题的维数；κ为自由液面曲率；σ为表面张力系数；δβ是β维的Dirac函数，由一维Dirac函数的乘积构成；n为自由液面的法向量。式(5)最右端积分项的积分区间为自由液面。模型尽管可以考虑流体黏性系数的不连续性，但依然假定流体黏度是常数，即

自由面边界条件获取途径为：在自由液面附近取一小块圆柱体，圆柱体上表面在气相中，而下表面在液相中，分别计算作用在圆柱体表面上的面力及圆柱体内部的体积力。根据达朗贝尔原理，考虑到体积力中惯性力作用，所以各力应该是平衡的。将圆柱体体积缩小并逐渐趋于零，由于体积力是较面积力的高阶小量，故可忽略所有体积力项。将方程在所考察点附近沿自由液面法线方向的平行和垂直方向分解，可得应力张量的垂直分量在自由液面两侧连续，而平行于法线方向的分量则应该满足

式中：[]表示相关变量的阶跃值。速度矢量在自由液面两侧也是连续的。

1.2 N-S方程求解

求解携带外力项的不可压缩流体的动量方程和连续性方程。如果底层网格是结构化正交网格，则采用projection method。首先重构全场密度：

式中：n为时间步数。

动量方程在固定的结构化正交网格上加以求解，而表面张力可以通过变动水面上的网格来计算，所以需要将表面张力扩散至变动水面周围的流体上。因界面由Dirac函数定义，所以这种扩散必然包含多维Dirac函数离散形式的构造。这种构造可从多方面进行，会有不同结果，但其最终形式必须满足一定的要求[12]，如：

为更新自由液面位置，在界面网格上定义真实的界面移动速度矢量。由于自由液面由流体组成且速度场在界面两侧连续，所以液面按照当地流体速度运动，即满足无滑移边界条件。为计算定义在界面网格上的速度，利用多维Dirac函数由流体速度插值得到界面移动速度，即

式中：小写和大写字母分别代表定义在流体结构化正交网格和界面网格上的变量。得到自由液面速度后，其位置更新可按下式进行：

式(11)的离散格式很多。离散格式关乎到浸没边界法的稳定性。采用恰当的离散格式能带来无条件稳定的优良性质。可采用显式向前差分格式。

流体的物性系数，如密度、黏度等，在一种相的流体中取值为常数，而在界面处存在阶跃值，需要在每一时刻利用自由液面位置，来确定全场的密度或黏度分布。

自由液面处密度场的梯度已知，而其他地方梯度值均为零，进而根据边界密度逐步更新全场密度值。但当界面卷曲较严重时，该方法不再采用，可以换用更为一般化的方法。由于式(4)解决了计算全场的密度梯度问题，所以可在其基础上继续取散度得到密度场的拉普拉斯算符：

式(12)右端由传统的中心差分格式计算。求解由此得到的泊松方程，从而得到全场密度。在三维条件下，为加速计算过程，考虑到密度梯度只是在界面附近才取非零值，所以将泊松方程的计算区域局限在界面附近，而边界均采用狄利克雷边界条件。采用格子Boltzmann方法来求解泊松方程，具体过程参考文献[13]。

表面张力计算参考拉普拉斯公式：

式中：f为表面张力面密度；σ为表面张力系数；R1和R2分别为曲面内过曲面上一点的曲线最大和最小曲率，即主曲率。表面张力的方向严格沿曲面上一点的法线方向。在离散网格上利用式(13)，采用三角形网格来离散自由液面：

式中：k(xl)为坐标为xl的点 Nl的平均主曲率的2倍，而定义在点Nl上的面积是泰森多边形的面积：

式中：点 q 是与点 l 共享一个三角形的另外两个节点；α和β为这些节点在三角形内的内角度数。从而法向平均主曲率可按下式计算：

法向平均主曲率乘以表面张力系数即可得到张力的面密度。采用平面三角形来离散曲面，即网格只能通过加密才能逐渐逼近真实的曲面；而在实际中无法还原曲面，也可采用曲面三角形来离散液面，按照以下方法来计算曲面上某一点的平均法向主曲率：

则作用在一小块三角形上的表面张力为：

在自由液面的模拟中，液面与计算区域的边界相交，从上面计算某一拉格朗日点的表面张力来看，该点需要完全被其他点所包围，而这一条件对于处在边界上的拉格朗日点来说是无法满足的。可以在边界外面再新增一排节点，这排节点并不会像内部节点一样参与拉氏力的计算及其扩散和位置的更新，其坐标根据具体所设置的边界条件人为地给定，比如对于滑移固壁、对称边界和周期边界来说，液面在边界处的法向量应该与边界平行。

模拟中速度场的局部差异性较大。随着计算的进行，网格质量开始下降，为使网格能够分布均匀，且三角形网格的边长能够小于预先设定好的上限值，存在下面3种情形需要处理(图1)。

预先设定三角形边长的取值范围，当其边长小于下限时，利用图1(a)所示方法来添加三角形；大于上限时，利用图1(b)所示方法删除三角形；特殊情形需重新划分三角形(图1(c))。上述网格重构可借助C语言中链表容器实现。

图 1 液面三角形网格重构Fig. 1 Reconstruction of liquid surface triangle meshes

2 水箱底部注水液面变动算例

2.1 计算条件说明

为验证上述模型的可行性和准确性，可尝试模拟包含自由液面和热传导的流场。为便于表述，以下物理量均表述为格子单位，时间以计算时步作为单位。计算区域仍然为正立方体，均匀离散成60×60×60的网格，流体运动黏度υ=0.1，热扩散系数D=0.1，自由液面表面张力系数σ=0.7。如图2所示，初始时刻自由液面处于平衡位置，其上部为气体，密度ρg=0.1，其下方为液体，密度ρf=1.0，全场流体处于静止状态，温度统一赋值为T0=1.0。计算区域的顶部为压力边界和温度边界，四周为周期边界，底部为壁面。在底面中央开一半径为10的小孔，从里面以竖向速度u=0.2注入温度为2.0T0的热流体，当计算进行到t=1 100后，将此小孔改为压力出口边界。预期在热流体的注入下，自由液面逐渐上升的同时会伴随一定波动；而当将速度入口改为压力出口后，液体将在内外压差作用下开始泄流，而自由液面则会逐渐下降。

2.2 自由液面模拟结果及分析

为清楚显示自由液面形态，单独将液面与通过计算区域中心且垂直于底面的平面交线绘制在图3中，给出的是几个不同时刻自由液面的位置和形态。

模拟前期，自由液面初始位置靠近底面，注入的液体来不及扩散，分布极不均匀；而在注入液体的推动下液面开始上涨，因而在中央部位向上鼓起。液面曲率越大，表面张力就越大，并且在距离小孔较近的中央部位液面上涨速度较四周快，液面向上鼓起成山包状(图3(a))，最终表面张力和液面中央部位的上下压差平衡，流体在表面张力作用下开始减速，并向四周扩散，液面持续上涨，之前形成的突起逐渐消失(图3(b))。由于立方体四周设置为周期边界，在边界上流速只可能沿竖向，所以速度场在遇到侧边界之前再次转向垂直，进而导致液面周围部分开始上涨。该突起的消失并非一次性单调完成，而是以在液面总体向上平移的基础上叠加一振幅逐渐衰减的波动形式完成。不久突起完全消失，液面基本上成一平面(图3(c))，当液面上涨到一定高度后，底面注入的液体在撞击到液面之前已经扩散得较为充分，在计算区域的水平截面上垂向速度分布得较为均匀，因此，液面不同部位的上涨速度差别不大，所以才会呈现出平面状(图3(c)。液面在停止注入液体的t=1 100时达到最高点(图3(d))，在将速度入口边界修改为压力出口边界后，液体开始自小孔流出，液面逐渐下降(图3(e))，下降前期液面依然成一平面，降至距离出口较近时，液面中央部位开始向下凹陷，形成一漏斗(图3(f))。

为分析液面的波动特征，将测点布置在液面正中央，并统计了该点纵坐标变化情况，结果如图4所示。在注水前期t＜100时，液面中央上升速度非常快，这是突起形成阶段。之后曲线在100＜t＜200时趋于水平，这是由于液面突起达到最高点后向下回弹所致。这段曲线体现该波动叠加在液面总体上涨趋势，液面上升趋势和突起向下的波动叠加在一起相互抵消，所以曲线成水平状。波动过程一直持续，波动幅值同时也按照某种规律衰减。

图 2 初始时刻液面位置和密度分布Fig. 2 Liquid level position and density distribution at initial time

图 3 不同时刻的液面位置Fig. 3 Liquid level position at different times

图 4 液面中央测点纵坐标随时间变化规律Fig. 4 Variation in vertical coordinates with time at central measuring point of liquid level

为定量分析模拟准确性，图5给出了液体所占据的体积随时间的变化。由于在注水期间入口速度均匀分布，流体不可压缩，根据连续性定理，液面以下液体体积的增长速度应该等于入口流量。同理在放水阶段液面下降的速度也应该等于出口的放水流量。从图5中可以看出，虽然在注水阶段液体体积的增长曲线并非严格成一条直线，但与理论直线的差距处于可接受范围。在放水阶段由于给定的是压力边界，出口速度并非均匀分布，所以理论曲线也并非是一条直线，得到的结果与理论解在总体趋势一致的基础上存在一定偏差。究其上述偏差的来源，可能一方面来自自由液面的离散处理，由于采用的是三角形平面面元来离散液面，实际上的液面应该是光滑曲面；另一方面来自自由液面与计算区域边界相交面的处理。前文已经阐述了如何处理液面与物理边界相交的问题，但实际上却无法完全做到。

2.3 温度场模拟结果及分析

为了比较不同扩散系数对模拟结果的影响，在保持其他参数不变条件下改变扩散系数，图6为对应于D=0.005和0.5的温度场。

图 5 液体体积的时间变化规律Fig. 5 Temporal variation law of liquid volume

图 6 扩散系数D=0.005和0.5时不同时刻的温度场云图Fig. 6 Temperature field cloud diagram at different times under diffusion coefficients D=0. 005 and D=0.5

扩散系数表示扩散强度与对流强度的比值，取值越小说明对流作用越明显。比如对于液面上升至顶端的t=1 100时刻，当扩散系数取值较小时，D=0.005，温度场受到对流控制，其分布仅局限在计算区域中央的狭长地带，并不像当D=0.5时分布得比较均匀。为定量分析温度场分布随时间以及扩散系数的变化，图7给出了不同扩散系数下计算区域中心测点温度的变化过程曲线。

温度变化来自于对流以及扩散两方面。液面的存在阻止了注入液体进一步向气体方向侵入，使温度无法仅通过对流尽快从液体区域向气体区域迁移；但是扩散作用则不受流场分布的影响，仅取决于温度场的局部差异。对于扩散系数取值较大的情况(图7)，温度分布主要受到扩散作用的控制，即使自由液面还没有上涨到中心测点，温度场不需要依靠对流机制仍然可以凭借扩散影响到中央测点，所以在图7中可以看到测点温度很快就有所上升。而对于扩散系数取值较小的情况，扩散作用影响有限，温度只能等到液面上涨到测点位置后凭借对流作用逐渐上升，所以测点温度在t=420液面抵达时才开始上升，而此前保持在初始值左右。不管抵达的速度快慢如何，测点温度一旦抵达其最大值附近便不再继续增加，而是在其最大值附近波动，波动原因是自由液面在上升过程中的几何波动。对于泄水过程，孔口由速度入口边界修改为压力出口边界后，因考虑的是不可压缩流体，压力波的传播速度为无穷大，所以立刻就可以观察到降压波。该降压波叠加在之前液面波动所导致的波动之上。对于扩散系数比较小的情况，降压波均出现在温度下降起始点之前，且扩散系数越小，温度下降起始点出现的时间越迟，这是因为温度下降主要受控于扩散作用。在本次模拟中，顶面边界的温度是给定的，并且其取值小于注入液体的温度，所以扩散系数越小，测点温度受到来自于顶面边界条件影响所需要的时间也就越长，而对于扩散系数较大的情形，扩散作用明显，测点温度下降与降压波几乎同时产生。

图 7 不同热扩散系数下区域中心测点温度变化Fig. 7 Temperature change curves of central measuring point of calculating region under different thermal diffusion coefficients