基于双向流固耦合的构架避雷针风振响应分析与减振优化设计

赵桂峰， 曹鹏毅， 石雨昊， 张 猛， 刘 冉

(郑州大学 土木工程学院，郑州 450001)

变电站是电网系统中接受、转换及分配电能的关键场所，保护其安全运行十分重要。变电站服役期间，雷击作为一种常见的自然灾害，严重威胁着站内电气设施的安全。为了防范由雷击产生的破坏事故，通常将避雷针安装在变电构架上部，形成构架避雷针结构[1]，如图1所示。

现代大型变电站中常用的圆钢管构架避雷针细且长，高度常在10～30 m，属于高耸结构[2]。由于其高、细、柔的结构特点，钢管避雷针对风荷载十分敏感，在服役期间会经历大量的风致振动引起的应力循环，甚至产生涡激共振现象，从而引起结构的失效[3]。近年来，国内已发生多起变电站构架避雷针风致坍塌事故[4-5]，严重影响了变电设施的正常运行，因此，针对构架避雷针的风振响应特性进行深入分析进而提出合理的减振控制措施，就显得十分必要。

图1 变电构架及避雷针结构Fig.1 Substation framework and lightning rod structure

针对变电站构架避雷针的风毁事故，已有研究者开展了一些研究：邵波等[6]利用Fluent软件模拟分析了某750 kV的变电构架避雷针结构在风场中的绕流现象，分析表明，避雷针各节段圆管结构均表现出涡脱现象,且圆管直径越小，涡脱现象越明显；孙涛等[7-8]通过对两起构架避雷针的断裂事故的分析，得出横风向涡激共振是引起构架钢管避雷针在较低风速下破坏的主要因素。由上可知，构架避雷针结构的横风向涡激效应在实际设计中不可忽略。

目前，针对工程结构或构件涡激振动理论和数值计算的研究已较多：李永乐等[9]以Fluent为工具，采用CFD(computational fluid dynamics)和CSD(computational structure dynamics)耦合的方法，分析了二维方柱绕流的涡激共振现象；唐友刚等[10]利用Fluent软件求解黏性Navier-Stokes方程、圆柱涡激振动的结构动力响应方程，运用动网格技术，实现流固耦合，模拟分析了圆柱单自由度和两自由度涡激振动；康庄等[11]基于软件OpenFOAM中的自编程PimpleDyFoam求解器，结合SST(shear stress transfer)k-ω湍流模型实现了对弹性支撑的双自由度圆柱涡激振动的数值模拟研究；Dahl等[12]通过试验的方法研究了雷诺数为1.1×104～6×104圆柱体两向自由度涡激振动，发现横向位移响应幅值最大可以达到1.35倍直径；Jauvtis等[13]通过试验的方法研究了不同质量比下的单自由度和双自由度的圆柱涡激振动情况；Singh等[14]采用有限元法对低雷诺数下的弹性支撑圆柱体进行了涡激振动响应数值研究。总结上述研究可知，当前采用单向流固耦合进行涡激振动数值分析的较多，主要原因是：所研究的固体结构在流场中的位移响应较小，因此可以忽略固体变形对流场产生的影响[15]。然而，对于本文拟研究的构架避雷针结构，由于其高、细、柔的结构特点，在风场中结构的变形不可忽略。为更精确地模拟避雷针结构在风场中的响应本文拟采用双向流固耦合[16]方法研究构架避雷针结构在不同风速下的顺风向和横风向风振响应，分析二者对结构的影响特点，进而针对性的提出减振优化设计建议，以期为工程应用提供参考。

1 双向流固耦合求解方法

1.1 数值模拟计算方程

1.1.1 流体计算方程

考虑空气与避雷针之间的流固耦合作用时，流体区域和固体区域应分别满足其基本控制方程，对于较小的风速情况，由于密度变化较小，一般采用控制方程为不可压缩N-S方程，在直角坐标系下，流体域的连续和动量方程表达如下：

连续性方程

(1)

式中，μx，μy，μz分别为流体在X，Y，Z3个方向上的速度分量。

动量方程

(2)

式中：f为流体在X,Y,Z3个方向上所受质量力，m/s2；p为流体压强；ρ为密度；v为运动黏性系数。

1.1.2 固体控制方程

对于固体部分由风速引起避雷针结构的振动响应，其控制方程为

(3)

式中：Ms为固体质量矩阵;Cs为固体阻尼矩阵；Ks为固体刚度矩阵；r为固体位移；τs为固体所受应力。

1.1.3 流固耦合方程

为保证流固耦合遵循最基本的守恒原则，所以在流固耦合交界面处，应满足

(4)

式中：τ为应力；d为位移；q为热流量；T为温度；下标f为流体；下标s为固体。

求解流固耦合问题需要分别列出两者的运动方程并联立。基于流体的连续性方程和动量方程，得到流体的运动方程为

(5)

式中：Pe为流体压力；Ue为结构位移；Mef为流体的质量矩阵；Kef为流体的刚度矩阵；Cef为流体的阻尼矩阵；ρRe为耦合质量矩阵。

当固体受到流体施加的压力后，振动方程为

(6)

式中：Fe为外部荷载矩阵；Re·Pe为流体压力荷载矩阵。通过联立式(5)和式(6)，得到流固耦合的方程为

(7)

式(7)表明，在流固耦合面上，各个节点具有同样的位移和压力自由度，当流体区域和固体区域耦合面上的接触点的解确定后，通过式(7)可以确定耦合面上的解向量，从而解决流固耦合动力学问题。

1.2 耦合算法

双向流固耦合一般分为弱耦合和强耦合，所谓的强耦合就是在建立的同一个方程中求解固体结构和流体域以及耦合作用，在同一时间步长内同时求解所有变量，强耦合计算结果精度较高，但是求解复杂，计算成本较高，在实际应用中较少。一般适用耦合计算的理论分析。弱耦合是将流体域和固体域分开求解，有明显的先后求解顺序，并通过流固耦合截面传递流体域和结构域的数据从而实现耦合计算。相对于强耦合，弱耦合可以充分利用现有通用的流体和结构软件，并且可以分别对每一个软件单独制定合适的求解方法。比较适用于对耦合场的数值计算。本文采用弱耦合算法，其一个时间步具体流程图如图2所示

2 流固耦合数值模型

2.1 工程背景

本论文以某500 kV变电站220 kV构架避雷针为工程背景，该变电构架避雷针是从2010年5月开始正式运行，在2014年12月冬季大风气候下发生了断裂坍塌事故，共服役4年7个月。其结构断裂事故现场如图3所示。为研究风荷载对构架避雷针的影响，本文选取此工程背景下的构架避雷针，进行建模分析。

图2 耦合算法流程图[17]Fig.2 Flow chart of coupling algorithm

图3 构架避雷针结构断裂事故现场Fig.3 Accident site of structure fracture of lightning rod of frame

2.2 数值模型

2.2.1 结构建模

模型的建立包括结构模型和流体模型，根据文献[18]可将构架避雷针的结构模型从上到下分为5部分，考虑到双向流固耦合分析计算量大，数据传递工作量较为频繁庞杂，因此本文模拟分析时采用1∶10的缩尺模型以简化计算工作量并提高计算效率。建模时通过改变材料弹性模量和密度来保证模型和原型结构抗弯刚度和质量缩尺比要求，其具体数值如表1所示。

表1 构架避雷针模型参数

构架避雷针的整体模型如图4所示。由伊运平的研究可知，构架避雷针的支座可近似为矩形截面，经计算截面抵抗矩可得，矩形长边的一面为弱轴面，短边的一面为强轴面。将模型导入ANSYS Workbench工作平台中的瞬态结构分析模块(Transient Structural)。在结构的底部采用固定端的方式施加约束，结构壁面设置为流固耦合面以实现和流体数据的相互传递。最后，在 ANSYS Workbench 中的ANSYS Meshing模块中采用MultiZone进行结构网格划分，网格总体尺寸设为4 mm，划分完后构架避雷针气弹模型的节点及网格数量：节点数为24 879，单元数为24 744。

图4 构架避雷针模型(mm)Fig.4 Model of frame lightning rod(mm)

图5和图6分别为相贯节点部位和各避雷针杆段连接部位网格示意图，从图5和图6中可以看出结构网格均为四边形网格，查看模型的单元质量因子，数值为0.98，说明有较好的网格质量。

图5 相贯节点部位网格示意图Fig.5 Grid diagram of intersecting node

图6 杆段连接部位网格示意图Fig.6 Grid diagram of connecting part of rod section

2.2.2 流场建模

流场二维模型如图7所示，流场域有5个边界条件，进出口边界，上下边界和避雷针外表面边界，避雷针形心到进口和上下表面的边界距离都为10D，到出口的距离为25D。用Hypermesh软件对模型进行3D网格划分，如图8所示。然后导入Workbench中的Fluent模块进行基本设置[19]：

计算方法——瞬态计算；

计算模型——Realizable k-ε湍流模型；

边界条件——流域进口为速度进口，速度方向垂直于进口边界表面，速度出口采用自然出流。对于进口边界的风速和湍流特性编制程序并采用UDF输入。

动网格——打开动网格选项，动网格是采用弹簧光顺法进行更新。

图7 流场二维计算域Fig.7 Two dimensional computing domain

图8 流体域网格模型Fig.8 Fluid domain mesh model

2.2.3 流固耦合求解设置

在Workbench平台中的System Coupling中进行流固耦合求解参数设置，利用Transient Structural中的流固耦合面和Fluent中的结构壁面进行耦合数据传递，设置求解时间为20 s，时间步长为0.01 s。最后设置求解顺序，先计算流体域再计算固体结构。

2.3 网格无关性的验证

为保证数值模拟的正确性，首先对原型构架避雷针气弹模型进行双向流固耦合计算，将计算结果与伊运平的试验数据进行对比来验证该计算方法的正确性。考虑到构架避雷针底部支座为矩形支座，因此需要考虑迎风面为强轴和弱轴两种情况分析其在不同风速下的风振响应。定义0°风向角为迎风面是强轴的情况，则90°风向角为迎风面是弱轴的情形。具体情形如图9所示，沿构架避雷针的强轴(S轴)和弱轴(W轴)两个方向布置测点。其具体情形如图10所示。模拟0°风向角下平均风风速为1 m/s，2 m/s，3 m/s，4 m/s，5 m/s，6 m/s，7 m/s，8 m/s，9 m/s，10 m/s共10种风速，根据缩尺模型的相似比算出对应的构架避雷针实际风速分别为3.16 m/s，6.35 m/s，9.49 m/s，12.65 m/s，15.81 m/s，18.97 m/s，22.14 m/s，25.30 m/s，28.46 m/s，31.62 m/s。下文表述的风速数值均按实际风速给出。

图9 构架避雷针风向角图示Fig.9 Wind direction angle of frame lightning rod

采用当地的基本风速23.83 m/s的工况进行网格无关性的检查，边界层厚度设置为1.8×10-5m。其网格的划分，如图11所示。

本文采用粗网格、中网格和细网格三种网格来进行网格无关性的验证，相关网格参数如表2所示。

表2 网格参数

图10 构架避雷针应变测点布置(mm)Fig.10 Layout of strain measuring points for frame lightning rod (mm)

图11 构架避雷针计算域网格Fig.11 Computational domain grid of frame lightning rod

图12为不同尺寸网格下构架避雷针各个监测点的应变变化，并与试验值进行对比。

从图12中可以看出，不同尺寸的网格下构架避雷针各个监测点的应变变化基本保持不变，通过和伊运平研究中的试验值对比分析可知：粗网格的最大偏差值为20%；中网格的最大偏差为8%；细网格的最大偏差为7%。考虑到时间的成本,最后决定采用中等的网格尺寸进行计算。

图12 不同尺寸网格下构架避雷针的应变Fig.12 Strain of lightning rod in frame with different size grid

3 构架避雷针的风振分析

3.1 结构基本风速下的受力特性分析

对构架避雷针气弹模型进行模态分析，在不施加外力和阻尼的情况下，仅对结构底部进行固定，计算结构的模态频率如表3所示。

表3 构架避雷针前4阶频率

结合伊运平研究计算得到的当地10 m高度处，重现期为 50年的基本风速为23.83 m/s，相对应于模型风速为7.5 m/s。通过双向流固耦合模拟计算原型构架避雷针在0°风向角下基本风速中各测点的应变时程，并对其进行傅里叶变换，可得到应变功率谱。图13给出了原型避雷针的顺风向测点的应变时程和应变功率谱。

图13 23.83 m/s风速下实际构架避雷针顺风向各监测点应变响应Fig.13 Strain response of lightning rod along wind direction of actual frame lightning rod under wind speed of 23.83 m/s

由图13中各测点应变功率谱对比可知，构架避雷针上部两段振动形式最为复杂，其振动响应包含前4阶振型成分，对应的自振频率分别为2.51 Hz，5.53 Hz，13.73 Hz，20.13 Hz，第Ⅲ段包含第1阶、第2阶、第4阶振型成分，下部两段主要以第一阶振型为主。所以避雷针结构的高频响应主要体现在结构的上部杆件，其振动由多阶振型共同控制，振动形式复杂，对结构响应不利。

对基本风速下构架避雷针的横风向风振响应进行监测分析，结果如图14所示。

图14 23.83 m/s风速下实际构架避雷针横风向各监测点应变响应Fig.14 Strain response of actual frame lightning rod across wind direction at 23.83 m/s wind speed

从图14中可以看出，构架避雷针的横风向响应主要表现为由卡门涡街引起的横风向涡激振动，结构的应变响应为零应变附近的震荡，结构响应频率基本为各阶基频。其中在基本风速下构架避雷针的首段振动幅度最为剧烈， 为受力最大位置。构架避雷针的第Ⅱ段振动形式最为复杂，其振动频率包含结构前4阶频率。

3.2 构架避雷针在0°和90°下的风振响应分析

图15给出0°风向角原型构架避雷针各测点在不同风速下的应变值均方根响应及绝对值最大值。图16给出顺风向与横风向应变响应的对比。

在0°风向角下，顺风向测点位于W轴上，横风向测点位于S轴上。从图15中可见，顺风向响应中，各测点的应变值随风速增大而增大，在各级风速下，W-3测点的应变值最大，其应变绝对值最大值可达到1.15倍的W-5测点的应变绝对值最大值。说明由于避雷针结构在Ⅲ段和Ⅳ段交界处界面突变较大，该部位结构的受力非常不利，在设计与运行阶段要重视该位置的校检与检验；由图16可知，该避雷针结构易发生横风向涡激振动的杆段是Ⅰ段、Ⅱ段和Ⅳ段，首次出现涡激振动的风速均为12.64 m/s，结构的横风向响应和顺风向响应处在同一数量级，对于Ⅰ段构架避雷针结构由横风向涡激振动引起的应变最大值甚至要大于顺风向的响应，尤其在锁定区间，当风速达到18.97 m/s的时候，S-1点应变绝对值最大值是W-1点应变绝对值最大值的5.12倍，因此构架避雷针的设计中横风向的风振响应同样要重视。从图中还可以看出，对Ⅰ段、Ⅱ段和Ⅳ段避雷针结构，涡激振动锁定区间均出现在12.64～25.30 m/s。当风速为18.97 m/s时，Ⅱ段和Ⅳ段构架避雷针的涡激振动幅值最为显著；对于Ⅰ段避雷针结构，在风速为22.13 m/s时，涡激幅值最为显著。对比各段避雷针发生涡激共振时的应变绝对值最大值，可知Ⅰ段和Ⅱ段避雷针涡激共振强度最大，受力最为不利。

图15 避雷针结构0°风向角各段底部监测点各风速下应变值Fig.15 Strain value of lightning rod structure at the bottom of each section at 0° wind direction angle under each wind speed

图16 避雷针结构0°风向角下顺风向与横风向响应对比Fig.16 Comparison of responses between downwind and crosswind of lightning rod structure under 0° wind direction

图17给出90°风向角原型构架避雷针各测点在不同风速下的应变值均方根响应及绝对值最大值。图18给出顺风向与横风向应变响应的对比。

在90°风向角下，构架避雷针的顺风向测点位于S轴上，横风向测点则位于W轴上。从图中可见90°风向角下构架避雷针的受力变化情况和0°相同。

综合0°风向角和90°风向角两种情况来看，当来流风速达到12.64 m/s的时候，构架避雷针Ⅰ段、Ⅱ段和Ⅳ段开始发生涡激共振现象。当风速达到18.97 m/s的时候，Ⅰ段、Ⅱ段和Ⅳ段构架避雷针涡激振动幅值最为显著。对于Ⅰ段和Ⅱ段构架避雷针，在风速为12.64～25.30 m/s的涡激振动锁定区间。其横风向的振动幅值要明显大于顺风向的幅值，为结构受力最不利位置。

图17 避雷针结构90°风向角各段底部监测点各风速下应变值Fig.17 Strain value of lightning rod structure at the bottom of each section at 90 ° wind direction angle under each wind speed

图18 避雷针结构90°风向角下顺风向与横风向响应对比Fig.18 Comparison of responses between downwind and crosswind of lightning rod structure under 90° wind direction

4 基于双向流固耦合分析的构架避雷针横风向抗风优化设计

由上述分析可知，该构架避雷针顶部的Ⅰ段、Ⅱ段最易发生横风向涡激振动，由此引起的振动响应对结构安全运行非常不利，有必要对横风向的涡激振动进行控制控制涡激共振的方法主要分为主动控制和被动控制两大类，其中主动控制是指在结构物投入使用期间人为的或者通过某种动力或操控装置扰乱流场；被动抑制措施是指在前期设计中改变结构物截面形状，或添加附属可绕流装置。考虑到构架避雷针结构服役环境，其风场难以人工操控，因此只能使用被动控制措施。而被动控制措施中外附螺旋导板由于抑振效率高，施工简单的原因得到广泛认可。实际工程中外附螺旋导板的结构参数(螺距、覆盖率等)均与结构的外径密切相关，不同尺寸的结构需要设计不同的导板参数，目前还未见有关于避雷针结构的外附螺旋导板以控制横风向涡激振动的研究，本文拟对此进行相关研究，以为该类结构的振动控制提供一定的理论参考。

4.1 建 模

本文采用螺头数为3，螺板高度为0.25D[20](其中，D为避雷针圆管直径)的外附螺旋导板，模型如图19所示，深色部分为添加的螺旋导板。左上角的为螺旋导板的放大图。定义螺距为两条侧板之间立管的轴向长度。通过改变螺距和导板覆盖率的方式来设置不同的工况进行对比分析。以各个工况下构架避雷针的应变幅值作为优化效果的衡量标准来寻找此工程背景下的最优螺旋导板布置参数。

图19 优化后避雷针构件Fig.19 Optimized lightning rod components

4.2 网格无关性验证

由于导板的尺寸比较小，网格的划分方法对数值模拟的精度影响比较大，因此在这里对添加了螺旋导板的构架避雷针进行网格无关性的验证。以基本风速下螺距为8D，盗版覆盖率为50%的工况为例。中细网格的尺寸见表2。验证结果图20所示。

图20 不同尺寸网格下外附螺旋导板构架避雷针应变变化Fig.20 Strain analysis of lightning rod with spiral guide plate under different grid sizes

从图20中可以看出，不同尺寸网格下构架避雷针的应变变化幅度极小，差异最大值为5%，其中大部分的差异值在1%左右。综合考虑到计算成本，本文采用中网格进行模拟计算。

4.3 不同螺距下的导板控制效率分析

控制螺旋导板覆盖率为50%，螺距为从大到小依10D，8D和5D(其中，D为螺旋导板所在位置构架避雷针的截面直径)。提取沿横风向5个测点的应变并将其与无外附螺旋导板的值进行对比。结果如图21和图22所示。其中，10D-S-1表示螺距为10D的构架避雷针中S-1测点。

图21 0°风向角各横风向监测点应变响应Fig.21 Strain response of each cross wind fixed point at 0 ° wind direction

图22 90°风向角各横风向监测点应变响应Fig.22 Strain response of each cross wind fixed point at 90° wind direction

从图20中结构的横风向响应来看，0°风向角下，构架避雷针的第Ⅰ、第Ⅱ和第Ⅳ段部分和90°风向角下构架避雷针的第Ⅰ、第Ⅱ和第Ⅴ段部分中外附螺距从10D减小到8D抑制效果增大，从8D～5D抑制效果减小。尤其对于0°风向角下的构架避雷针的第Ⅰ、第Ⅳ段部分和90°风向角下构架避雷针的第Ⅰ、第Ⅱ、第Ⅳ段部分。外附螺距为8D的优化后结构不仅没有出现涡激共振现象，而且最大的横风向响应也减小了60%以上。结合0°和90°风向角下构架避雷针的风振结果来看，螺距为8D的外附导板优化措施可以有效的减小构架避雷针的横风向涡激共振。尤其对于构架避雷针上部两段部分，优化效果显著。

4.4 不同导板覆盖率下的导板控制效率分析

设置螺旋导板覆盖率为0(无导板)，30%和50%三种情况进行对比分析，结果如图23和图24所示。

图23 0°风向角各横风向监测点应变响应Fig.23 Strain response of each cross wind fixed point at 0 ° wind direction

图24 90°风向角各横风向监测点应变响应Fig.24 Strain response of each cross wind fixed point at 90 ° wind direction

从图23和图24中可以看出，在0°风向角下，构架避雷针的第Ⅰ、第Ⅲ段部分和90°风向角下构架避雷针的第Ⅰ、第Ⅱ、第Ⅳ段部分中螺旋导板覆盖率从0增加到30%时抑制效果在增加，覆盖率从30%增加到50%时抑制效果在减小。外附覆盖率为30%的螺旋导板的优化构架避雷针在0°风向角下的第Ⅰ、第Ⅳ段部分和90°风向角下的第Ⅰ，第Ⅱ、第Ⅳ段部分最大的横风向响应减小了60%以上。从图中可以看出，螺旋导板覆盖率对构架避雷针涡激响应的影响效果没有螺距的影响效果显著。

4.5 优化前后结构顶端位移的变化分析

构架避雷针作为一种细长构件，结构的顶端位移变化情况也是一种极具说服力的数据。提取添加螺旋导板前后构架避雷针的顶端位移变化情况并绘制成图25～图26所示。

从图25和图26中可以看出，无论是0°还是90°，结构的顶部位移变化情况和应变变化情况基本保持一致，在导板的覆盖率为50%不变的情况下，当螺距从5D增大到8D时，构架避雷针的横风向响应逐渐减小，当螺距从8D增大到10D时，构架避雷针的横风向响应逐渐增大；当保持外附螺旋导板的螺距为8D不变时，螺旋导板的覆盖率从0增加到30%的期间，构架避雷针的横风向顶端位移逐渐减小，当导板的覆盖率从30%增加到50%的期间，构架避雷针的横风向顶端位移逐渐增大。综上所述，当外附螺距为8D，覆盖率为30%的螺旋导板时，构架避雷针的横风涡激振动得到明显的抑制。

图25 不同螺距下构架避雷针顶部位移Fig.25 Top displacement of frame lightning rod under different pitch

图26 不同覆盖率下构架避雷针顶部位移Fig.26 Top displacement of frame lightning rod under different coverage

5 结 论

本文以某500 kV高压变电站典型圆钢管构架避雷针结构为例，采用双向流固耦合方法模拟分析了避雷针的顺风向及横风向风振响应，并采用外附螺旋导板的方法对其横风向风振响应进行了减震优化分析，主要结论如下：

(1)采用双向流固耦合方法可以有效模拟圆钢管构架避雷针的风振响应特点。分析表明，构架避雷针的顺风向和横风向风振响应较为突出，且受多阶振型影响较大，尤其是避雷针顶部杆件的高频振动响应显著，在某些条件下甚至占据主导地位。

(2)无论是0°风向角还是90°风向角，结构涡激共振的锁定区间都在12.64～25.30 m/s。其中当风速为18.97 m/s时，结构的涡激振动幅值达到最大，该峰值甚至超过设计基本风速(23.83 m/s)作用下结构的顺风向风振响应值。

(3)在易发生涡激振动的杆段设置螺旋导板可有效抑制杆段横风向涡激振动的产生，从而减小对结构不利的风振响应；但螺旋导板的结构参数与覆盖率需要通过优化设计以确定其合理值。针对本文分析的避雷针结构，经分析确定当外附螺旋导板螺距为8D，覆盖率为30%、板高0.25D时对结构的横风向涡激振动控制效果最好，一方面抑制了涡激振动的发生；另一方面使结构的最大横风向响应减少了60%以上。