基于流固耦合的南水北调渠道结构动力破坏机理研究

梁钟元，谢晨龙，许新勇，3

(1.水利部南水北调规划设计管理局，北京 100038；2.华北水利水电大学水利学院，河南 郑州 450046；3.河南省水工结构安全工程技术研究中心，河南 郑州 450046)

0 引 言

南水北调中线工程是我国水资源优化配置的重大工程，全长1 432 km，是典型的长距离线性输水工程，跨越河南、河北、天津、北京4省市，其中部分渠道穿越地震高烈度区，一旦发生破坏，水体漫溢，将对工程沿线造成无法估量的损失。以往渠道工程地震计算时较少考虑水体的影响，但中线工程渠道水体流量大、水位深，其影响不容忽视[1]。

目前，国内外学者对流固耦合和结构抗震方面均做了大量研究，但对综合考虑流固耦合与地震动分析的研究相对较少。王铭明等[2]对重力坝坝体-库水耦合系统进行了地震计算，将结果与Westergaard公式解析解比较发现，流固耦合计算有上部偏大，下部偏小的趋势；刘云贺等[3]应用振型叠加反应谱法对南水北调中线渡槽工程进行了抗震分析，结果表明渡槽的动位移、动应力和静动综合反应满足设计要求；刘晶波等[4]基于波动方程推导了粘弹性边界方程，通过算例证明了三维粘弹性人工边界的良好精度和稳定性。以上文献均主要集中在结构抗震的某一方面，对综合考虑流固耦合效应和合理的地基边界条件等影响下结构地震响应方面的研究较少，尤其对渠道工程相关研究更少。

本文以南水北调中线潮河六标渠道工程为例，建立水体-结构-地基耦联体系的渠道模型，采用Fluid 30流体单元模拟流固耦合作用，引入粘弹性动力人工边界，进行渠道衬砌结构地震动力破坏机理的研究，拟为南水北调中线渠道工程的安全运行维护和提升抗震措施提供有益的理论支撑和借鉴。

1 流固耦合方法和Fluid 30单元

流固耦合问题研究要同时考虑动力学方程、Navier-Stokes方程及流体连续性方程。假设流体具有可压缩和非粘滞性，流体密度和声压一致。通过对以上方程简化，可得水体单元波动方程为

(1)

式中，c为流体中的声速，在水中取c=1 640 m/s；p为声压，流场中声压随位置和时间发生变化，在声学中，可以通过声压来求速度等物理量；2为Laplace算符，

为解决边界阻尼造成的能量损失以及更完整地描述流固声耦合问题，需要引入能量损失项和流体作用于结构表面的声压荷载项。流固耦合方程可表示为

(2)

式中，vol为积分域；δp=δp(x，y，z，t)为虚的声压改变量(声压的变分)；s为流体与结构的接触面；{n}为垂直于边界面的法向量。

为体现渠道的流固耦合作用，本文选用ANSYS软件中的声学流体单元Fluid30模拟水体[5]，该单元可通过波动方程和界面流固耦合方程实现水体与结构的相互作用，其单元结构形式如图1所示。

图1 Fluid30单元示意

2 工程概况和数值模型

选取南水北调中线潮河六标高填方段渠道工程为实例进行研究，渠道结构的断面尺寸及构成材料如图2所示，建立渠道的水体-衬砌的流固耦合数值仿真模型如图3所示。数值仿真模型的整体结构包括渠道衬砌板、渠堤土、地基及水体，其中衬砌结构自上而下由10 cm厚混凝土面板、复合土工膜、2 cm厚保温板、20 cm厚反滤料层组成。数值模型的渠道衬砌结构采用Solid 65单元离散，渠堤和地基采用Solid185单元离散，渠内水体采用Fluid30单元离散，整体模型单元56 896个，节点62 730个，其中水体单元17 632个，节点19 890个。渠内水位为加大水位7.645 m。渠道结构各项材料参数见表1。

图2 渠道结构断面示意(单位：m)

图3 渠道-水体的流固耦合仿真模型

表1 渠道结构各项材料参数

动力计算时，为了对远域地基的阻尼效应进行更好的模拟，本文模型动力边界采用粘弹性地基边界[4，6-7]，通过MATLAB语言编程输入三维弹簧单元COMBIN14构建，其实现方法如图4所示。

图4 粘弹性边界实现流程及示意

图5 地震加速度时程曲线

3 地震波的选取与输入

该工程的地基场地类别为Ⅰ类，场地基本地震烈度为Ⅶ度，工程场地基岩相应的罕遇地震水平峰值加速度为0.1g，竖向地震峰值加速度取水平向的2/3。根据频谱特性相近原则，选取人工波进行计算，地震波特征周期为0.4 s，地震峰值加速度为0.1g，波长时间为30 s，地震动输入间隔为0.02 s。取沿渠道横向为X向，沿渠道竖向为Y向，沿渠道纵向为Z向。地震波三向加速度的归一化时程曲线如图5所示。

根据粘弹性边界的地基放大效应，同时考虑边界节点上外行波的吸收与入射波的输入，将加速度削减至原来的1/2后，采用Seismosignal地震波处理软件将地震加速度时程进行积分[4]，获得速度和位移的时程，实现等效应力的转化，完成地震波的输入。

4 模态计算

为研究渠道地震动时的结构振动特性和水体的流固耦合效应，提供渠道工程地震动分析的基础，采用Fluid30单元模拟水体，进行了渠道与水体的模态分析。将Fluid30单元方法与传统的Westergaad公式方法的计算结果进行了对比分析，结果见表2。从表2两种不同方法计算的自振频率和振型描述来看，第1阶模型自振频率的Fluid30单元方法相对于Westergaad公式方法，增幅达到了23.73%，其后高阶频率的计算结果，两者相差逐渐缩小，但总体上来说Fluid30单元方法比Westergaad公式方法的频率略大，周期较小。取模态分析的前10阶计算结果，两者的均值相差为10.93%。从Fluid30单元方法与Westergaad公式方法的模型振型分布描述来看，两者具有一致性、对称性。Fluid30单元方法在进行动力计算时，考虑了渠道水体的可压缩性，能更好地反映渠道-水体之间的相互作用，和传统Westergaad公式方法的结果对比也说明其考虑水体的效果合理。因此在动力计算时，可采用Fluid30单元方法来模拟渠道-水体的流固耦合效应，该方法可为水工流固耦合问题提供一种新的求解思路。

表2 自振频率和振型描述

5 渠道地震动破坏机理

为研究渠道结构地震动破坏的机理，考虑渠道-水体的流固耦合效应，同时引入粘弹性动力边界，进行地震动时程计算分析。地震作用与持时密切相关，在地震波动过程中，加速度的变化将造成结构受力状态不时改变，这是时程计算的显著特征。为探究渠道衬砌结构不同位置对地震动力响应的差别，选取渠底中心A、渠坡中心B、渠坡衬砌板顶端C为特征点，进行衬砌地震动态响应分析；同时，为研究地震过程中衬砌板上下表面的动力响应差别，在A、B、C点对应位置的衬砌下表面选取D、E、F3个特征点一并分析，特征点位置如图3所示。

渠道衬砌板上表面各特征点地震动时程相对位移变化规律如图6所示，由图6可知：①地震时程达到10s左右时，各个特征点的位移响应均达到最大值，这与纵向地震峰值加速度出现的时间基本一致；②对比各个特征点位移的变化规律，发现衬砌板上表面的位移响应基本一致，可认为衬砌结构在地震过程中，基本做整体移动。

表3给出了衬砌板上表面各个特征点的最大位移响应与自由场位移，结构上的最大位移响应基本是自由场输入的2倍，验证了粘弹性边界地基对地震波的放大效应[7]；同时，地基辐射阻尼对结构地震响应影响明显。因此，在实际工程的抗震分析时对其影响应予以适当考虑。

表3 衬砌板上表面各特征点的最大位移 cm

渠道衬砌板结构各特征点的应力时程响应如图7所示。从地震持时过程中衬砌板特征点的应力时程变化可知，上、下表面各点的第一主应力变化趋势基本相同。衬砌板上表面拉应力值普遍大于下表面，但在接近分缝的B、E两点，规律相反，下表面所受拉应力会显著大于上表面。

渠道衬砌板和底部的地震动应力最大值如图8所示，地震持时达到10.5 s时，渠道结构的第一主应力达到最大值，位置出现在底板齿墙与渠坡衬砌板下表面交接处，最大值为5.73 MPa。由图8可知，衬砌板和底板上表面应力分布呈现中间较大、两侧横缝分缝处较小的趋势，下表面则为每条纵缝周围局部应力较大，中间逐渐减小。由衬砌板的应力分布整体趋势可知，渠坡衬砌板与底板相交的齿墙和每条纵缝附近是应力值较大区域，动力作用下存在破坏的可能性。

图8 地震动持时10.5 s时衬砌板和底板上、下表面应力最大值(单位：Pa)

6 结 论

本文针对南水北调中线渠道工程地震动破坏问题，考虑水体流固耦合效应和粘弹性动力边界，进行了渠道结构地震动损伤破坏研究，可得到如下结论：

(1)Fluid 30流体单元所计算的结构频率值较传统Westergaard方法的计算值平均增加10.93%，由此可见，动力计算时，水体的流固耦合效应不容忽视。

(2)渠道衬砌板之间的地震动位移相对较小，可视为衬砌结构整体移动；衬砌板上表面拉应力值普遍大于下表面，但在衬砌板分缝处，规律相反，下表面所受拉应力会显著大于上表面。

(3)渠道衬砌板的地震动应力最大值出现在齿墙与渠坡衬砌板交接处，该处为渠道底板横截面、齿墙及衬砌板相交处，截面型式和连接角度突变，为整体结构的薄弱区域，日常运行工作中，应加强该部位混凝土和止水材料的观测和维护。