基于流固耦合的闸下淹没出流流场瞬态分析

赵 勇，郑圣义，蔡卫江，荣 红，许 栋

(1.南瑞集团(国网电力科学研究院)有限公司，江苏南京211106；2.河海大学能源与电气学院，江苏南京211100)

水工钢闸门的安全关系到人民生命、财产安全，而水流诱发的闸门剧烈振动很容易造成闸门破坏[1-2]。闸门开启泄流时，由于周围流场运动因素不稳定，闸门上出现脉动荷载，从而引发闸门振动。据不完全统计[3]，我国失事的闸门中有20多起是因水流诱发振动，流激振动应予以高度重视。不同水流特性和结构特性情况下，闸门振动的性质大不一样。特定条件下，工作闸门若是局部开启，闸后将发生淹没水跃[4-5]。淹没水跃循环往复地冲击门叶结构，旋滚涡流容易造成结构空蚀以及振动，恶劣情况下闸门会因动力失稳而破坏。通过对国内外文献的查阅总结[6- 8]，应用有限元数值理论，构建流固耦合有限元模型，分析比较各工况下弧形闸门结构动特性，总结结构动特性在各工况下的变化规律，获得对结构振动影响较大的开度区间。对于水工钢闸门设计有一定的参考价值。

现阶段国内对弧形闸门淹没出流下瞬态流速变化、动水压力变化和湍流动能变化等方面的研究存在着大量空白。本文首次基于ANSYS FLOTRAN多场耦合分析平台，对水体—弧形闸门耦合模型进行了瞬态求解，研究了不同开度下流场瞬态流速、近壁水体动压力和湍流动能随时间变化变化规律。

1 水体—闸门耦合理论基础与基本方程

水—闸门耦合问题属于流固耦合问题流体动力作用下有限位移的长期问题，在泄流过程中，水体与闸门面板间有明确的耦合界面，水体位移比较有限，且泄流过程时间较长。水—闸门耦合(流固耦合)基本方程包括流体控制方程、耦合界面节点动压力方程、等效节点载荷方程和结构控制方程[9]等。

(1)流体控制方程

(1)

(2)耦合界面节点动压力方程

(2)

式中，Np=f+g+h，其中，f，g，h为三部分节点动压力矩阵的行数；{Qp}为Sp面处的节点动压力；{Qn}为耦合界面的节点动压力；{Qu}为其他存在节点上的动压力。推导得耦合界面Sn上的节点动压力的基本方程

(3)

(3)等效节点载荷方程。耦合界面上的节点动压力{Qn}可分解出耦合界面Snk上第k个单元S个节点处的动压力，即

(4)

式中，{q}(k)为第k个单元s个节点处的动压力；[A](k)为位号矩阵；SE为耦合界面上的流体单元个数。附加动压力(结构坐标系下)的等效节点载荷为

(5)

(4)结构控制方程

(6)

2 模型建立

2.1 工程实例

本文以黄河某水电站胸墙式泄洪闸水体—弧形闸门为研究对象。针对该工作门开启泄流情况，应用数值模拟软件ANSYS，建立水体—弧门耦合有限元模型，研究淹没出流下流场的瞬态特性。某胸墙式闸门为深孔闸门，正常水位下门后有一定深度的水体。水库正常蓄水位134.00 m，下游平均水位125.30 m，高于溢流堰顶，水头差8.70 m。闸门底槛高程120.192 m，闸门运行的水力特点是低水头闸下淹没出流，要求操作方式为动水启闭。实际运行中，闸门常局部开启泄流，闸门泄流时门底后出现淹没水跃，振动问题较为严重。

2.2 有限元模型及边界条件

计算选取上游段水体为孔口正前方水体，赋予其自由液面以下的相应水压力，不影响计算结果。整个流体区域长20.6 m，宽8.904 m；门前水体长3.5 m，高4.155 m；门后水体长17.1 m，高11.420 m。工作闸门尺寸8.904 m×5.977 m，面板前沿距离支铰中心9.0 m，支铰中心高程126.0 m，门底高程120.192 m，闸门设计水头13.808 m。坐标原点取在两支铰中心线的对称中心，x轴为水流方向，y轴沿铅直方向，z轴沿主横梁轴向[10]。工程实例的流体为水，不考虑温度，只需给出水的密度和粘度。钢闸门为Q235钢，弹性模量E=2.06×105MPa，泊松比μ=0.3，容重γ=78.5 kN/m3。弧形闸门是典型的空间薄壁复杂结构体系，本文首次采用实体单元模拟面板。水体以及弧形闸门各构件模拟单元的选择为：流体(水体)为FLUID142单元，面板、主梁、纵梁、边梁、小横梁、支臂、筋板为等SHELL63单元；支臂斜撑杆为BEAM188单元；面板、圆柱铰、支座与支铰间的连接块为SOLID187单元。网格采取不同区域不同网格密度，在流固耦合交界面附近区域的水体网格较密，远离交界面区域的水体网格相对稀疏。水体—弧形闸门耦合有限元模型如图1所示。

图1 水体-弧形闸门耦合有限元模型(开度0.20)

流体边界条件：假设底面和侧面为无滑移固壁，即x、y、z向流动速度均为零；同时考虑计算简便，上下游上体上表面设置为固壁边界并施加滑移条件，确保有水流方向的速度。上游进口为压力入口，出口为自由出流。

结构边界条件：闸门在支铰孔和吊耳孔处受x、y、z方向位移约束，x轴沿水平方向，y轴沿竖直方向，z轴沿支铰中心连线方向。本文选择RNGk-ε模型模流态[11]，相比于标准k-ε模型，RNGk-ε模型在ε方程中增添了一项用以反映主流时均应变率，从而对耗散率方程的模拟实现了优化，湍流的各向异性得以表达，非常适合用于模拟流线的弯曲流动。

3 水体—弧形闸门耦合流场瞬态分析

3.1 瞬态流速和动水压力分析

图2所示为0.20开度下1.00 s时刻对应的流场三维流速矢量。从图2可以看出该开度下最大流速出现在闸门的底缘附近，过流断面中心处流速相对左右侧壁稍大一些。为便于观察，选取各流场对称面流速矢量，并将闸门底缘附近结果进行适当放大。部分开度下瞬态流速矢量如图3、4、5、6所示。

图2 0.2开度下1.00 s流场流速矢量(单位：m/s)

通过对各开度下流场瞬态流速和动水压力的计算分析可知：各开度下水流在流经弧形闸门底缘时，水流的流态发生急剧变化，由于过流断面的缩减，水流流速增加，且流经门底时流向发生偏转，门底后方水域均出现了逆向回流。在计算时间内，各开度下门底附近水流的最大流速随时间有不同幅度的增加，在0～0.16 s时间内，流速急剧增加，开度为0.05～0.20最为明显，0.05时瞬态流速相对最高，呈现近乎射流状态；0.05～0.20时流速变化剧烈，流态相对恶劣。

从图3、4、5、6可以看出，近壁动水压力大小分布不再满足“梯形递增”的形式，闸门底缘处水体压力不是最大的，门体附近最大水压力集中在距离闸门底缘约1/3处。闸门底缘前后断面流体运动规律符合水力学伯努利能量方程，即流速高处压力低，流速低处压力高。在计算时间内，各开度下门前水体(除底缘附近水体)压力随时间非线性地递增，开度越小门前水体压力增加越多，且离上游耦合面下1/3处越近增加越明显，最大增加约0.3 MPa，增幅为35.1%(开度0.05)。在计算时间内，各开度下门后水体压力随时间非线性地递减，随着开度的递增，门后水体压力减小的幅度逐渐变小。其中，开度为0.05时水体压力最大降幅为112.4%，而开度为0.75时水压最大降幅仅为24.8%。反映了开度越小，门后水体压力下降越明显。

图3 0.05开度下瞬态流速矢量(单位：m/s)

图4 0.05开度下动水压力云图(单位：Pa)

图5 0.375开度下瞬态流速矢量(单位：m/s)

图6 0.375开度下动水压力云图(单位：Pa)

3.2 闸门底缘近壁水体动压力分析

为了解闸门底缘处水体压力瞬态变化情况，在闸门底缘处选取4个近壁水体节点，生成节点动压力时程变化曲线。门底前缘和门底后缘各取2个节点，底缘处近壁节点选取情况如图7所示。图8为各开度下4个节点动压力时程变化曲线。

图7 闸门底缘处近壁水体节点选取

从图7的节点动压力时程变化可以看出：①各开度下不同时刻，水体节点压力不同。固定开度下，压力大小关系为节点1>节点2>节点3>节点4，说明在门底前缘上方压力高于下方压力，在门底后缘上方压力低于下方压力，门底前缘压力要高于后缘压力。②对比各开度下各节点压力峰值，发现其峰值均随开度的增加而增加，从0.05开度到0.75开度，各节点的压力峰值的增幅依次为70.8%、88.1%、117.1%和131.3%，说明闸下泄流时过流面积越小流速越高水压越低，过流面积越大流速越小水压越高，符合伯努利方程；同时反映出随着开度的增加，闸门底缘近壁水体中，门底后方水体压力变化最显著。③开度0.05～0.20时，门底后缘节点3和节点4压力先后出现负值，其大小在-6 kPa左右，说明该开度区间下门底产生负压。当开度大于0.375时负压不再出现。反映了小开度下淹没出流造成的闸门底缘涡流空化现象，使得闸门底缘存在着下吸力。④在计算时间内，各节点压力值变化方式不同。在0～0.2 s时间内，开度为0.05～0.20时，各节点压力值出现“陡增陡减”变化；而0.375～0.75开度下，各节点压力值出现“陡增缓减”的变化，反映出开度在0.05～0.20区间，作用在闸门底缘的脉动压力变化较剧烈。

图8 各开度下门底近壁节点压力时程变化

4 结 论

通过对不同闸门开度下的流场流速和近壁水体动压力的瞬态计算分析，得到以下结论：①开度在0.05～0.20区间，门底流速急剧增加流态相对剧烈；②闸门底缘处水压力并非最大，门体附近最大水压力集中在上游面板距离底缘约1/3处；③门后淹没水跃对闸门产生较强的脉动压力，开度在0.05～0.20区间，门底后缘水体节点压力出现负值，表示闸门底缘存在着下吸力，其大小在-6 kPa左右；④开度在0.05～0.20区间下，各节点压力值出现“陡增陡减”的变化，反映作用在闸门底缘的脉动压力变化比较剧烈。