高耸烟囱风致响应精细化计算方法

孙梦然，刘仰昭，戴靠山，丁志斌，尹业先; SUN Mengran; LIU Yangzhao; DAI Kaoshan; DING Zhibin; YIN Yexian

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

高耸烟囱风致响应精细化计算方法 PDF

- ORCID：
孙梦然 ¹
✉
- ORCID：
刘仰昭 ¹
✉
- ORCID：
戴靠山 ^1,2,3
- ORCID：
丁志斌 ¹
- ORCID：
尹业先 ⁴

1. 四川大学，建筑与环境学院，成都 610065； 2. 四川大学，深地科学与工程教育部重点实验室，成都 610065； 3. 四川大学，破坏力学与防灾减灾四川省重点实验室，成都 610065； 4. 山东电力建设第三工程公司，山东青岛 266100

中图分类号： TU311.3

最近更新：2024-03-07

DOI：10.11835/j.issn.2096-6717.2022.003

摘要

高耸烟囱的风致响应可分为顺风向响应和横风向响应，其中顺风向响应以大气脉动风引起的抖振响应为主，横风向响应以Karman旋涡脱落引起的涡激振动为主，准确地预测和评估这两种风致响应对其抗风设计和结构安全性至关重要。在Tamura提出的二维平面尾流振子模型的基础上进一步推导，将该模型成功运用在三维结构上，提出可用于实际工程结构的有限元迭代计算方法，为高耸烟囱类结构横风向涡振响应的计算提供了新的方法。此外，基于结构的固有模态坐标，建立了适用于高耸烟囱耦合抖振响应分析的有限元CQC频域计算方法，并将频域计算结果与时域计算结果对比。结果表明：有限元迭代计算方法可以有效地计算三维烟囱的涡振响应，烟囱抖振响应频域计算和时域计算结果吻合良好。

关键词

高耸烟囱; 结构风工程; 涡激共振; 抖振; 有限元方法

高耸烟囱被广泛地运用于石油、化工、电力等工业领域。结构柔、质量轻、阻尼小是超高、大长细比的烟囱所具有的特点，这些都使得高耸烟囱对风荷载极其敏感^[

1]。对于烟囱这种高耸塔筒型结构，最为重要的风致响应是大气紊流作用下的抖振响应和结构尾流中旋涡交替脱落引起的涡激振动^{[参考文献 2

百度学术}2]。大气中由脉动风引起的抖振响应在结构抖振现象中占主要地位^{[参考文献 3-4}3-4]。涡激振动是由钝体结构尾流中旋涡的交替脱落引起的^{[参考文献 2

百度学术}2]。在一定风速范围内，结构尾流中旋涡脱落的频率将不再随风速变化，而是接近结构横风向的振动频率，形成所谓的风速锁定区间^{[参考文献 5-9}5-9]，这将引起结构较大振幅的运动，从而发生涡激共振^{[参考文献 2

百度学术}2]。成功预测这两种风致响应是烟囱这类高耸结构抗风设计的关键之一。

Vickery等^[

10-11]在均匀流场中对圆形断面振动柱体的横风向气动力的测量表明，一般存在两种作用力，第1种产生于旋涡自身的脱落，无论圆柱是否振动都存在，带有强迫性质；第2种与圆柱的运动有关，带有自激性质。在振幅很小时，可以在线性随机振动理论的框架内处理圆柱对旋涡脱落力的响应，即对有强迫性质的气动力的响应，但在更大的振幅时，运动诱发的力对气动作用力的贡献变得更为重要，必须在任何现实的模型中加以考虑。对于这种运动诱发的力，许多数学模型已经被提出，其中最经典的是由Hartlen等^{[参考文献 12

百度学术}12]提出的升力振子模型和Tamura等^{[参考文献 13-14}13-14]提出的尾流振子模型。虽然Hartlen等^{[参考文献 12

百度学术}12]提出的升力振子模型对模拟结构的涡振响应具有较高的潜力，但该模型是一个基于瑞利（Rayleigh）微分方程构建出的经验性模型，其各个经验参数都不具备实际的物理意义，而Tamura等^{[参考文献 13-14}13-14]提出的尾流振子模型很好地解决了这一缺陷。Tamura等^{[参考文献 13-14}13-14]针对圆形断面连续系统的涡激振动提出的尾流振子模型包含一个非线性尾流振荡器。该模型将尾流看作一个具有质量和刚度的振子，且进一步考虑了尾流振子的上下摆动(即角位移α)，以及尾流有效长度的周期性变化，可以较好地模拟二维流场中圆柱的涡振响应，是一种为数不多的各参数均具有明确物理意义的非定常数学模型。但令人遗憾的是，尾流振子模型（Tamura）并未被成功地运用在实际工程中，其主要原因为该模型是基于二维平面提出的，直接计算三维结构将面临复杂的积分运算，难以编程。Davenport在准定常的假定下推导了自然风的紊流在细长结构单位展长上产生的抖振力公式^{[参考文献 3

百度学术}3]。结构的抖振响应计算可分为频域法和时域法两大类。传统的抖振频域计算是基于单模态叠加的SRSS法，该方法忽略了模态间的气动耦合，用于计算高耸烟囱这类高柔结构的误差较大。而有限元CQC方法可以考虑多模态间气动耦合效应，与SRSS方法相比，是一种更适用于高耸烟囱的精确方法。脉动风时程对抖振响应的时域计算影响重大，CDRFG方法是由Aboshosha等^{[参考文献 15

百度学术}15]提出的一种连续离散流随机生成技术，该方法改进了不同频率内的湍流流速相关性，得到的湍流谱与目标Von Karman谱具有很高的相似度。李永乐等^{[参考文献 16

百度学术}16]指出时域结果和频域结果的一致性是计算结果可靠性的一种验证。

笔者在Tamura等^[

13-14]提出的二维平面尾流振子模型的基础上进一步推导，考虑了风速和烟囱断面直径的变化，将该模型成功运用在三维结构上，提出了可用于实际工程结构的有限元迭代计算方法，为高耸烟囱这类结构横风向涡振响应的计算提供了一种新的视野。基于结构的固有模态坐标，建立了适用于高耸烟囱耦合抖振响应分析的有限元CQC频域计算方法，并和时域法对比计算了烟囱的抖振响应，两种方法的结果表现出良好的一致性。

1 工程背景

以某大型烟囱为背景。烟囱模型如图1所示，高216.5 m，底部外径24 m，至110 m标高处变为18.2 m，直到标高180 m处，增大到18.3 m到顶。底部壁厚65 cm，顶部壁厚35 cm，最薄处壁厚30 cm。底部有两个开口，一面尺寸为8 m×9 m,另一面尺寸为3.5 m×4 m。标高27.4 m处有两个南北对称分布的入口烟道，开口7.6 m×13 m。底部北面开口和入口烟道开口周围大约2 m区域，壁厚增加到1 m。

图1 烟囱基本结构示意图

Fig. 1 Schematic diagram of basic structure of chimney

烟囱外筒采用ANSYS有限元软件建立其集中质量模型。在ANSYS有限元软件中，烟囱从0到210 m标高按每5 m一段划分单元，210 m到216.5 m划分为一个单元，因此整个外筒划分为43个单元，共44个节点。模型采用的单元是beam4梁单元，通过定义每一个单元不同的断面面积和惯性矩来体现出烟囱外筒断面和壁厚随烟囱高度的变化。烟囱的内筒在标高205 m处与外筒固结，因此将内筒简化为一个附加的集中质量添加在205 m处的节点上。烟囱ANSYS模型的模态分析结果如图2所示。经计算，该模型的X方向一阶自振频率为0.330 23 Hz，与同结构ABQUAS壳单元模型(图3)的0.315 67 Hz差距约为4.5%；Y方向一阶自振频率为0.339 71 Hz，与壳单元模型的0.337 65 Hz差距约为0.61%，均小于一般工程要求的5%误差,可见该梁单元模型能较好模拟实际结构。

图2 ANSYS梁单元模型及模态信息

Fig. 2 ANSYS beam element model and modal information

图3 ABQUAS壳单元模型及模态信息

Fig. 3 ABQUAS shell element model and modal information

2 涡激共振精细化预测计算

2.1 二维典型断面涡激荷载

Tamura尾流振子模型可以较好地模拟二维流场中圆柱的涡振响应，且模型中各参数均具有明确的物理意义。在考虑了风速和烟囱断面随烟囱高度的变化后，烟囱整体是一个处于三维流场中的连续模型，但对于烟囱的每一个节点对应的断面，将其假设是处于二维流场之中，这样每一断面的涡激荷载可以基于经典的Tamura尾流振子模型来考虑^[

13]。烟囱断面Tamura尾流振子模型如图4所示。

图4 烟囱断面Tamura尾流振子模型示意图

Fig. 4 Schematic diagram of Tamura wake oscillator model for chimney section

图4中，C.R.为振子的扭转中心；C.G.为振子的重心；α为尾流振子的角位移；h为尾流振子的宽度。任一断面处的尾流振子模型公式为

\{\begin{array}{l} \ddot{Y} + 2 η \dot{Y} + Y_{} = - m^{*} V_{ω N}^{2} (f α + f \frac{\dot{Y}}{V_{ω N}} + C_{D} \frac{\dot{Y}}{V_{ω N}}) \\ \overset{}{\ddot{α}} - 2 ζ Ω_{K N} (1 - \frac{4 f^{2}}{C_{Y O}^{2}} α^{2}) \overset{}{\dot{α}} + Ω_{K N}^{2} α = - χ \ddot{Y} - \frac{Ω_{K N}^{2}}{V_{ω N}} \dot{Y_{}} \end{array}

（1）

式中：η为结构的阻尼比；Y表示的是任一断面处的无量纲位移，它是断面实际位移y与断面直径D的比值，即Y=y/D。Y是一种简谐位移，可以表示成 $Y = \hat{Y} • s i n (ω_{N} t)$ ，其中 $\hat{Y}$ 表示的是任一断面的无量纲位移振幅；涡激荷载主要激发烟囱的一阶振型^[

11]，ω_N可以取烟囱第一振型的固有圆频率。V_ωN表示任一断面的无量纲风速，即

V_{ω N} = U / (ω_{N} D)

，U表示烟囱任一断面处的来流风速。Ω_KN表示任一断面的无量纲风速V_ωN与涡激共振起振风速V_ωK的比值，即

Ω_{K N} = V_{ω N} / V_{ω K}

。涡激共振起振风速V_ωK是一个固定值，它仅与圆形断面的斯托罗哈数S_t有关，即

V_{ω K} = 1 / (2 π S_{t})

，其中S_t取0.2。m^*代表质量比，即

m^{*} = ρ D^{2} / (2 m)

，ρ代表空气密度，m代表任一断面的单位长度质量。经验参数的取值为f=1.16；ζ=0.038；C_YO=0.4；h^*=1.24；χ=0.625；C_D=1.2^{[参考文献 13

百度学术}13]。

从脉动频率的物理意义来讲，尾流振子模型所描述的总脉动气动力系数主要包含与柱体后侧旋涡脱落频率同频脉动的部分（记为 $C_{Y} (ω_{K})$ ）、与柱体振动频率同频脉动的部分（记为 $C_{Y} (ω_{N})$ ），以及以某些高次谐波为主的随机脉动部分。其中，只有脉动频率与结构固有频率相近的部分 $C_{Y} (ω_{N})$ 才会对柱体的风致响应产生决定性的影响，其余脉动部分的影响极小，基本可以忽略。所以，着重对与 $C_{Y} (ω_{N})$ 相关的关键参数进行讨论。

对于一个在横风向以固定振幅做简谐振动（即位移 $Y = \hat{Y} s i n (ω_{N} t)$ ）的柱体，气动力系数中与柱体振动位移同频脉动的部分可以表示为

C_{Y} (ω_{N}) = {\hat{C}}_{Y} (ω_{N}) s i n (ω_{N} t + φ)

（2）

式中： $φ$ 为 $C_{Y} (ω_{N})$ 领先于柱体振动位移 $Y$ 的相位角。进一步地，可以将 $C_{Y} (ω_{N})$ 分解为与柱体振动位移Y同相和90 $°$ 异相（即与振动速度 $\dot{Y}$ 同相）的两个部分，并将这两个部分的脉动幅值分别记为 ${\hat{C}}_{Y I} (ω_{N})$ 和 ${\hat{C}}_{Y O} (ω_{N})$ 。它们可以通过式（3）、式（4）与 ${\hat{C}}_{Y} (ω_{N})$ 和 $φ$ 联系在一起。

{\hat{C}}_{Y} (ω_{N}) = \sqrt[]{[{\hat{C}}_{Y I} (ω_{N} {)]}^{2} + [{\hat{C}}_{Y O} (ω_{N} {)]}^{2}}

（3）

φ = a r c t a n [{\hat{C}}_{Y O} (ω_{N}) / {\hat{C}}_{Y I} (ω_{N})]

（4）

以上对脉动气动力系数 $C_{Y} (ω_{N})$ 的分解可以被赋予一种新的物理意义，即与柱体振动位移同相的部分可以被等效地视为气动刚度力贡献的部分，其脉动幅值 ${\hat{C}}_{Y I} (ω_{N})$ 可以被视为等效气动刚度系数；而与振动位移90 $^{\circ}$ 异相（即与振动速度 $\dot{Y}$ 同相）的部分则可以被等效地视为气动阻尼力贡献的部分，其脉动幅值 ${\hat{C}}_{Y O}^{} (ω_{N})$ 则可以被称作等效气动阻尼系数。对于Tamura尾流振子模型，关键参数 ${\hat{C}}_{Y I}^{} (ω_{N})$ 和 ${\hat{C}}_{Y O}^{} (ω_{N})$ 可以由式(5)计算得到。

\{\begin{array}{l} {\hat{C}}_{Y I} (ω_{N}) = - f \hat{α} c o s θ \\ {\hat{C}}_{Y O} (ω_{N}) = - [f (\frac{\hat{Y}}{V_{ω N}}) + C_{D} (\frac{\hat{Y}}{V_{ω N}}) + f \hat{α} s i n θ] \end{array}

（5）

式中： $\hat{α}$ 为尾流振子角位移幅值；θ为尾流振子的角位移α和柱体位移Y的相位差，二者可以由式（6）求得。

\{\begin{array}{l} 4 ζ^{2} Ω_{K N}^{2} • {\hat{α}}^{2} {(1 - f^{2} {\hat{α}}^{2} / {C_{Y O}^{2}}^{2} + (Ω_{K N}^{2} - 1) • {\hat{α}}^{2} - [χ^{2} + (Ω_{K N}^{2} / V_{ω N})^{2}] \hat{Y})}^{2} = 0 \\ θ = a r c t a n [\frac{2 ζ Ω_{K N}^{} (1 - f^{2} {\hat{α}}^{2} / C_{Y O}^{2}) • χ - (Ω_{K N}^{2} - 1) • (Ω_{K N}^{2} / V_{ω N})}{2 ζ Ω_{K N}^{} (1 - f^{2} {\hat{α}}^{2} / C_{Y O}^{2}) • (Ω_{K N}^{2} / V_{ω N}) + (Ω_{K N}^{2} - 1) • χ}] \end{array}

（6）

值得注意的是，对于以高耸烟囱为代表的浸没于空气来流中的工程结构，气动阻尼力贡献的部分能较为明显地影响结构最终的表观阻尼，而气动刚度力贡献部分对整个结构的影响却很微小，基本可以忽略。所以，结构任一断面处单位长度的涡激力可以仅按气动阻尼力部分进行近似考虑，见式(7)。

\begin{array}{l} F_{i} (t) = - \frac{1}{2} ρ U^{2} D \cdot \\ [f (\frac{\hat{Y}}{V_{ω N}}) + C_{D} (\frac{\hat{Y}}{V_{ω N}}) + f \hat{α} s i n θ] \cdot c o s (ω_{N} t) \end{array}

（7）

2.2 三维连续结构振动方程

式(7)仅适用于二维断面，对于三维烟囱结构，需要考虑风速和烟囱断面的变化。烟囱涡振三维模型如图5所示。烟囱断面尺寸随烟囱高度而变化。来流风沿X轴方向作用于高度为H的高耸烟囱结构上，风速沿铅垂高度分布服从幂指数律，即

U (z) = U_{R} {(\frac{z}{z_{R}})}^{φ}

（8）

式中：U_R为参考高度Z_R处的参考风速；U（z）为烟囱任意高度z处平均来流风速；φ为地表粗糙度系数。

图5 三维烟囱涡振示意图

Fig. 5 Schematic diagram of three-dimensional vortex-induced vibration of chimney

假定烟囱结构第一阶弯曲振型的振型函数为φ（z），那么应用振型坐标变换可以求得广义位移y（z,t），即

y (z, t) = φ (z) q (t)

（9）

式中：q（t）为第一阶弯曲振型广义坐标。

引入以下模态质量、模态刚度和模态阻尼：

M = \int_{0}^{H} m (z) φ^{2} (z) d z

（10）

K = ω_{N}^{2} M

（11）

N = 2 ξ \sqrt[]{M \cdot K}

（12）

式中：m（z）为结构质量分布密度；ω为第一阶模态的自振圆频率；ξ为第一阶模态的阻尼比。

可以得到第一阶弯曲模态的运动方程为

M (\ddot{q} + 2 η ω_{N} \dot{q} + ω_{N}^{2} q) = Q (t)

（13）

式中的Q(t)为模态广义力，可由式（14）求得。

Q (t) = \int_{0}^{H} F_{y} (z, t) φ (z) d z

（14）

式中F_y（z,t）为高度z处单位长度塔柱所受到的Y轴向横风向涡激力，由式（7）可以得到

F_{y} (z, t) = \frac{1}{2} ρ U {(z)}^{2} D (z) \hat{C} (z) • c o s (ω_{N} • t)

(15)

式中：D(z)为烟囱断面尺寸随高度z变化的函数； $\hat{C} (z)$ 为气动阻尼系数幅值随高度z变化的函数。

\begin{array}{l} \hat{C} (z) = \\ - [f (\frac{\hat{Y} (z)}{V_{ω N} (z)}) + C_{D} (\frac{\hat{Y} (z)}{V_{ω N} (z)}) + f \hat{α} (z) s i n θ] \end{array}

（16）

若取烟囱10 m高度处的风速为基本风速，记为U_R，则烟囱任意高度处平均风速U(z)为

U (z) = U_{R} {(\frac{z}{10})}^{φ}

（17）

将式(15)和式(17)带入式(14),可以得到模态广义力的表达式

Q (t) = \int_{0}^{H} \frac{1}{2} ρ U_{R} {[(\frac{z}{10})]}^{φ}^{2} D (z) \hat{C} (z) • c o s (ω_{N} t) φ (z) d z

（18）

将式(18)代入式(13)得

\begin{array}{l} M \overset{• •}{q} + 2 η ω_{N} M \overset{•}{q} + ω_{N}^{2} M q = \\ \int_{0}^{H} \frac{1}{2} ρ {(U_{R} {(\frac{z}{10})}^{φ})}^{2} D (z) \hat{C} (z) • c o s (ω_{N} t) φ (z) d z \end{array}

（19）

2.3 基于有限元的迭代计算方法

从式(19)可以看出，直接求解该公式来算得三维烟囱的涡振响应十分困难，为了避免传统数值计算中复杂的求解运算，采用基于有限元的迭代计算方法进行求解。迭代计算方法是利用ANSYS有限元软件和Matlab编程结合的方式实现的。迭代方法的流程如图6所示。

图6 基于有限元的迭代方法流程

Fig. 6 The process of iterative method based on finite element

1）预设参考高度（10 m）处无量纲风速V_ωN和顶部无量纲振幅 $\hat{Y}$ 可能区间，计算整个区间内所有 $\hat{C}$ 值，得到 $\hat{C}$ 数值表；

2）对于每一个确定的风速，假定其顶部无量纲振幅 ${\hat{Y}}_{假定}$ ；

3）通过结构的一阶模态振型计算出各节点的假定无量纲振幅，从表中查找每个节点对应的 $\hat{C}$ 值；

4）由式(15)计算所有节点的涡激荷载时程，导入ANSYS中计算模型实际顶部无量纲振幅 $\hat{Y}$ ，并与 ${\hat{Y}}_{假定}$ 对比；

5）如果 $\hat{Y}$ 与 ${\hat{Y}}_{假定}$ 的误差小于设定阈值，则认为找到此风速下烟囱的涡振振幅的一个解（由于“迟滞”现象的存在，部分风速下可能出现多个解的情况），继续寻找下一个风速的涡振振幅，否则重复步骤2）~4）。

通过以上可以看出，基于有限元的迭代计算方法步骤简单，逻辑清晰，在计算处于三维流场的烟囱的涡振响应时，可以避免求解复杂的耦合方程。

2.4 方法验证

为了检验迭代计算方法的可靠性，以Feng^[

17]风洞试验中的模型为算例，分别用迭代计算方法和Runge-Kutta数值方法^{[参考文献 18

百度学术}18]计算该模型的涡振响应，将两种方法的结果和该模型的风洞试验数据进行对比。

风洞试验模型如图7所示，该模型为一个弹性支承刚性圆柱，模型长0.685 8 m，断面直径0.076 2 m，结构基本参数为：m=0.949 kg，k=4.1 N/m，c=0.004 1 kg/s，圆柱两端分别连接一个并联的弹簧和阻尼。利用迭代计算方法时，弹性支承刚性圆柱模型在ANSYS有限元软件中采用beam4单元建立，每个质量点约束其y方向和z方向的自由度。弹性支承采用combine14单元，一端与圆柱连接，另一端约束全部自由度为固定端。

图7 弹性支承刚性圆柱模型

Fig. 7 A flexibly mounted rigid circular cylinder

结果如图8所示。当无量纲风速在[0.8,1.1]时，出现了所谓的“锁定区间”，模型出现涡激振动。当无量纲风速在[0.9,1.0]时，模型的涡激振动出现了一个特殊的迟滞现象(图中红色箭头所示)，也称为“滞回区间”。在这个区间中，风速从小到大和从大到小两种情况下，模型的锁定频率、锁定区间和幅值响应是不同的。这是因为在这个风速区间中，结构的涡激振动出现了高、低两种幅值稳态振动的可能，结构最终会进入哪种幅值的振动状态取决于模型前一时刻所受的激励。因为迭代法忽略了气动刚度的影响，而Runge-Kutta方法考虑了气动刚度，所以两种理论方法的计算结果存在一些差距。值得注意的是，有限元迭代方法忽略了涡激共振气动力中具有强迫性质的部分而只考虑了具有自激性质的部分。Staubli^[

19]指出，在涡振锁定区间中部，具有自激性质的气动激励部分影响很大，具有强迫性质的部分影响较小；而在涡振锁定区间边缘，具有强迫性质的气动激励部分影响有所增大。这导致迭代法计算的位移与风洞试验结果在涡振锁定区间边缘有所区别。另外，精确确定尾流振子模型(Tamura)中各经验参数需要通过测力试验与流体可视化试验的联合测试，整个测试过程复杂繁琐，故采用Tamura等^{[参考文献 13

百度学术}13]给出的建议取值，该值与Feng^{[参考文献 17

百度学术}17]试验的实际值有所区别，这是迭代法计算结果与风洞试验结果存在差距的主要原因。但总的来说，两种理论计算方法的结果与风洞试验结果在涡振响应峰值处的无量纲位移及对应的无量纲风速已经非常接近，表明迭代法具有足够的精度且运用迭代法计算得到的涡振响应具有可靠性。

图8 弹性支承刚性圆柱涡振响应对比

Fig. 8 Comparison of vortex-induced vibration response of the flexibly mounted rigid circular cylinder

2.5 烟囱涡振响应计算结果

在考虑了烟囱各断面处风速的变化以及断面直径的变化后，烟囱整体是一个处于三维流场的连续系统。基于Tamura尾流振子模型，烟囱的涡激荷载和运动方程的推导见式(8)~式(19)。运用迭代计算方法，在结构阻尼比取0.5%时，烟囱结构的涡振响应如图9所示。

图9 烟囱涡振曲线

Fig. 9 Vortex-induced vibration curve of chimney

从图9的结果可以看出，烟囱涡振响应的最大无量纲位移幅值可达0.082(实际位移幅值1.5 m)，此时发生涡激共振的无量纲风速（10 m高度处）为0.474（实际风速22.9 m/s），这说明烟囱这种高柔结构的涡振响应对结构的安全存在着一定的隐患，应该采取减振措施来减小烟囱的涡振响应。另外，三维流场中烟囱的涡振响应和二维流场中弹性支承刚性圆柱的涡振响应有着相似的特性，都存在着涡振的滞回区间（图中红色虚线段），这种涡振的滞回现象具有一定的研究意义。

图10是有限元迭代法计算的1%阻尼比下烟囱横风向涡振响应结果和按中国《烟囱工程技术标准》（GB/T 50051—2021）^[

20]及美国规范Code Requirements for Reinforced Concrete Chimneys and Commentary（ACI 307-08）^{[参考文献 21

百度学术}21]计算的结果。可以看出，有限元迭代法的结果和按美国规范计算的结果十分接近，而按中国规范计算的结果明显小于按美国规范和有限元迭代法计算的结果。这是因为中美两国规范中对结构的阻尼比规定不同，且烟囱横风向响应对结构的阻尼比十分敏感，中国规范规定结构阻尼比为5%，而美国规范和有限元迭代法的结构阻尼比为1%。中国规范是基于比较简单的卢曼方法，美国规范则主要参考了Vickery等^{[参考文献 10-11}10-11]的模型，Vickery等^{[参考文献 10-11}10-11]的模型是在大量试验和实测数据基础上建立的经验公式，其中很多参数并不具有实际物理意义。有限元迭代法是基于Tamura尾流振子模型，该模型中的参数均具有实际物理意义，较全面地计入了横风向涡激共振中涉及的流-固耦合效应。

图10 烟囱顶部横风向响应对比

Fig. 10 Comparison of cross wind response at chimney top

3 抖振响应精细化预测计算

3.1 频域计算

烟囱的抖振响应频域分析采用的是有限元CQC分析方法，相比于采用单一模态响应进行SRSS组合的传统方法，该方法可以考虑自然风的任意风谱和空间相关性以及结构抖振响应的多模态和模态耦合效应，并且该方法计算效率较高。根据Davenport提出的准定常理论，自然风的紊流在单位长度结构断面上产生的抖振力的公式为^[

\{\begin{array}{l} L_{b} = \frac{1}{2} ρ U^{2} D [2 C_{L} χ_{L u} \frac{u}{U} + (C_{L}^{'} + C_{D}) χ_{L w} \frac{w}{U}] \\ D_{b} = \frac{1}{2} ρ U^{2} D [2 C_{D} χ_{D u} \frac{u}{U} + C_{D}^{'} χ_{D w} \frac{w}{U}] \\ M_{b} = \frac{1}{2} ρ U^{2} D^{2} [2 C_{M} χ_{M u} \frac{u}{U} + C_{M}^{'} χ_{M w} \frac{w}{U}] \end{array}

（20）

式中：C_L为升力系数；C_D为阻力系数；C_M为扭矩系数，这3种静风力系数参考长度均为烟囱断面的直径。 $C_{L}^{'}$ 、 $C_{D}^{'}$ 、 $C_{M}^{'}$ 分别是3种静风力系数对风向角的导数；u为紊流脉动风速的纵向分量；w为紊流脉动风速的横向分量；χ_L_u、χ_L_w、χ_D_u、χ_D_w、χ_M_u、χ_M_w为气动导纳函数。对于烟囱这种圆形断面结构，C_L、C_M、 $C_{L}^{'}$ 、 $C_{D}^{'}$ 、 $C_{M}^{'}$ 可以取0。

由模态叠加，烟囱结构的抖振响应X可表示为

X = Ψ q

（21）

式中：Ψ为正规振型矩阵；q为广义模态坐标矩阵。

在广义模态坐标下，烟囱结构的运动控制方程可以表示为

M \ddot{q} + C \dot{q} + K q = Q_{s e} + Q_{b}

（22）

式中：M为广义质量矩阵；C为广义阻尼矩阵；K为广义刚度矩阵。广义自激力向量Q_se和广义抖振力向量Q_b分别见式(23)、式(24)。

Q_{s e} = \frac{1}{2} ρ U^{2} (A_{s} q + \frac{B}{U} A_{d} \dot{q})

（23）

Q_{b} = \frac{1}{2} ρ U^{2} (A_{b u} \frac{u}{U} + A_{b w} \frac{w}{U})

（24）

式中：A_s为气动刚度矩阵；A_d为气动阻尼矩阵；A_bu和A_bw为结构的总抖振力气动矩阵；u为节点紊流脉动风速沿纵向的向量；w为节点紊流脉动风速沿横向的向量。

式(20)所表示的气动抖振力可以简写为式（25）所示形式。

P_{b} = \frac{1}{2} ρ U (C_{b u} u + C_{b w} w)

（25）

式中：

P_{b} = \{\begin{array}{l} L_{b} \\ D_{b} \\ M_{b} \end{array}\}

C_{b u} = B \{\begin{array}{l} 2 C_{L} \\ 2 C_{D} \\ 2 B C_{M} \end{array}\}

C_{b w} = B c o s θ \{\begin{array}{l} C_{L}^{'} + C_{D} \\ C_{D}^{'} \\ B C_{M}^{'} \end{array}\}

在单元坐标系中，单元节点上的等效抖振力为

Q_{b}^{e} = \int_{L} B^{T} P_{b} d x = \frac{1}{2} ρ U (\int_{L} B^{T} C_{b u} A d x u^{e} + \int_{L} B^{T} C_{b w} A d x w^{e})

= \frac{1}{2} ρ U (A_{b u}^{e} u^{e} + A_{b w}^{e} w^{e})

（26）

A_{b u}^{e} = \int_{L} B^{T} C_{b u} A d x

（27）

A_{b w}^{e} = \int_{L} B^{T} C_{b w} A d x

（28）

式中： $A_{b u}^{e}$ 为纵向脉动风速的单元抖振力气动矩阵； $A_{b w}^{e}$ 为横向脉动风速的单元抖振力气动矩阵；B为插值函数矩阵。

将单元的节点等效抖振力 $Q_{b}^{e}$ 从单元的局部坐标系坐标转换到整体坐标系并进行组装，便可得到总的抖振力气动矩阵。

由随机振动理论，广义模态响应向量q和节点位移向量X的功率谱密度为

S_{q} (ω) = H^{*} (ω) S_{Q b} (ω) H^{T} (ω)

（29）

S_{X} (ω) = Ψ H^{*} (ω) S_{Q b} (ω) H^{T} (ω) Ψ^{T}

（30）

H (ω) = {[- ω^{2} M + i ω (C - \frac{1}{2} ρ U b A_{d}) + K - \frac{1}{2} ρ U^{2} A_{s}]}^{- 1}

（31）

式中：H(ω)为频率响应函数矩阵；H^*(ω)为对频率响应函数矩阵的共轭；H^T(ω)为对频率响应函数矩阵的转置； $S_{Q b} (ω)$ 为广义抖振力的功率谱密度矩阵。

广义抖振力的功率谱密度矩阵可表示为

S_{Q b} (ω) = \frac{1}{4} ρ^{2} U^{2} (A_{b u}^{*} S_{u u} A_{b u}^{T} + A_{b u}^{*} S_{u w} A_{b w}^{T} +

A_{b w}^{*} S_{w u} A_{b u}^{T} + A_{b w}^{*} S_{w w} A_{b w}^{T})

（32）

式中：S_u_u和S_w_w分别为脉动风速向量u和w的功率谱密度矩阵； $S_{u w} = S_{w u}^{*}$ 为脉动风速向量u与w的交叉谱密度矩阵。

由式（29）、式（30），功率谱密度矩阵S_q和S_X中的每一个元素可分别由式（33）、式（34）算得。

S_{q_{i j}}^{} (ω) = \sum_{k = 1}^{m} \sum_{l = 1}^{m} H_{i k}^{*} (ω) S_{Q b k l}^{} (ω) H_{j l}^{} (ω)

（33）

S_{X_{i}}^{} (ω) = \sum_{k = 1}^{m} \sum_{l = 1}^{m} ψ_{i k} S_{q k l}^{} (ω) ψ_{i l}

（34）

进而求得相应的方差为

σ_{q_{i i}}^{2} = \int_{0}^{\infty} S_{q_{i i}}^{} (ω) d ω

（35）

σ_{X_{i}}^{2} = \int_{0}^{\infty} S_{X_{i}}^{} (ω) d ω

（36）

烟囱的抖振计算采用的是Von Karman风谱。以往的研究表明交叉风谱对抖振分析结果的影响并不明显，因此忽略S_u_w和S_w_u的影响。按照美国荷载规范Minimum Design Loads and Associated Criteria for Buildings and Other Structures（ASCE/SEI7-16）^[

22]的D类地貌类型（平原、直接暴露于从开阔水面上吹来的风和无障碍地海岸，包括泥平地、盐碱地和不间断的冰地）选取参数，如表1所示。

表1 抖振响应计算参数

Table 1 Parameters of buffeting response

基本风速/(m/s)	空气密度/(kg/m³)	阻尼比	名义湍流度	阻力系数
57	1.225	0.1%~2%	0.15	0.7
地表粗糙高度	粗糙度系数	湍流高度指数	沿长度尺度变化指数	湍流积分长度尺度/m
0.01	1/9	-0.167	0.125	198.12

3.2 时域计算

Aboshosha等^[

15]提出了一种连续离散流随机生成技术(Consistent Discrete Random Flow Generation：CDRFG)，该方法改进了不同频率内的湍流流速相关性，得到的湍流谱与目标Von Karman谱具有很高的相似度。比较该方法产生的基础力矩、顶部加速度以及风洞试验产生的数据可以发现，该方法与风洞试验结果非常接近。采用CDRFG方法生成了基本风速为57 m/s的烟囱脉动风速时程，时长为10 min，间隔为0.01 s，脉动风速时程和功率谱分别如图11~图14所示，具体参数如表1所示。按照式(20)中的Davenport准定常抖振力公式，计算得到烟囱抖振力时程，进而将抖振力时程输入烟囱ANSYS模型中进行抖振时域计算。

（a） 10 m高度处

（b）烟囱顶部

图11 顺风向脉动风速时程

Fig. 11 Time history of turbulent flow in downwind direction

（a） 10 m高度处

（b）烟囱顶部

图12 顺风向烟囱脉动风功率谱

Fig. 12 Power spectrum density of turbulent flow in downwind direction

（a） 10 m高度处

（b）烟囱顶部

图13 横风向脉动风速时程

Fig. 13 Time history of turbulent flow in crosswind direction

（a） 10 m高度处

（b）烟囱顶部

图14 横风向烟囱脉动风功率谱

Fig. 14 Power spectrum density of turbulent flow in crosswind direction

3.3 抖振响应计算结果

从施工建设到运行使用，烟囱在不同阶段下结构的阻尼存在差别。通过时域分析和频域分析两种方法，考虑了不同风速下（结构阻尼比为0.5%）和不同结构阻尼比下（基本风速为57 m/s）烟囱顶部的抖振响应。因为CDRFG方法生成的脉动风速时程具有随机性，同一基本风速下生成的多个脉动风时程彼此间存在差异，所以每一基本风速下烟囱抖振响应时域计算的结果值存在一定的浮动范围。对此，在参考点每一基本风速（U_R=10、20、30、40、50、60 m/s）下生成多条风速时程进行时域计算，以每个基本风速下时域计算结果置信度为90%的区间和平均值与频域计算结果对比，其结果分别如图15和图16所示。由图可见，时域分析和频域分析的结果具有良好的一致性，且随着风速的增大，该烟囱抖振响应时域计算结果的随机性在逐步增大；而随着该烟囱阻尼比的减小，抖振响应时域计算结果的随机性也在增大。在美国荷载规范Minimum Design Loads and Associated Criteria for Buildings and Other Structures（ASC/SEIE7-16）^[

22]D类地貌条件下，不同基本风速的烟囱顶部的抖振响应都较小。不同结构阻尼比下的结果反映出烟囱的抖振响应对结构阻尼比十分敏感，这间接说明使用调频减振装置来控制烟囱的抖振响应是一种有效的减振措施。另外值得注意的是，虽然缠绕螺旋线（板）等气动措施一般可以有效抑制圆形断面的涡激振动，但运用于实际烟囱时，该类措施在超临界和跨临界雷诺数环境下的减振效果相比于亚临界雷诺数环境将会降低^{[参考文献 23-25}23-25]，以及可能会使烟囱的阻力系数增大，从而增大烟囱抖振响应的风险^{[参考文献 26-27}26-27]。

（a）顺风向

（b）横风向

图15 不同风速下烟囱顶部位移RMS值

Fig. 15 RMS values of the buffeting displacement at the chimney top under different wind speeds

（a）顺风向

（b）横风向

图16 不同阻尼比下烟囱顶部位移RMS值

Fig. 16 RMS values of the buffeting displacement at the chimney top under different damping ratios

图17是按中美两种规范和有限元CQC方法计算的烟囱顺风向响应的结果。可以看出，有限元CQC方法计算的结果和按照中美两种规范计算的结果吻合良好。其中，按美国规范计算的结果较大是因为按照美国规范Code Requirements for Reinforced Concrete Chimneys and Commentary（ACI 307-08）所取的阻力系数在烟囱顶部一段会出现突变，从0.65变为1，而有限元CQC方法的阻力系数统一取0.7。按中美两国规范计算的结果存在差异有两个原因：首先，中美两国规范的计算方法不同，美国规范计算顺风向风荷载的方法是将平均风荷载和脉动风荷载相加。脉动风荷载是通过平均风荷载产生的基础弯矩进一步计算而来。而中国规范计算脉动风荷载的方法是在平均风荷载的基础上乘风振系数。其次，相比美国规范，中国规范中所采用的10 m高度处的湍流强度、峰值因子以及体型系数等都较小。

图17 烟囱顶部顺风向响应对比

Fig. 17 Comparison of downwind response at chimney top

3.4 顺风向和横风向的效应组合

虽然顺风向风荷载、横风向风振等效风荷载一般同时出现，但结合实际情况和工程经验^[

21]，沿顺风向的风振响应与横向风振响应的相关性较小。按美国规范Code Requirements for Reinforced Concrete Chimneys and Commentary（ACI 307-08）将顺风向响应和横风向响应进行组合，烟囱在设计风速(57 m/s)和最大涡激振动幅值对应风速(22.9 m/s)各高度的顺风向弯矩、横风向弯矩和组合弯矩分别如图18、图19所示。可以看出，在设计风速下，由于此时烟囱只有抖振响应，烟囱的顺风向弯矩远大于横风向弯矩，顺风向弯矩对组合弯矩的贡献也较大，烟囱整体的风致响应以顺风向抖振响应为主导；而在最大涡激振动幅值对应风速下，抖振响应和涡振响应同时作用于烟囱，但因为发生了涡激振动，烟囱的横风向弯矩远大于顺风向弯矩，组合弯矩也基本以横风向弯矩的贡献为主，此时烟囱整体的风致响应以横风向涡振响应为主导。

图18 57 m/s风速下烟囱弯矩沿高度变化

Fig. 18 Variation of chimney bending moment along height at wind speed of 57 m/s

图19 22.9 m/s风速下烟囱弯矩沿高度变化

Fig. 19 Variation of chimney bending moment along height at wind speed of 22.9 m/s

4 结论

基于Tamura尾流振子模型，针对圆形断面结构涡振响应计算提出一种新的迭代计算方法。通过与前人试验测得的风洞数据的对比，证明了迭代计算方法的精度和可靠性。基于结构的固有模态坐标，建立了适用于高耸烟囱耦合抖振响应分析的有限元CQC频域计算方法。基于以上两种方法，以某大型烟囱为算例，计算了烟囱结构的横风向涡振响应和顺风向抖振响应。根据上述研究，获得下列结论：

1）有限元迭代计算方法具有足够的精度和可靠性，将Tamura尾流振子模型成功地运用在三维结构上，与中美两国烟囱设计规范相比，更加充分地考虑了发生涡激共振时高耸烟囱的流固耦合效应，可以有效地计算出实际工程结构的涡振响应，为高耸烟囱这类结构横风向涡振响应的计算提供了一种新的视野。

2）从烟囱整体的涡振响应可以看出该算例中的烟囱受涡振影响较大，存在一定安全隐患，需要采用一定的减振措施降低其涡振响应的幅值。三维流场中烟囱的涡振响应存在着涡振的滞回区间，这种涡振的滞回现象具有一定的研究意义。

3）对于横风向响应，运用有限元迭代法计算的结果和按美国规范Code Requirements for Reinforced Concrete Chimneys and Commentary（ACI 307-08）计算的结果较为接近，而按中国《烟囱工程技术标准》（GB/T 50051—2021）计算的结果则明显小于前两者；对于顺风向响应，运用有限元CQC方法计算的结果和按中美两国烟囱设计规范计算的结果均较为接近。

4）虽然顺风向风荷载、横风向风振等效风荷载一般同时出现，但是，结合实际情况和工程经验，沿顺风向的风振响应与横向风振响应的相关性较小。当发生横风向涡振响应时，烟囱整体风致响应以横风向涡振响应为主；而当风速并未处于涡振风速区间时，烟囱整体风致响应以顺风向抖振响应为主。

参考文献

陈鑫. 高耸钢烟囱风振控制理论与试验研究[D]. 南京: 东南大学, 2012. [百度学术]

CHEN X. Theoretical and experimental study on vibration control of high-rise steel chimneys under wind load. Nanjing: Southeast University, 2012. (in Chinese) [百度学术]

刘仰昭. 均匀流中方形柱体横风向气动特性及模型化研究[D]. 成都: 西南交通大学, 2018. [百度学术]

LIU Y Z. Study on aerodynamic characteristics of a transversely oscillating square cylinder in smooth flow and its modeling approach [D]. Chengdu: Southwest Jiaotong University, 2018. (in Chinese). [百度学术]

CHEN C, MANNINI C, BARTOLI G, et al. Experimental study and mathematical modeling on the unsteady galloping of a bridge deck with open cross section [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2020, 203: 104170. [百度学术]

陈政清. 桥梁风工程[M]. 北京: 人民交通出版社, 2005: 64-133. [百度学术]

CHEN Z Q. Wind engineering of bridge [M]. Beijing: China Communications Press, 2005: 64-133. (in Chinese) [百度学术]

张勇, 曹素功, 马如进, 等. 大跨度悬索桥涡激振动动态监控预测[J]. 振动与冲击, 2020, 39(8): 143-150. [百度学术]

ZHANG Y, CAO S G, MA R J, et al. Dynamic monitoring and prediction for vortex induced vibration of a sea crossing bridge [J]. Journal of Vibration and Shock, 2020, 39(8): 143-150. (in Chinese) [百度学术]

廖海黎, 李明水, 马存明, 等. 桥梁风工程2019年度研究进展[J]. 土木与环境工程学报(中英文), 2020, 42(5): 56-66. [百度学术]

LIAO H L, LI M S, MA C M, et al. State-of-the-art review of bridge wind engineering in 2019 [J]. Journal of Civil and Environmental Engineering, 2020, 42(5): 56-66. (in Chinese) [百度学术]

MANNINI C. Incorporation of turbulence in a nonlinear wake-oscillator model for the prediction of unsteady galloping response [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2020, 200: 104141. [百度学术]

LIU Y Z, MA C M, LI Q S, et al. A new modeling approach for transversely oscillating square-section cylinders [J]. Journal of Fluids and Structures, 2018, 81: 492-513. [百度学术]

张亮亮, 吴蕊恒, 倪志军, 等. 全风向角下二维切角方形桥塔气动措施数值模拟[J]. 土木与环境工程学报(中英文), 2019, 41(2): 116-121. [百度学术]

ZHANG L L, WU R H, NI Z J, et al. Numerical simulation on aerodynamic measures of 2D-corner-cutoff square cylinder under all yaw wind angles [J]. Journal of Civil and Environmental Engineering, 2019, 41(2): 116-121. (in Chinese) [百度学术]

VICKERY B J, BASU R I. Across-wind vibrations of structures of circular cross-section. Part I. Development of a mathematical model for two-dimensional conditions [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1983, 12(1): 49-73. [百度学术]

VICKERY B J, BASU R I. The response of reinforced concrete chimneys to vortex shedding [J]. Engineering Structures, 1984, 6(4): 324-333. [百度学术]

HARTLEN R T, CURRIE I G. Lift-oscillator model of vortex-induced vibration [J]. Journal of the Engineering Mechanics Division, 1970, 96(5): 577-591. [百度学术]

TAMURA Y, AMANO A. Mathematical model for vortex-induced oscillations of continuous systems with circular cross section [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1983, 14(1/2/3): 431-442. [百度学术]

TAMURA Y. Mathematical models for understanding phenomena: Vortex-induced vibrations [J]. Japan Architectural Review, 2020, 3(4): 398-422. [百度学术]

ABOSHOSHA H, ELSHAER A, BITSUAMLAK G T, et al. Consistent inflow turbulence generator for LES evaluation of wind-induced responses for tall buildings [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2015, 142: 198-216. [百度学术]

李永乐, 廖海黎, 强士中. 桥梁抖振时域和频域分析的一致性研究[J]. 工程力学, 2005, 22(2): 179-183. [百度学术]

LI Y L, LIAO H L, QIANG S Z. Bridge buffeting analysis in time and frequency domains [J]. Engineering Mechanics, 2005, 22(2): 179-183. (in Chinese) [百度学术]

FENG C C. The measurement of vortex-induced effects on flow past stationary and oscillating circular D-section cylinders [D]. Vancouver: University of British Columbia, 1968. [百度学术]

Mannini C, Massai T, Marra A M, et al. Modelling the interaction of VIV and galloping for rectangular cylinders[C]//International Conference on Wind Engineering - ICWE14, 2015. [百度学术]

STAUBLI T. Calculation of the vibration of an elastically mounted cylinder using experimental data from forced oscillation [J]. Journal of Fluids Engineering, 1983, 105(2): 225-229. [百度学术]

烟囱工程技术标准: GB/T 50051—2021 [S]. 北京: 中国计划出版社, 2021. [百度学术]

Technical standard for chimney engineering: GB/T 50051—2021 [S]. Beijing: China Planning Press, 2021.(in Chinese) [百度学术]

Code requirements for reinforced concrete chimneys and commentary: ACI 307-08 [S]. Farmington Hills, MI: Concrete Institute, 2008. [百度学术]

American Society of Civil Engineers. Minimum design loads and associated criteria for buildings and other structures: ASCE/SEI 7-16 [S]. Reston, VA: American Society of Civil Engineers, 2016. [百度学术]

GARTSHORE I S, KHANNA J, LACCINOLE S. The effectiveness of vortex spoilers on a circular cylinder in smooth and turbulent flow[M]//Wind Engineering. Amsterdam: Elsevier, 1980: 1371-1379. [百度学术]

RUSCHEWEYH H. Straked in-line steel stacks with low mass-damping parameter [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1981, 8(1/2): 203-210. [百度学术]

杜晓庆, 林伟群, 施春林, 等. 高雷诺数下并列双圆柱绕流的大涡模拟[J]. 哈尔滨工业大学学报, 2019, 51(6): 193-200. [百度学术]

DU X Q, LIN W Q, SHI C L, et al. Large eddy simulation of flow around two side-by-side circular cylinders at a high Reynolds number [J]. Journal of Harbin Institute of Technology, 2019, 51(6): 193-200. (in Chinese) [百度学术]

刘庆宽, 卢照亮, 田凯强, 等. 螺旋线对斜拉桥斜拉索高雷诺数风致振动影响的试验研究[J]. 振动与冲击, 2018, 37(14): 175-179. [百度学术]

LIU Q K, LU Z L, TIAN K Q, et al. Experiments on the effect of helical line on the stay-cable vibration at high Reynolds number [J]. Journal of Vibration and Shock, 2018, 37(14): 175-179. (in Chinese) [百度学术]

郑云飞, 刘庆宽, 战启芳, 等. 螺旋线参数对斜拉索气动特性影响的试验研究[J]. 工程力学, 2020, 37(Sup1): 301-306. [百度学术]

ZHENG Y F, LIU Q K, ZHAN Q F, et al. Experimental study on helical line paraments' effect on aerodynamic characteristics of stay cables [J]. Engineering Mechanics, 2020, 37(Sup1): 301-306. (in Chinese) [百度学术]