基于流固耦合的平板闸门动水关闭过程数值模拟

戴冰清 茅泽育

摘要： 闸门动水关闭是一个复杂的物理过程，在此过程中的水动力特性及闸门本身结构特性易受到多种因素的综合影响。针对平板闸门动态关闭过程，采用动网格技术、RNG kε紊流模型及VOF自由水面处理技术三者相结合的方法建立了三维非定常水气两相流模型，在此基础上，进一步结合流固耦合方法针对不同因素综合影响下闸门的形变、等效应力的变化规律进行数值模拟研究，综合探究来流流量、闸门底缘型式、闸门关闭速度等因素对平板闸门闭门过程的水动力特性及闸门结构特性的影响。研究结果表明：闸门最大形变出现在底缘处，最大等效应力出现在门顶处;合理减小闸门关闭速度及上游来流流量、采用底缘前倾角较大的闸门等措施，有利于改善闸孔出流条件，减小闸门门体所受负压的区域，并减小闸门形变及等效应力。研究成果可为闸门的优化设计及确定合理运行工况提供科学依据。

关 键 词： 平板闸门; 水动力特性; 结构特性; 数值模拟; 流固耦合; 动网格

中图法分类号：  TV135.4

文献标志码：  A

DOI： 10.16232/j.cnki.1001-4179.2021.08.033

0 引 言

闸门是水利枢纽安全运行的关键部分之一，它能否安全运行直接关系着泄水建筑物运行的技术可行性和安全可靠性[1]，也是一直受到关注的问题。闸门动水关闭是一个复杂的物理过程，同时受到上游水头、水流流速、下闸速度、闸门体型和底缘型式等多因素的综合影响。国内外对闸门失事原因的调查结果表明，横向流、折冲水流及水跃冲击等水流会对闸门的结构造成一定的破坏。因此需要综合探究闸门闭门过程的水动力特性及闸门结构特性的变化规律，为工程上闸门的优化设计及确定合理运行工况提供科学依据。

近年來，学者们分别采用物理模型试验、原型观测、数值计算等方法研究了闸门关闭过程中的水动力学特性[2]以及影响水流水力特性的因素。如Willi等[3]采用半经验法研究了溢洪道形状及闸门开度对溢洪道中水流水力特征的影响;Marcou等[4]采用Lattice Boltzmann模型，研究了闸门关闭过程中的明渠水流特性;Masoud[5]认为水闸的出水性能与水流的弗劳德数，和上游水深与闸门开度的比值相关;Marzieyh[6]等采用多普勒测速计测得了闸门附近水流速度的分布情况。南水北调工程建设期间，方神光等[7]模拟计算了同时改变渠首流量、闸门开度和分水口流量时渠道中的非恒定过渡过程;吴保生等[8]根据渠道水位、流量等信息来优化渠道自动化控制系统。还有一些学者采用数值模拟方法进行了初步探索，刘昉等[9]将闸门关闭过程简化为闸门的单向运动，再采用动网格方法进行数值模拟;文林森等[10]基于VOF方法，对某水电站事故闸门闭门过程进行了初步数值模拟研究;郭桂祯等[11]就流体对结构自振特性的影响规律进行了初步研究。

可以看出，前人对闸门动水关闭过程的研究大多仅限于水流的水动力学特性本身，而很少或没有考虑水流对闸门门体本身结构特性的影响。实际工程中，闸门动水关闭过程中水气两相流交替变化复杂，水在流经闸门底缘过程中由于固体边界的影响，必然形成绕流弯曲流场，水流分离现象不可避免造成闸门底部压力梯度变化大，致使整个流域具有复杂的动边界问题。

鉴于此，本文运用FLUENT软件，结合动网格技术，综合考虑了不同入口流量、不同底缘型式、不同闸门关闭速度等因素的综合影响，采用有限元数值模拟计算方法分析平板闸门闭门过程中的水动力学特性变化。同时结合流固耦合方法，对闸门门体的形变、等效应力变化规律进行研究分析。

1 数值模型及理论基础

有限元方法具有计算精度高、能适应较为复杂的几何形状条件等优点，被视为一种具有广阔应用前景的工程分析手段。本文采用平板闸门有限元法展开分析，整个过程可以分为前处理（包括建立模型及划分网格）、求解、后处理（包括采集处理并分析结果）3个部分。

1.1 水流控制方程

采用雷诺时均模型[12]，不可压缩水流的连续性方程和动量方程分别如下：

ρ t +div ρ u  =0 （1）

ρu  t +div ρu u  =- p x + τxx x + τyx y + τzx z +Fx （2）

ρv  t +div ρv u  =- p y + τxy x + τyy y + τzy z +Fy （3）

ρw  t +div ρw u  =- p z + τxz x + τyz y + τzz z +Fz （4）

式中：p是压强，τxx等是因分子黏性作用及湍流黏性影响而产生的作用在微元体表面上的黏性应力τ的分量;Fx、Fy、Fz是微元体上的体积力，本文模型中体积力只有重力，z轴竖直向上，Fx=0、Fy=0、Fz=-ρg。

1.2 紊流模型

采用RNG k-ε紊流模型计算结果具有更高的可信度[13]。其中k方程及ε方程分别为

t  ρk +  xi  ρkui =  xj  αkueff k xj  +Gk-ρε  （5）

t  ρε +  xi  ρεui =  xj  αεueff ε xj  +C*1ε ε k Gk-C2ερ ε2 k   （6）

μt=ρCμ k2 ε

Gk=μt  ui xj + uj xi   ui xj

C*1ε=C1ε- η 1-η/η0  1+βη3   （7）

η= 2Eij·Eij 1/2 k ε  （8）

Eij= 1 2   ui xj + uj xi   （9）

式中：Gk为压力生成项;ueff为有效黏性系数;Eij为时均应变率。

Cμ=0.084 5，C1ε=1.42，C2ε=1.68

，αk=αε=1.393，η0=4.377，β=0.012。其他物理量含义参见文献[14]。

1.3 VOF模型

水气界面的自由表面跟踪采用VOF方法[14]。假设在一个单元中，第q相流体体积分数为αq，则有3种情况：αq=0表示单元网格中无第q相流体;αq=1表示单元网格中全为第q相流体;0<αq<1时，表示通过求解水的容积分数来追踪水气交界面。

1.4 动网格技术

闸门动水关闭过程的复杂性体现在水气两相流交替变化。当水流流经闸门底缘时，由于固体边界的影响，将形成绕流弯曲流场。闸门底缘处产生水流分离现象，从而导致闸门底部压力梯度变大，整个流域具有复杂的动边界问题。对于这类问题，需要动网格模型实现网格划分[15]。

在进行初始模型建立及网格划分时，同时采用边界型函数或者UDF（用户自定义函数），用以定义指定的边界运动方式。将运动函数定义在所划分的局部网格面或网格区域上，如果流场中包含运动与不运动两种区域，则需要将它们组合在初始网格中以进行识别。用于动网格的更新有弹簧光顺模型、动态层模型及局部重划模型。

1.5 流固耦合理论

依据数据传递方式的不同，流固耦合可以分为单向流固耦合和双向流固耦合。由于平板闸门整体刚度较大、位移小，相对于流场中的变化可忽略不计[16]，因此采用单向流固耦合方法进行不同工况下平板闸门变形及等效应力分析。单向流固耦合适用于相互作用中結构位移场变化较小、对流场影响可以忽略的耦合计算，并且仅将流场结果数据传递给固体结构仿真计算。先进行流场计算，然后通过耦合面把流场结果数据传递给结构场作为荷载或边界条件，最后求解结构位移场[17]。

2 模型验证

本文采用加拿大Mica电站进水口闸门水力学试验资料对本文构建的数值模型进行验证[18]。Mica电站进水口采用上游底缘倾角为30°的平板闸门，闸门门宽6.7 m，门高7.188 m，最大门厚1.143 m，启闭机设计容量为248.1 t，闸门事故紧急闭门速度为6.1 m/min，事故闸门上有长为14.5 m的斜喇叭进口，纵向高程差为67 m。电站进水口及引水管道系统的水工模型和事故闸门的几何比尺λl=18，流量比尺λq=1 374.6，压力比尺λp=18，时间比尺λt=4.24。在闸门门体上设置了20个测点。在试验工况H=71.5 m，Q=335 m3/s且闸门关闭速度为6.1 m/min时，试验得到闸门底缘所受压强随闸门相对开度变化的曲线。

图1是考虑流固耦合作用后数值模拟计算结果与模型结果的对比。图1（a）表示闸门底缘测点10所受压强的试验值和数值模拟值，两者变化规律基本一致。在水流阻滞作用下，闸门底缘所受压强快速上升;随着相对开度减小到0.8，水流逐步脱离闸门底缘，底缘所受压强持续减小并在相对开度为0.2～0.3时出现最小值;随着开度继续减小，底缘所受压强逐渐升高，关闭至底部时接近上游静水压强。图1（b）为根据20个测点实测得到的压强值绘制的闸门压强分布图，与图1（c）中模拟闸门压强等值线图对比发现，闸门典型开度下的试验与模拟所得的压强分布规律基本一致，靠近底缘处开始呈逆压梯度分布，下游面板处压强降低为零值附近。验证结果表明本文采用的数值模型可以满足平板闸门闭门过程的水动力特性三维模拟。

3 数值试验的模型建立及边界条件

3.1 计算区域网格划分

将模型简化为进口段、闸室段和出口段组成的三维模型。计算尺寸：宽2 m，高2 m，闸前区域长3 m，闸后区域长15 m。闸门门宽2.5 m，门高2.2 m，最大厚度0.2 m，采用钢制平面闸门型式（见图2）。计算区域以四面体为单元进行初始网格划分，根据接近度和曲率，选择最高关联度，设置局部网格加密。以底缘前倾角为60°型式的闸门为例，初始状态计算区域有限元计算网格总计约43万个，平均网格质量为0.755，如图3所示。

在初始网格划分基础上，首先应用ICEM CFD软件进行初始模型建立及网格划分，同时采用UDF（用户自定义函数）定义指定的边界运动方式为闸门以不同速度匀速下落，并且将运动函数定义在所划分的局部网格面或网格区域上。

本文采用动网格中的动态分层法及局部重划法，以适应闸门几何和网格结构的复杂性[19]。其中，动态分层法的中心思想是：在边界发生运动时，如果紧邻边界的网格层高度增加到一定程度，就将其划分为2个网格层;如果网格层高度降低到一定程度，就将紧邻边界的2个网格层合并为1层。而局部重划法，则是用于将计算过程中畸变率或尺寸变化过大的网格集中在一起进行局部网格的重新划分。若重划后的网格可以满足畸变率及尺寸要求，则用新网格代替原来的网格，否则放弃重新划分的结果[20]。将这两种方法进行结合，可以更大程度地保证网格质量及计算精度，使闸门区网格更新后不发生明显的畸变现象。

3.2 数值计算方法

设置计算域左端进口条件为压力进口，水深为1.5 m;出口为压力出口，参考大气压位置设置在自由水面处，并设置参考操作工质密度为1.225 kg/m3。流场与闸门界面设置单向流固耦合边界条件，壁面默认为无滑移边界条件。采用控制变量法进行工况计算方案设计，计算工况如表1所列。

计算方法：选择基于压力求解器的瞬态计算方法;选择用于非稳态可压缩或不可压缩流体流场中求解压力速度耦合关系的PISO算法;选择PRESTO！格式空间离散化压力方程;采用二阶迎风格式动量方程。其余参数保持默认值，并设置时间步长为0.005 s，计算1 000步，瞬态场每一时间步的最大迭代次数为10。先将静态闸门全开时的状态按上述设置进行计算，直到整个自由水面以下的计算流域充满水（见图4），再将最后时刻的case和data值导入动网格计算case中作为初始化数据，开始计算。

3.3 流固耦合计算

闸门材料为Q345C钢，其弹性模量为2.1×105 MPa，泊松比为0.3，密度为7.85×103 kg/m3。由材料力学知识可知，由于金属材料Q345C钢为延性材料，通常以屈服形式失效，符合米塞斯屈服准则，一般采用第四强度理论（又称形状改变比能理论，不论什么应力状态，只要形状改变比能达到极限值，便引起屈服）进行校核。根据第四强度理论[20]，等效应力计算公式为

σ = 1  2    σx-σy 2+ σy-σz 2+ σz-σx 2+6 τ2xy+τ2yz+τ2zx     （10）

实际工程中闸门通过启闭机启闭，仅在空间z轴上有位移，因此将闸门上端设置为固定端。将平底底缘及前倾角底缘型式闸门以四面体为单元划分网格。在设置的耦合面上加载来自流体的作用力，求解闸门总变形及等效应力。

4  平板闸门关闭过程的数值模拟及水动力特性分析

将闸门上端设置为固定端，底缘为自由端，闸门可视为悬臂梁模型，分析研究流固耦合作用下闸门关闭过程中的水动力特性。

4.1 结果计算

对于闸底底坎为平顶堰，当相对开度e/H≤0.65时，为闸孔出流;e/H>0.65时，为堰流。对于底坎为平顶堰的闸孔出流，有以下计算公式[21]：

Qv=μbe 2gH0  （11）

式中：e为闸门开度;H为闸门底坎顶到自由表面深度;μ为闸孔自由出流流量系数，μ=0.60-0.76 e H ;b为闸孔宽度;H0为闸孔全水头，H0=H+α0v0/（2g）。质量流量与体积流量的换算，Qm=ρQv。

理论计算值与数值模拟（工况6）得到的闸孔出流质量流量值随闸门开度的变化曲线如图5所示。显然，计算值与数值模拟值误差较小。

4.2 不同工况下的压强分布规律

本节控制闸门关闭速度为0.3 m/s，探究其他因素对闸门关闭过程的影响规律。图6是平底底缘与前倾角60°底缘的闸门在相对开度分别为0.65和0.30的情况下计算区域内的压强分布。显然，闸门底缘处，压强呈大梯度逆压梯度分布。开度较大时闸后的回流现象并不明显，形成的漩涡对闸门底缘的压强分布影响较小，底缘没有形成明显的负压区;当开度较小时，水流与闸底分离产生脱流现象，闸门后水流漩涡区明显增大。这一区域处于水气两相流混合交替状态，闸门底缘产生射流扰动作用，需要补入气体，吸气作用强烈，造成图示闸门底缘处的负压分布，从而使闸门底缘压强呈明显的逆梯度分布特征。

由图6还可以看到：当Q=30 m3/s，v=0.3 m/s時，动水压强峰值出现在闸门上端，负压出现在闸门底缘。平底底缘型式闸门的关闭过程中，较小开度时水流分离点一般出现在上游，在闸门下部底缘产生空腔时无法及时从下游补充空气，闸门底缘产生极不稳定的负压，工程中此时的动水垂直力会呈现下吸力特性。而前倾角底缘型式的闸门设计更为合理，流线型轮廓的水流条件较好。闸门底缘的水流分离现象不明显，几乎没有形成负压区域，工程上可以减少底缘的空化现象。前倾角底缘闸门所受到的压强峰值较小，且底缘倾角为60°及45°时，在闸门处基本不产生负压。因此，可以通过改变闸门底缘型式以改善闸孔出流条件，减小过流对闸门产生的影响，并且宜选用较大倾角。

由图7（a）和图8（a）可见，最大压强随开度减小而增大。来流流量分别为18，30，36 m3/s时，最大负压分别为29.80，61.60，83.16 kPa。

由图7（b）和图8（b）可见，闸门关闭前期，所受到的负压处于较为稳定的小范围不规律变动状态。原因主要有：① 在来流流量不变的情况下，随着闸门的关闭，过流断面面积减小，流速增大，弗劳德数增大，水流本身处于强烈的紊动状态，各个物理量及补气条件受影响，则闸门底缘所受负压产生不规律变动;② 由于计算机内存及计算时间限制，网格的数量及质量对计算精度产生了一些影响;③ 在模型建立时，对边界条件及初始条件的简化在一定程度上也影响着计算精度。相对开度减小到0.4时，闸门所受最大负压（图中的最小压强）急剧增大。

来流流量越大，闸门受到的负压也越大，负压峰值分别为-8.7，-16.9，-23.4 kPa。原因为：来流流量越大，断面流速越大，在闸孔处断面收缩时水流分离形成的涡区较大，水流分离现象越剧烈，而闸后补气越困难，则造成底缘处承受负压越大。

4.3 不同闸门关闭速度对闸门的影响

表2是在3个典型开度下，以平底底缘闸为研究对象，采用不同关闭速度时闸门所受的最大及最小压强值。同一开度下，关闭速度为0.4 m/s时，闸门所受压强最大，0.2 m/s时最小。开度越大，关闭速度对压强大小影响越大。各个关闭速度下的闸门所受负压峰值相差不大。因此减小闸门的关闭速度，可以减小闸门承受的最大压强，但对闸门底缘承受的负压影响不大。

5 平板闸门结构特性数值模拟及分析

5.1 平板闸门变形情况分析

由图9～10可知：当Q=30 m3/s，v=0.3 m/s，e/H=0.2时，从顶端至底部闸门总变形逐渐增大，在底缘处结果偏离初始状态距离最远。闸门的最大形变量随开度减小而增大。开度相同时，来流流量越大，闸门最大形变量越大。例如平底底缘闸门在开度为0.2时，3种工况下闸门的最大形变值分别达到了3.25×10-4，7.96×10-4，10.96×10-4m。

同一工况下，前倾角底缘型式闸门的最大变形区域及最大变形量均小于平底底缘型式闸门，说明前倾角底缘型式闸门具有更好的结构稳定性，且倾角为60°时，闸门产生的变形量相对较小。

5.2 平板闸门等效应力分析

由图11～12可见：当Q=30 m3/s，v=0.3 m/s，e/H=0.2时，随着开度的减小，闸门等效应力值呈增大趋势。主要原因有：① 动水压力随开度的减小逐渐增大;② 动水压力对闸门的作用面积逐渐增大。最大等效应力值出现在闸门顶部，并且往闸门底部逐渐减小。在相同开度时，来流流量越大，闸门最大等效应力值越大。例如，在开度0.2时，3种工况下，平底底缘闸门的最大等效应力值分别达到了5.85，14.93，20.68 MPa。

60°前倾角底缘型式闸门的最大等效应力值始终小于平底底缘型式闸门，更不容易被破坏。这与前文中前倾角底缘型式闸门的水流条件较好、不易形成负压区域等优点相结合，在工程中建议使用较大前倾角底缘型式闸门而非平底式。

6 结 论

闸门动水关闭过程中涉及到复杂的水气两相流，整个计算域具有复杂的动边界问题。为了研究平板闸门闭门过程的水动力特性及闸门结构特性，本文采用动网格技术，应用RNG kε紊流模型、VOF方法及流固耦合方法，建立了平板闸门动水关闭过程的三维非定常水气两相流模型并应用试验资料进行了验证，进一步进行了数值模拟，主要研究结果如下。

（1） 闸后出现的漩涡区可以减小水雾对闸门的侵蚀作用。在闸门底缘处，呈大压强梯度变化并呈逆压梯度分布现象。

（2） 闸门最大形变出现在底缘处，最大等效应力出现在门顶处。因此工程上可以针对闸门的不同部位采取措施，进一步提高闸门的性能。

（3） 合理减小闸门的关闭速度可以减小闸门承受的最大压强，但对闸门底缘承受的负压无明显影响。入口流量对流域及闸门也有一定影响，入口流量越大，闸门承受更大的负压值，产生较大形变及等效应力值。因此也需要合理控制上游来流流量使闸门安全运行。

（4） 前倾角底缘型式闸门具有水流条件较好、不易形成负压区域、闸门形变及等效应力较小等优点，在工程中建议使用较大前倾角（如60°）底缘型式闸门而非平底式。可改善闸孔出流条件，减小过流对闸门产生的影响。

参考文献：

[1]  黄金林.平面闸门底缘型式及选择[J].中国水利，2004（14）：44-45.

[2] 赵梦丽.泄洪洞事故闸门动水闭门水力及爬振特性研究[D].天津：天津大学，2017.

[3] WILLI H H，ROGER B.Plane gate on standard spillway[J].Journal of Hydraulic Engineering，1988，114（11）：1390-1397.

[4] MARCOU O，CHOPARD B，S EL YACOUBI，et al.Lattice boltzmann model for the simulation of flows in open channels with application to flows in a submerged sluice gate[J].Journal of Irrigation and Drainage Engineering，2010，136（12）：809-822.

[5] MASOUD G.Flow through side sluice gate[J].Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE），2003，129（6）：458-463.

[6] MARZIEYH E，MANOUCHEHR H，SAYED SAEID E.Flow characteristics of a sharp-crested side sluice gate[J].Journal of Irrigation and Drainage Engineering，2014，141（7）：6014007.

[7] 方神光，李玉榮，吴保生.大型输水渠道闸前常水位的研究[J].水科学进展，2008（1）：68-71.

[8] 吴保生，尚毅梓，崔兴华，等.渠道自动化控制系统及其运行设计[J].水科学进展，2008（5）：746-755.

[9] 刘昉，赵梦丽，冷东升，等.不同底缘形式的平板闸门水力特性数值模拟[J].水利水电科技进展，2017，37（5）：46-50，77.

[10]  文林森，王才欢，杨伟，等.水工附环闸门闭门过程水力特性数值模拟研究[J].长江科学院院报，2017，34（10）：68-73.

[11] 郭桂祯.平板闸门垂向流激振动特性与数值计算研究[D].天津：天津大学，2011.

[12] 王福军.计算流体动力学分析[M].北京：清华大学出版社，2004.

[13] 胡揭玄.横风下公路桥梁与车辆系统的气动特性的数值模拟及试验研究[D].长沙：长沙理工大学，2012.

[14] 王福军.计算流体动力学分析—CFD软件原理与应用[M].北京：清华大学出版社，2004.

[15] 郭泽宇.基于自适应动网格技术的类鱼游动水动力及流场特性的数值研究[C]∥中国力学学会庆祝中国力学学会成立60周年大会论文集，2017.

[16] 李明.弧形闸门动力特性及流激振动数值模拟[D].长沙：长沙理工大学，2013.

[17] 闻德荪，李兆年，黄正华.工程流体力学（水力学）[M].北京：高等教育出版社，2004.

[18] 中国水利水电科学研究院.Mica电站进水口闸门水力学试验研究[R].北京：中国水利水电科学研究院，2011.

[19] 张斌，杨涛，丰志伟，等.网格质量反馈的弹性体动网格改进[J].国防科技大学学报，2018，40（1）：10-16.

[20] 甘慶明，付芳琴，张琴，等.基于动网格的往复式压缩机进气阀数值模拟[J].石油机械，2019，47（9）：105-110，117.

[21] 吴持恭.水力学[M].北京：高等教育出版社，2005.

（编辑：胡旭东）

引用本文：

戴冰清，茅泽育.基于流固耦合的平板闸门动水关闭过程数值模拟

[J].人民长江，2021，52（8）：214-221.

Numerical simulation of plane gate closing under running water based on

fluid-structure interaction

DAI Bingqing，MAO Zeyu

（ Department of Hydraulic Engineering，Tsinghua University，Beijing 100084，China ）

Abstract：

The gate closing under running water is a complex physical process.In this process，the hydrodynamic characteristics and the structural characteristics of the gate itself are affected by many factors.Aiming at the dynamic closing process of the plane gate，a 3D unsteady water-gas two-phase flow model was established by combining dynamic mesh technology，turbulence k-ε model and VOF free surface treatment technology.On this basis，combined with the method of fluid-structure interaction，the variation law of gate deformation and equivalent stress under the influence of different factors was numerically simulated，and the influence of inflow flow，gate bottom edge type，gate closing speed and other factors on the hydrodynamic characteristics and gate structure characteristics of plane gate closing process were comprehensively explored.The results showed that the maximum deformation of the gate occurred at the bottom edge，and the maximum equivalent stress occurred at the top of the gate.Reasonable reduction of gate closing speed and upstream inflow flow，and the gates with large bottom rake angle were conducive to improving the outflow conditions downstream，reducing the negative pressure area of gate body，and reducing the gate deformation and equivalent stress.The above research results can provide a scientific basis for the optimal design of the gate and the determination of reasonable operating conditions.

Key words：

plane gate;hydrodynamic characteristics;structural characteristics;numerical simulation;fluid-structure interaction;dynamic mesh