基于SPH方法的二维宽顶堰溢流数值模拟
2017-06-05顾声龙吴玉帅解宏伟袁晓伟
顾声龙 吴玉帅 解宏伟 袁晓伟
摘要：光滑粒子水动力学方法（SPH）是一种基于纯拉格朗日思想的无网格粒子方法，在众多科学和工程领域得到了初步应用。采用SPH方法对二维宽顶堰溢流进行了数值模拟，并利用机群进行并行计算来提高计算效率，将模拟数据与试验数据进行了对比。结果表明，SPH方法可以模拟出二维宽顶堰自由出流和淹沒出流时水流流态、水跌、回流区、水跃现象，计算出溢流流量、流速等水力参数，并且采用数据拟合得出与理论相符的断面流速分布图。结果表明了SPH方法可以较好模拟二维宽顶堰溢流中的流速分布以及水面曲线。
关键词：SPH；二维宽顶堰溢流；数值模拟；自由面计算；并行计算
中图分类号：TV652.1 文献标志码：A 文章编号：1672-1683（2017）00-0127-05
近年来SPH方法在自由表面流动模拟中取得了较大进展，能够模拟传统有网格算法较难处理的自由面翻卷破碎、液舱晃荡、流固耦合等现象，得益于其拉格朗日粒子特性；SPH方法中粒子不仅用于估算场变与近似控制方程，也代表介质系统，具有诸如密度、压力、速度、内能等宏观物理量，相当于物质点。SPH方法既具有了拉格朗日方法描述物质界面准确的优势，又兼备欧拉方法容易计算的长处。因此，SPH法已用来研究渍坝、自由表面流动、多相流等问题。
本文采用SPH方法来模拟二维宽顶堰溢流现象，利用宽顶堰溢流实验实测数据，采用集群并行计算方式，来研究SPH方法在宽顶堰自由出流和淹没出流水力特性方面的适用性及模拟的精度。
1SPH理论
1.1SPH方法的基本方程
SPH方法，本质上是把真正的空间近似划分，最小的划分空间由一系列表面绝对光滑的粒子填充。这些空间粒子常常用坐标、速度、质量、内能等参量表征。推导SPH的基本方程，分为两步：第一步，给出SPH方法的积分表示法，即在整个变量场中，对任意位置应用核近似法进行积分近似；第二步给出其粒子近似值，对给出的函数积分表达式，通过对在其光滑长度范围内的粒子叠加求近似，得到其粒子的近似公式。
1.2SPH方法核近似与粒子近似
在SPH方法中，函数f（x）的积分表示式如下：
（1）
（2）
（3）式中：h为光滑长度，用来定义光滑函数影响区域的光滑长度。由于W不是狄拉克函数，故式（3）的积分表示式只能是近似式，上式即为核近似方程的初始形式。
（4）
（5）
式（5）即为SPH方法的粒子近似法。在粒子i处的函数的粒子近似式最终可写为
（6）
式（6）说明了粒子i处的任一函数值可通过应用光滑函数对其紧支域内所有粒子相对应的函数值进行加权平均近似求得。
1.3控制方程
将描述流体运动的质量守恒方程和动量守恒方程（Navier-Stokes方程）表示为
（7）
（8）
（9）
（10）
2数值模拟
2.1模型建立
本文选取宽顶堰溢流为试验模型，见图1，参数：堰高P=0.08 m，堰宽B=0.15 m，堰厚δ=0.4m，堰上水头为H=0.065 m，δ/H=6。二维计算模型的几何和物理参数与该试验模型相同，由于SPH方法中入流问题仍是一个未解决的难题；故本文入流方式采用“大水箱”，该方法可以在一定范围内保持模型的“恒水位，对计算结果影响不大。计算模型共计20 256个粒子，边界粒子2 125个，本文利用集群进行并行计算；模型参数：核心函数为五次样条核函数；时间步长算法为预估校正法；密度过滤为谢巴德滤波器方法；黏性处理方法为层流和亚粒子尺度紊流；状态方程为微可压；边界类型为边界排斥力型；时间步长为1×10-5s，并且可变；柯朗弗里德里希一列维数为0.2；光滑长度系数为0.92；模拟溢流时长为10 s，集群共运行1.5 h。由于试验模型的二维性较强以及构建三维计算模型的不易之处，所以本文采用二维模型；分析模拟结果时，将试验模型数据转换为单宽数据与模拟结果进行对比即可。
2.2结果及分析
2.2.1水流形态
本文试验模型中hs为下游水位高出堰项距离，H0=H+α0v20/2g；由图2可知，宽顶堰溢流水面有两次跌落；堰坎首端水面跌落是由于水流经过堰坎时，在纵向受到边界的约束，过流断面面积减小，流速增大，势能减小。水面最大跌落处形成收缩断面C-C；而后，由于堰顶阻力，使水面形成壅水曲线，逐渐接近堰顶断面的临界水深。下游水位较低，在堰坎末端再次出现跌落。当下游水位大于堰高后，随着下游水位超过堰项水深幅度的增大，堰项在收缩断面后发生波状水跃，但这时下游水深并不影响收缩断面水深；当hs大于堰项收缩断面水深的跃后共轭水深时，成为淹没出流，下游水深将影响堰顶水位变化，堰上水位被壅高，堰顶呈缓流，见图3。实验研究得出hs/H0=0.8是宽顶堰是否淹没的判别条件，在以下内容中会通过试验数据得出这一比值。
图4（a）-图4（d）四幅图为计算模型1.04 s、2.66 s、3.98 s、9.98 s四个时刻的水流形态。二维宽顶堰溢流模型初始状态中水面线高度为0.145m，在水流达到稳定状态前，堰前水面线的变化可明显体现出宽顶堰对水流的影响，见图4（a）中断面6.92 m处。当计算模型达到稳定状态后，水面线在宽项堰前基本保持平行而且有壅高现象，宽项堰阻碍作用得到了体现，见图4（b）和图4（c）中断面6.8～7.0 m处；堰项水流出现了第一次跌落之后，随后的水面线与堰顶基本保持平行，计算模型在出堰后成功地展示了出堰流的第二次跌落和宽项堰末端有回流区现象的存在，见图4（b）和图4（c）中断面7.2～7.6 m处。这些特性与已知的试验模型宽顶堰堰流特性非常相符。并且从图4（a）-图4（c），可明显看出自由出流的水跃现象；图4（d）呈现出淹没出流时的水流流态，由图4（d）可得，淹没出流的水面曲线变化平缓近乎和堰项平等，在堰出口处，下游水深稍有回升，所以堰下游水位稍高于堰顶水位；在断面7.5 m处还有波状水跃，由于计算模型入流方式为“大水箱”方式，在9.98 s时，水位有一定的下降，故波状水跃现象呈现的并不十分明显。
2.2.2流量及流速
试验模型流量通过直角形三角薄壁量水堰来测量，公式如下：
Q=1.365H2.48 （11）
其中：H=H1-H2。在与本文二维计算模型相同水头下，试验模型数据为：Q=6.5×10-4m3/s：单宽流量Q单=Q/B=4.3×10-3m3/s；下游水位超过堰项水深hs=0.056 m，堰前断面平均流速v0=Q/A=0.19 m/s，故hs/H0=0.84，据此可说明发生了淹没出流。
计算模型流量采用积分求和方式，步骤如下。（1）首先从结果文件中选取计算断面的水平流速，选取方法为：先从结果文件中读取计算断面上水平速度以及竖直坐标，读取依据为水平坐标是否等于设定的计算断面位置；（2）选取数据之后，将水平速度和竖直坐标按竖直坐标进行排序，然后根据公式（12）进行计算，见图（5）。
（12）
流量计算断面选为x=6.75 m，距离堰前端△x=0.25 m，△x/H=3.8；计算结果见图（6）。
二维计算模型中，当t<1.2 s，模型刚开始运行，还未达到稳定状态，其流量变化幅度比较大，即图6中N区；当1.2 s
图7展示了堰前断面6.75 m在t=2.66 s流速变化情况，由图可知，在初始状态中入口处水流速度很稳定，所以从渠底到堰高这段距离内，流速分布比较均匀；从堰高至水流表面，水流速度变化幅度很大。对数据进行拟合，可得变化曲线为三次方程。由于计算粒子数目有限并且一些参数设定带来的计算误差等原因，并没有得到流速分布的对数方程；不过断面上流速的变化趋势与对数方程很相似。图8展示了堰上断面7.2 m在t=2.66 s的流速变化情况，该断面上流速明显比6.75 m断面上流速大，主要是由于过流断面面积变小。
图9展示了在t=2 66 s，从断面6.75 m到断面7.5 m上平均流速变化情况。由图可得，断面6.75 m至断面7.0 m属于堰前，流速变化不大，分布比较均匀，有稍微的增长趋势；断面7.0m至断面7.4 m属于堰上断面，流速快速增长，变化趋势为先增大后减小，最大值在断面7.2 m上；原因是随着堰上断面面积变小，水流动能增加，势能減小，并且断面7.2 m为收缩断面，收缩断面之后水深有一定的增长，水面线与堰项基本保持平行。断面7.4 m至断面7.5 m属于堰后，流速变大，是由于水流势能减小，动能增加。从图7、图8、图9可看出，流速存在一定的震荡现象，主要原因是：（1）计算粒子数目不足，导致结果数据过少；（2）并没有实现恒定流的入流方式，使水头在逐渐渐小。（3）模型参数设置上还未达到最佳，还需进一步调试。
3结语
本文采用SPH方法对宽项堰溢流进行了数值模拟，得出了宽项堰溢流流态、流量、流速并且观察到了水面跌落、回流区、水跃等现象，SPH方法模拟数值结果与试验数据相近以及所得到的撞击、翻卷及碎波的形态、大小和位置都与实验观察吻合，证明了SPH方法在模拟自由水面出流方面的独特优势。由于“入流”问题的存在，对模拟结果带来了一定的误差，在断面流速变化图上可以看到；若模型保持水位恒定，模拟所得的断面流速变化曲线将满足对数方程。本文利用集群进行并行计算，计算效率高，节约了计算时间。
在下一步的研究中，建立三维模型和解决入流问题，使模拟结果与实际更相近以及水流形态与实验观察更吻合，使SPH方法在水动力学方面得到了更广泛的应用。
