考虑流固耦合效应的弧形钢闸门地震动力响应研究

刘 丰, 王正中, 厉宽中， 徐 超, 仵 凡, 张欢龙, 张学东

(1.西北农林科技大学 水利与建筑工程学院， 陕西 杨凌 712100； 2.中国电建集团华东勘测设计研究院有限公司， 浙江 杭州 311122； 3.黄河勘测规划设计研究院有限公司， 河南 郑州 450003)

1 研究背景

近年来我国兴建了一大批高坝水库，金属结构制造水平逐步提升，弧形钢闸门朝着高水头、大孔口的方向发展，但作为水工枢纽的调节咽喉，考虑动力响应分析的结构设计理论方法仍亟待完善[1]。西南地区属于地震高发带，在汶川地震中，许多闸门结构系统遭到严重破坏[2]。闸门结构破坏和变形不仅影响闸门正常使用功能，严重时还会导致大坝丧失调蓄水流功能，进而引发溃坝等次生灾害发生，为此，亟需对弧形钢闸门在地震作用下的动力响应规律进行研究分析，建立弧形钢闸门在地震作用下的安全评价准则。

弧形钢闸门动力响应分析一般可采取实地原型监测、等效模型试验、仿真数值模拟等方法。原型监测需要实地考察研究，无法对闸门结构在地震等灾害下的安全性进行预测；模型试验方面，其制作结构复杂、难度高、代价大，周期长且尺寸效应明显，难以满足结构分析要求；有限单元法等数值模拟方法具有可操作性强、便捷的优点，然而精确化精细化数值模拟、大量的计算设备资源及时间成本是该方法面临的主要问题[3-4]。随着计算机在数值计算速率瓶颈上的突破，数值模拟方法将成为工程应用中的主要设计分析方法。当前已有不少学者采用数值方法对该问题进行了研究，其中，孔剑等[5]和李坤等[6]基于ANSYS软件采用了振型分解反应谱法对闸门结构抗震效应及自振特性开展了研究，研究结果表明该方法能够反映闸门在地震荷载作用下的真实受力状态，该闸门结构满足抗震要求。吴一红等[7]以伽辽金方法导出流体与结构的系统有限元方程，对一弧形闸门动力特性展开分析讨论。Brusewicz等[8]和李桑军等[9]对弧形钢闸门进行干湿模态分析，认为在对闸门动力特性分析时有必要考虑流体影响，动力特性分析对闸门安全评价有重要意义[11-12]。Faridmehr等[10]在ABAQUS/Explicit中建立三维弧形闸门模型，考虑流固耦合作用对结构动力响应展开分析，发现与静力分析结果差距很大，闸门结构分析时流固耦合作用不可忽视，利用流固耦合数值方法对闸门结构动力响应进行研究有重要意义。Fenves等[13]和Gogoi等[14]建立坝-水体耦合体系，对耦合体系材料特性展开参数分析，结果表明材料属性和水体可压缩性对耦合系统有重要影响。Buldgen等[15]结合数值模拟与模型试验对经典Westergaard地震动水压力求解公式展开研究，结果表明由于该公式假定结构完全刚性而忽略与流体域之间的耦合作用[16]，附加质量法求解结果比考虑流固耦合数值解更为保守。针对弧形闸门在地震动荷载作用下造成的结构破坏，现有理论分析方法难以反映真实振动耦合工况[17-18]。当前，闸门结构动力设计还没有统一的设计标准[19]，需要借助各种科学技术手段对其展开研究。上述研究多关注于闸门结构计算，然而，对闸门在地震过程中与流体作用下动应力、位移等参数的分析均未涉及，均侧重于结构静态响应分析，但是地震荷载激励下作用在闸门迎水面的水体参与闸门运动的耦合过程不可忽略[20-21]，闸门结构的振动运动反作用于水体，使其流场发生改变，流场改变后的水体又影响闸门结构的阻尼力、弹性力、惯性力等动力特性，从而影响闸门结构的动力响应特性，但现有研究中缺少考虑地震激励引发动水荷载作用下流体与固体闸门结构间闸门动力响应的分析[22-23]。

为此，本文建立了闸门-水体三维耦合数值模型，在地震荷载激励下考虑流固耦合效应，对地震动荷载修正，选用具有代表性的重力坝-库水模型验证了本文计算方法的有效性和准确性，从动力方面入手，对弧形钢闸门在地震激励下位移和应力的动力响应以及时空分布特征进行分析研究，对弧门地震动水压力计算分析，并与规范设计结果进行对比，以期为地震激励下弧形闸门结构设计提供指导。

2 资料与研究方法

2.1 弧形闸门概况

选用某工程表孔弧形钢闸门进行考虑流固耦合效应的地震动力响应分析，该闸门半径为21 m，支铰距闸门底垂直距离为8.97 m，面板宽为12 m。弧形闸门是主要由门叶和支臂组成的空间薄壁结构体系，采用有限元分析软件ANSYS建立弧形钢闸门三维有限元模型，对弧形钢闸门在设计水位下的地震响应进行分析。门叶结构采用能反映空间应力的壳单元shell93，材料为Q345B；支臂截面为箱型截面，采用壳单元shell181，钢材型号为Q235B，有限单元网格划分单元总数14 922个。图1为闸门有限元模型图，图中沿水流方向为X轴，闸门中轴线为Y轴，垂直水流方向为Z轴(主横梁轴向)，其正方向按右手坐标系确定。在支铰处只释放Z方向的转动约束，面板底缘与门槽接触处施加Y方向位移约束，施加Z向的位移约束于面板两侧，设计工况运行水位为16.6 m。

图1 弧形钢闸门有限元模型

2.2 研究方法

闸门在运行过程中受到过闸水流动荷载，需要考虑水体域的可压缩性对流场的影响以及对闸门动力响应的影响，虽然针对实际求解时讨论的三元流动问题，为了简化求解，对流体假设为不可压缩的理想流体，但仍然在大多数的情况下不能得到解析解，只能采用数值求解的方法。闸门流固耦合的数值计算通过流体域控制方程提供压力，固体结构依据压力计算得到变形位移，流固耦合方程在二者之间完成计算数值交换[24]。

2.2.1 流固耦合有限元控制方程 在流固耦合有限元分析方法中，对流体控制方程假设流体为均匀、无粘性且无漩的理想流体，同时还仅限于对线性小变形的情况[25-27]，其流体有限元采用压力表达的控制方程为式(1)：

(1)

式中：u、v、w为各自方向的位移，m；ρ为流体密度，kg/m3；p为动压力，kPa。虽然可以列出物理控制表达方程，但是对于三维问题的微分方程的求解十分困难，采用Galerkin对整个流场域求解也是很困难，尤其是在流固耦合交界面求解，所以需要把流体域划分为流体元，然后集总整个流体域离散的运动方程，将每个流体元的求解结果Gauss数值积分加权集总，得到控制方程系数矩阵，进而求解得到整个流体域的运动方程(2)：

(2)

式中：H、A、E、B分别为矩阵函数系数；r为位移量，m；q0为激励矢量，N。同样也需要对固体结构进行有限元离散，固体结构离散后的运动方程见式(3)：

(3)

式中：Ms、Cs、Ks分别为结构的质量矩阵、阻尼矩阵、刚度矩阵；fp为交界面处的流体作用力,N；fo为除fp外的其他外界激励矢量，N。公式(3)代入各流体元贡献集总fp得到公式(4)与流体接触结构运动方程：

(4)

公式(2)与(4)即为流固耦合系统的运动控制方程。对于时域求解，上式可以归并为公式(5)：

(5)

2.2.2 考虑流固耦合效应的动力响应分析方法 流固耦合计算分析需要首先在几何模型软件中生成结构模型和流场模型，利用ANSYS APDL编程建立弧形闸门力学模型，在Workbench中快速生成外部流场，导入计算平台生成有限元模型。然后设置闸门与水体接触面为流固耦合交界面，随后输入外部激励荷载，以驱动模型耦合迭代。FLUENT中流体域在外部激励的作用下产生压力传递给固体，ANSYS中固体将受到的压力以边界形变反馈给流体域，如此往复迭代直至达到计算收敛条件为一个计算子步，此间动网格技术对流体域网格更新修正以保证模型网格不发生畸变。经完成设定计算时长，固体、流体计算模块输出计算结果，流固耦合系统迭代计算流程如图2。

图2 弧形钢闸门流固耦合系统迭代计算流程

3 流固耦合模型动水压力验证

3.1 地震输入动荷载修正

地震动荷载输入有加速度输入模型和位移输入模型，其中位移输入模型更适用于动力响应空间效应明显的结构，而弧形钢闸门正是多点支撑的空间薄壁结构，因此，首先需要将地震加速度时程曲线积分转换成位移曲线，然后再加载至闸门结构。图3为地震加速度历时曲线，地震响应分析场地为Ⅱ类场地，地震波类型为EI波。地震波作用时间子步长0.02 s，作用时间总共10 s，设计最大地震加速度为0.2g，对应抗震设防烈度为8度，地震作用荷载方向为顺水流方向。地震波在实际采集中，受到各种人为误差的影响会造成通过加速度积分得到的位移曲线出现漂移的现象，所以其位移时程曲线需要先经过基线法校正[28]。图4为经校正后的地震波位移时程曲线，由图4可知，整个历时中位移数值围绕在基准值0上下波动，没有漂移现象出现。

图3 地震波加速度时程曲线 图4 地震波位移时程修正曲线

3.2 地震动水压力验证

有限单元法较之有限体积法在数值模拟计算中得益于连续性和小变形使其计算精度更高，计算结果也更为可信，流固耦合计算涉及流体计算和固体结构计算，且固体结构计算结果受到流域计算结果的影响，两个计算域之间相较于常规分析有数据之间的传递，数据的多次传递误差也会加大，因此有必要对流体计算结果验证以保证闸门计算结果的准确性。本文选用经典的重力坝模型用以验证数值模拟水动力结果，数值计算模型以文献[29]中的重力坝-库水地震试验为参照，图5为对应的数值模型，图中水体高度H=95 m，长度L=300 m，流体域上表面为开放边界，不可压缩，坝体与水体接触边界设置为流固耦合边界，流体域其余壁面为固壁边界。

图5 重力坝-库水数值模型

在Workbench设置流固耦合数值分析模块，通过对水体-重力坝双向流固耦合迭代耦合求解，得到地震动水压力沿坝高分布结果。文中设置11个控制点于坝体，在深水处加密。图6为流固耦合数值解与文献[29]结果以及规范[19]结果的对比，流固耦合数值解与实验结果基本吻合，其平均误差不超过3%，说明本文数值模型有较高的精度和准确度，能用于后续相关分析研究。

图6 重力坝动水压力分布图

4 弧形钢闸门地震响应分析

4.1 弧形钢闸门流固耦合模型

闸门-水体系统地震动力响应分析与场地地震分析的本质类似，就是考虑流固耦合效应的外源系统波动效应[30]，因此在选取计算域时需要截取有效计算区域。图7为闸门-水体耦合模型，计算模型按设计工况分析计算，图中H为设计水位高程16.6 m，计算水体长度L选取大于3倍水体的高度，为60 m，设置水体-结构接触面为流固耦合面，水体顶部为开放区域，其余壁面均为固壁边界。为查看闸门不同位置的动参数响应情况，选取9个闸门关键构件予以监测，分别为闸门顶部面板M1、面板中下部M5，上主横梁翼缘M4、下主横梁翼缘M6、上主横梁腹板M7、纵梁腹板M8、下主横梁腹板M9以及上支臂M2和下支臂M3。

图7 闸门-水体耦合模型

4.2 动位移响应

图8显示了最终时刻闸门考虑流固耦合效应的地震动位移响应云图。图8显示，弧形闸门门叶位移主要发生在下半部，位移最大值集中于门叶中下部，受到约束的支铰端位移最小。结构位移最大值为17.32 mm，规范[19]设计允许值为l/600=19685/600=32.81 mm，满足刚度要求。而且，闸门结构动位移响应区域主要集中在闸门顶部以及上下主梁之间。

图9展示的是闸门在地震荷载作用下整体最大位移点位移时历曲线。由图9分析可知在地震动荷载作用下闸门整体历时位移响应在0.16 s时达到峰值27.5 mm，随后闸门结构动位移响应振幅减小，最终达到稳定17.3 mm，与云图位移数值一致，地震动位移响应峰值与稳定值之间振动幅度数值相差37.1%。从而说明闸门动位移响应峰值与云图显示结果有差异。

图8 闸门整体位移云图

图9 闸门整体最大位移点位移时程曲线

图10为闸门结构板、梁、柱构件在地震动荷载作用下的位移时程曲线，纵轴为位移d(m)，横轴为时间t(s)，其中点M5、M8和点M4、M7以及点M6、M9两点之间位移基本一致，故位移曲线绘于同一张图上。图8显示位移云图中最大位移17.3 mm，位于上下主横梁之间，与总位移历时曲线峰值27.5 mm不符，由图10通过对各构件的位移分析可知，闸门在地震动荷载作用的过程中最大动位移响应发生在面板顶部M1处，数值为27.5 mm，发生的时间为0.16 s前后0.02秒的子步时间范围内，与动位移响应峰值相对应；面板顶部M1位移波动幅值为各监测构件最大，在±10 mm数值范围内波动，其次为面板中部M5和纵梁腹板M8，上支臂M2和下支臂M3径向位移波动最小。由此可知位移响应云图并不能反映闸门在动力响应历时过程中的真实力学特征，结构响应极值有可能发生在响应历时的某一时刻；面板顶部M1出现结构动位响应峰值，说明面板顶部区域为动力响应薄弱环节。

图10 闸门各构件位移时程曲线

分析面板顶部M1、面板中部M5以及梁构件M4～M9位移历时曲线可以发现，最大动位移响应发生在闸门面板顶部，随后减小，在上主横梁处达到极小值，越过上主横梁极值点处后动位移数值又开始逐渐增大，随后动位移数值在面板中部出现极值，然后动位移数值开始慢慢减小又在下主横梁处达到极小值，其中下主横梁处极值大于上主横梁处达到的极值，闸门面板顶部的极值大于面板中部出现的极值。闸门面板位移数值沿门高方向分布呈正弦曲线分布。

4.3 动应力响应

图11为闸门最终时刻整体应力云图。由图11可知，总体应力响应最大值为483.72 MPa，位于支臂与下主横梁翼缘与支臂连接区域，有应力集中点，应当予以剔除。对闸门主要位置构件进行强度验算时，取用1.5的构件验算安全系数。弧形闸门应力分布空间效应强，所以对其应力验算时以第4强度理论进行验算，数值计算结果为Mises应力，闸门构件验算结果如表1所示。由表1可知，闸门下主横梁翼缘M6验算空间应力数值最大，即闸门下主横梁为结构强度控制的主要受力构件，但数值均在安全容许设计值范围之内，因此仅就最终计算结果判定该闸门强度满足设计要求。

图11 闸门整体应力分布云图

表1 闸门主要构件强度评判

图12为闸门动力响应板、梁、柱构件应力历时曲线。分析图11、12可知，闸门应力云图显示应力最大值主要集中在上、下主横梁区域，其次为支臂应力数值较大，面板上部应力数值较小；各构件的动应力时程曲线结果显示在地震作用开始的1 s时域内动应力振动幅值较大，伴随着地震作用衰减其动应力响应逐渐趋于平稳；对闸门应力强度校核时所使用的数值为闸门动应力响应终点时刻数值，然而在构件M4、M6动应力时程曲线变化过程中，上主横梁翼缘M4在0.06 s时数值达到222.1 MPa，下主横梁翼缘M6在0.02 s时数值达到258 MPa。从而说明闸门结构动力响应主要集中在初始响应1 s时域内的闸门上、下主横梁区域；在该时域内，构件M6动应力数值已经超过考虑安全系数的材料容许应力，有强度破坏的可能，结合应力分布云图不难看出，上下主横梁翼缘与支臂交接处为动力响应薄弱区域，应当采取有效构造措施进行加固处理，增强局部强度。

图12 闸门各构件Mises应力时程曲线

表2为考虑流固耦合效应的闸门各构件地震动力响应分析的最大计算数值结果与规范[19]计算结果的比较。

表2 闸门各构件位移及应力的规范计算结果与动力数值结果对比

由表2可以看出，考虑流固耦合作用的弧形钢闸门动力响应与规范结果计算差值较大，其中面板顶部M1处动位移与静位移比值为2.19，面板中部区域M5应力比值达到2.44，结构响应两者比值均值为1.52，表中比值数值整体上部区域大于下部，闸门底部区域结构响应增幅较小。由此可见静动力位移差异最大点在门顶面板中点处，静动力应力差异最大点在门叶中部为了进一步分析地震激励下弧形钢闸门地震动水压力，分别以考虑流固耦合和规范[19]方法计算不同水深地震动水压力，结果见表3和图13。其中规范中动水压力采用Westergaard动力法计算，该方法用于计算作用垂直刚性坝面于无限长流体的动水压力，公式如下：

图13 流固耦合方法与规范方法计算的不同水体高度闸门动水压力曲线

表3 流固耦合方法与规范方法计算的不同水体高度闸门动水压力比较

(6)

式中：Ph为水深h所对应的动水压力，kPa；ag为地震加速度，m/s2；ρw为水体密度，kg/m3；H为总水深，m；h计算水深，m。

Westergaard公式多用于垂直刚性面，对于弧门根据闸门规范[19]倾斜迎水面需要考虑折减系数k(k=θ/90，θ为倾斜面与迎水面夹角)，算例中闸门面板和水面之间的夹角θ=75.9°。图13中数值模拟结果与考虑折减系数后的Westergaard解沿弧形闸门高度的分布规律呈抛物线趋势，计算结果显示随水深增大动水压力增加速率加大，且数值模拟数值增加幅度大于公式计算值，结合表3可知在水深12.6 m处两者比值达到最大1.13，数值计算最大地震动水压力占静水压力比值为14.74%，结合文献[31]、[19]，对比分析图6与13可以发现，弧形闸门地震动水压力分布曲线较竖直壁面情况向左移动。不难看出，与规范计算结构响应相比较，流固耦合作用下闸门结构应力、位移响应增加明显，耦合作用对结构响应影响显著；现行规范计算结果难以真实反映地震动水压力在弧门的空间分布特征，更不能展现地震荷载于弧门的历时特性。

由上述分析可知，有必要考虑流固耦合作用对闸门安全的评价，且从图10与图12分析结果可以得知，结构动力响应最大值会出现在响应时域上某一时刻，且数值超过了设计安全容许值，所以有必要通过动力响应时程分析对闸门结构的安全性进行分析评价和研究。

5 结 论

在考虑地震激励下的弧形钢闸门结构设计中，结构动力分析显得尤为重要。本文中对考虑流固耦合响应的弧形钢闸门地震动力响应分析进行研究，得到的主要结论如下：

(1)所建流固耦合模型对水动力的计算求解计算结果与实验结果对比显示其最大误差不超过3%，说明本文的流固耦合模型流体域动水压力计算能够达到较高的精度，可以应用于弧门结构动力响应研究。

(2)通过流固耦合数值模型对该弧形钢闸门各构件进行动力响应监测，得到闸门结构最大动位移发生在门顶面板中心，为27.5 mm，小于现行规范容许值32.81 mm，满足刚度要求；最大动应力发生在下主横梁跨中区域，为258 MPa，不满足取用安全系数强度要求应力值230 MPa；而规范法计算的最大动位移位于竖梁腹板，为18.12 mm，最大应力位于下主横梁翼缘，为198.7 MPa，均满足强度和刚度要求，表明现行规范计算结果偏危险。

(3)考虑流固耦合作用对闸门结构动力响应、弧门地震动水压力分布特征进行分析，结果表明，考虑流固耦合的动力计算结果是规范最大位移计算结果的2.19倍，最大应力的2.44倍，且在下主梁应力响应处动力响应结果已超出规范设计容许值；规范法对于高水头地震动水压力的计算结果偏大，对表孔低水头的计算结果又偏小，难以真实地反映地震动水压力于闸门面板的分布特征。闸门结构耦合时程分析更有利于结构安全评价与研究，因此考虑流固耦合作用对于闸门结构的分析研究有重要意义，能为闸门动力设计问题提供可借鉴的数值模型。