文章目录
1. 涡激振动简介
涡激振动是结构在一定条件下，由于旋涡周期性脱落带来的横风向周期性激励与结构自振频率一致导致的共振现象。
发生涡激振动需要满足两个条件：
1、一个是在结构高度范围内存在共振区，即结构顶点风速大于共振风速（也称为临界风速）；
2、发生强风共振，即结构处于跨临界区（临界区的定义可以参考任意一本风工程教材）。
本文主要针对混凝土塔架在混凝土段吊装完毕以后，未张拉预应力钢绞线前可能会发生的涡激振动进行计算。
2. EuroCode（欧标）中涡激振动的计算流程
EuroCode中计算涡激振动方法可以参考EuroCode1-4 附录E.1。
2.1 发生涡激共振的判定（检验是否需要进一步计算）
按照EC1-4，满足下式可以不再考虑涡激振动：
v c r i t , i > 1.25 ⋅ v m (1) v_{crit,i}>1.25 \cdot v_{m} \tag{1} vcrit,i>1.25⋅vm(1)
其中， v m v_{m} vm为验算高度处（一般取结构顶部）10min平均风速，风电行业的取值一般由风资源直接提供，不考虑其他修正； v c r i t , j v_{crit,j} vcrit,j为临界风速，它是导致漩涡脱落频率与结构自振频率相等的风速，表达式如下：
v c r i t , i = d ⋅ n i , y S t (2) v_{crit,i}=\frac{d\cdot n_{i,y}}{S_{t}} \tag{2} vcrit,i=Std⋅ni,y(2)
S t S_{t} St——斯托拉哈数，取0.18；
d d d——取3/4高度处塔筒外直径；
n i , y n_{i,y} ni,y——结构 i i i阶固有频率。
2.2 惯性力的计算
本节先看两个无量纲常数——斯科拉顿数和雷诺数，然后给出惯性力及其依赖参数的计算方法。
2.2.1 斯科拉顿数
S c = 2 ⋅ δ s ⋅ m i , e ρ ⋅ d 2 (3) S_{c}=\frac{2\cdot \delta_{s}\cdot m_{i,e}} {\rho \cdot d^2} \tag{3} Sc=ρ⋅d22⋅δs⋅mi,e(3)
δ s \delta_{s} δs——结构阻尼，取值0.03（按照EuroCode1-4 Table F.2取值）；
m i , e m_{i,e} mi,e——模态 i i i的单位长度等效质量；
ρ \rho ρ——空气密度；
2.2.2 雷诺数
R e ( v c r i t , i ) = d ⋅ v c r i t , i ν (4) Re(v_{crit,i}) = \frac{d\cdot v_{crit,i}}{\nu} \tag{4} Re(vcrit,i)=νd⋅vcrit,i(4)
ν \nu ν——空气运动粘滞系数( ν ≈ 1.5 × 1 0 − 6 m 2 / s \nu \approx1.5\times 10^{-6} m^2/s ν≈1.5×10−6m2/s)
2.2.3 单位长度上的惯性力
F w ( s ) = m ( s ) ⋅ ( 2 ⋅ π ⋅ n i , y ) 2 ⋅ ϕ i , y ( s ) ⋅ y F , m a x (5) F_{w}(s) = m(s)\cdot (2\cdot \pi \cdot n_{i,y})^2 \cdot \phi_{i,y}(s) \cdot y_{F,max} \tag{5} Fw(s)=m(s)⋅(2⋅π⋅ni,y)2⋅ϕi,y(s)⋅yF,max(5)
最大位移的计算
y F , m a x d = 1 S t 2 ⋅ 1 S c ⋅ K ⋅ K w ⋅ c l a t (6) \frac {y_{F,max}}{d}=\frac {1}{S^2_{t}} \cdot \frac {1}{S_{c}}\cdot K \cdot K_{w}\cdot c_{lat} \tag{6} dyF,max=St21⋅Sc1⋅K⋅Kw⋅clat(6)
EuroCode提供了两种计算最大位移的方式，这里给出的是普遍性的计算方法，另一方法在本节末尾2.2.5小节中给出。
其中，
c l a t = { c l a t , 0 , v c r i t , i / v m , L j ⩽ 0.83 ( 3 − 2.4 ⋅ v c r i t , i v m , L j ) ⋅ c l a t , 0 , 0.83 < v c r i t , i / v m , L j < 1.25 0 , 1.25 ⩽ v c r i t , i / v m , L j (7) c_{lat}=\left\{\begin{matrix} c_{lat,0} \ , \qquad \qquad \qquad v_{crit,i}/v_{m,Lj}\leqslant 0.83\\ (3-2.4\cdot \frac{v_{crit,i}}{v_{m,Lj}})\cdot c_{lat,0} \ , \qquad 0.83<v_{crit,i}/v_{m,Lj}<1.25\\ 0 \ , \qquad \qquad \qquad \quad \ 1.25\leqslant v_{crit,i}/v_{m,Lj} \end{matrix}\right. \tag{7} clat=⎩⎨⎧clat,0 ,vcrit,i/vm,Lj⩽0.83(3−2.4⋅vm,Ljvcrit,i)⋅clat,0 ,0.83<vcrit,i/vm,Lj<1.250 , 1.25⩽vcrit,i/vm,Lj(7)
对于塔架（圆柱筒结构） c l a t , 0 c_{lat,0} clat,0根据雷诺数不同取值如下：
a) c l a t , 0 c_{lat,0} clat,0 取值
v m , L j v_{m,Lj} vm,Lj是相关长度中心处的平均风速，相关长度定义如下：
b) 相关长度定义图示
相关长度 L j L_{j} Lj取值如下：
L j / d = { 6 , y F , m a x / d < 0.1 4.8 + 12 ⋅ y F , m a x d , 0.1 ⩽ y F , m a x / d ⩽ 0.6 12 , 0.6 < y F , m a x / d (8) L_{j}/d=\left\{\begin{matrix} 6, \qquad\qquad\quad\quad y_{F,max}/d<0.1\\ 4.8 + 12\cdot \frac{y_{F,max}}{d}, \qquad 0.1\leqslant y_{F,max}/d\leqslant 0.6\\ 12, \qquad\quad\quad\quad\quad 0.6<\quad y_{F,max}/d \end{matrix}\right.\tag{8} Lj/d=⎩⎨⎧6,yF,max/d<0.14.8+12⋅dyF,max,0.1⩽yF,max/d⩽0.612,0.6<yF,max/d(8)
因此，需要先对 y F , m a x y_{F,max} yF,max进行假设，然后再试算，直到得出符合假设的 y F , m a x y_{F,max} yF,max。
2.2.4 有效长度相关系数（ K w K_{w} Kw）与振型系数（ K K K）
K w = ∑ j = 1 n ∫ L j ∣ ϕ i , y ( s ) ∣ d s ∑ j = 1 m ∫ l j ∣ ϕ i , y ( s ) ∣ d s ⩽ 0.6 (9) K_{w}=\frac{\sum_{j=1}^{n}\int_{L_{j}} \left | \phi_{i,y}(s) \right |ds}{\sum_{j=1}^{m}\int_{l_{j}}\left | \phi_{i,y}(s) \right |ds}\leqslant 0.6 \tag{9} Kw=∑j=1m∫lj∣ϕi,y(s)∣ds∑j=1n∫Lj∣ϕi,y(s)∣ds⩽0.6(9)
注意 ( 9 ) (9) (9)式中分子分母积分域并不相同，分子是在相关长度上进行积分，分母是在整个结构上进行节分。
K = ∑ j = 1 m ∫ l j ∣ ϕ i , y ( s ) ∣ d s 4 π ⋅ ∑ j = 1 m ∫ l j ϕ i , y 2 ( s ) d s (10) K=\frac{\sum_{j=1}^{m}\int_{l_{j}} \left | \phi_{i,y}(s) \right |ds}{4\pi \cdot \sum_{j=1}^{m}\int_{l_{j}} \phi^{2}_{i,y}(s)ds} \tag{10} K=4π⋅∑j=1m∫ljϕi,y2(s)ds∑j=1m∫lj∣ϕi,y(s)∣ds(10)
对于塔架（悬臂结构），其 K w K_{w} Kw和 K K K的取值如下：
K w = 3 ⋅ L j / d λ ⋅ [ 1 − L j / d λ + ( L j / d λ ) 2 / 3 ] (11) K_{w}=3\cdot \frac{L_{j}/d}{\lambda}\cdot\left [ 1-\frac{L_{j}/d}{\lambda}+ (\frac{L_{j}/d}{\lambda})^2/3 \right ] \tag{11} Kw=3⋅λLj/d⋅[1−λLj/d+(λLj/d)2/3](11)
其中， λ = l / d \lambda=l/d λ=l/d。
K = 0.13 (12) K=0.13 \tag{12} K=0.13(12)
2.2.5 最大位移的另一种计算方法
先挖个坑，回头慢慢填
2.3 塔架结构固有频率及振型的计算
上一节惯性力的计算中需要提供结构的频率以及振型，本节给出塔架的频率及振型的计算，相关内容可以参考EuroCode1-4 附录F。结构频率以及振型的求解也可以借助有限元方法。
2.3.1 自振频率的计算
n 1 = ε 1 ⋅ d h e f f 2 ⋅ W s W t (13) n_{1}=\frac{\varepsilon_{1}\cdot d}{h^{2}_{eff}}\cdot \sqrt{\frac{W_{s}}{W_{t}}} \tag{13} n1=heff2ε1⋅d⋅WtWs (13)
h e f f = h 1 + h 2 3 (14) h_{eff}=h_{1}+\frac{h_{2}}{3} \tag{14} heff=h1+3h2(14)
相关几何参数取值如下图：
c) 几何参数取值图示
注意：图中 h 3 = h 1 / 3 h_{3}=h_{1}/3 h3=h1/3 。
上式中d为塔架顶部的直径， W s W_{s} Ws是对刚度有贡献的（塔架）结构部分的质量， W t W_{t} Wt是（塔架）结构总质量， ε 1 \varepsilon_{1} ε1的取值为1000（钢制塔架）或者700（混凝土塔架）。
本文主要针对混凝土塔架在混凝土段吊装完毕以后，预应力张拉之前的施工阶段的计算，所以 ε 1 \varepsilon_{1} ε1的取值主要为700。钢制塔架的涡激振动计算可以直接按照DIN EN4133进行计算。
2.3.2 振型的计算
对于塔架，其一阶振型可按下式计算：
ϕ 1 ( z ) = ( z h ) ζ (15) \phi_{1}(z)=(\frac{z}{h})^{\zeta} \tag{15} ϕ1(z)=(hz)ζ(15)
其中， ζ \zeta ζ取2.0 ， h h h为结构高度。
2.3.3 等效质量的计算
m e = ∫ 0 l m ( s ) ⋅ ϕ 1 2 ( s ) d s ∫ 0 l ϕ 1 2 ( s ) d s (16) m_{e}=\frac{\int_{0}^{l}m(s)\cdot \phi ^2_{1}(s)ds}{\int_{0}^{l} \phi ^2_{1}(s)ds}\tag{16} me=∫0lϕ12(s)ds∫0lm(s)⋅ϕ12(s)ds(16)
对于悬臂结构，规范建议可以只计算上部1/3高度（图 c) 中 h 3 h_{3} h3）的质量。
2.4 荷载循环次数的计算
本节就一个公式，如下：
N = 2 T ⋅ n y ⋅ ε 0 ⋅ ( v c r i t v 0 ) 2 ⋅ e x p ( − ( v c r i t v 0 ) 2 ) (17) N = 2T\cdot n_{y}\cdot \varepsilon_{0}\cdot (\frac{v_{crit}}{v_{0}})^2\cdot exp(-(\frac{v_{crit}}{v_{0}})^2)\tag{17} N=2T⋅ny⋅ε0⋅(v0vcrit)2⋅exp(−(v0vcrit)2)(17)
根据2.3节算出的惯性力可以进一步求得涡激振动引起的应力变化，进而根据SN曲线得出在给定风工况下涡激振动对应的需用循环次数 [ N ] [N] [N]，若 N ⩽ [ N ] N\leqslant [N] N⩽[N]则涡激振动满足要求，否则不满足。
d) 混凝土SN曲线
上面的混凝土SN曲线是根据ModelCode1990得到的，当前ModelCode已经更新到2010版本，但是2010版本没有明确指出 S c d , m i n < 0 S_{cd,min}<0 Scd,min<0时的曲线，所以可以暂且采用1990版本给出的结果。
从一般性体型信息文件中提取计算涡激振动所需参数的代码如下（未完待续）：
import openpyxl as opx
import numpy as np
def write_VIV_inp_para(filepath, sheetname, sec_num, rho_c=2650):
wb = opx.load_workbook(filepath)
ws = wb[sheetname]
part_num = sec_num - 1 # 混凝土段的节段数量
H = float(ws.cell(part_num + 3, 2).value)
H_s = list() # 质点高度组成的列表
Delta_s = list() # 积分节点
M_s = list() # 质点质量
Phi_s = list() # 振型
for i in range(part_num):
hc = 0.5 * (float(ws.cell(i + 3, 1).value) + float(ws.cell(i + 3, 2).value))
H_s.append(hc)
height = float(ws.cell(i + 3, 2).value) - float(ws.cell(i + 3, 1).value)
do_bot = float(ws.cell(i + 3, 3).value)
do_top = float(ws.cell(i + 3, 4).value)
di_bot = do_bot - 2 * float(ws.cell(i + 3, 5).value)
di_top = do_bot - 2 * float(ws.cell(i + 3, 6).value)
Vol_o = 3.1416 / 12 * height * (do_bot ** 2 + do_top ** 2 + do_top * do_bot)
Vol_i = 3.1416 / 12 * height * (di_bot ** 2 + di_top ** 2 + di_top * di_bot)
Vol = Vol_o - Vol_i
m_s = rho_c * Vol
M_s.append(m_s)
phi_s = (hc/H)**2
Phi_s.append(phi_s)
Delta_s.append((float(ws.cell(i + 3, 1).value), float(ws.cell(i + 3, 2).value)))
return H_s, M_s, Phi_s, Delta_s
def calcu_me(M_s, Phi_s, Delta_s):
ds = np.array([(i[1]-i[0]) for i in Delta_s])
M_s_np = np.array(M_s)
Phi_s_np = np.array(Phi_s)
me = sum(M_s_np*Phi_s_np**2*ds) / sum(Phi_s_np**2*ds)
return me
3.GB（国标）中涡激振动的计算流程
国标的计算相对简单，先把公式列下来，回头有空慢慢码字
w k = β z μ s μ z w 0 w_{k}=\beta_{z} \mu_{s} \mu_{z} w_{0} wk=βzμsμzw0
β z = 1 + 2 g I 10 B z 1 + R 2 \beta_{z}=1+2gI_{10}B_{z}\sqrt{1+R^2} βz=1+2gI10Bz1+R2
R = π 6 ζ 1 x 1 2 ( 1 + x 1 2 ) 4 / 3 R=\sqrt{\frac{\pi}{6\zeta_{1}} \frac{x^{2}_{1}}{(1+x^{2}_{1})^{4/3}}} R=6ζ1π(1+x12)4/3x12
x 1 = m a x { 30 T 1 k w w 0 , 5 } x_{1} = max\left \{ \frac{30T_{1}} {\sqrt{k_{w}w_{0}}},5 \right \} x1=max{kww0 30T1,5}
T 1 = 0.41 + 0.001 z 2 d T_{1}=0.41+0.001\frac {z^2}{d} T1=0.41+0.001dz2
B z = k H α 1 ρ x ρ y ϕ 1 ( z ) μ z B_{z}=kH^{\alpha_{1} }\rho_{x} \rho_{y} \frac{\phi_{1}(z)}{\mu_{z}} Bz=kHα1ρxρyμzϕ1(z)
ρ x = 10 H + 60 e − H / 60 − 60 H \rho_{x}=\frac{10\sqrt{H+60e^{-H/60}-60}}{H} ρx=H10H+60e−H/60−60
ρ y = 10 B + 60 e − H / 60 − 60 B \rho_{y}=\frac{10\sqrt{B+60e^{-H/60}-60}}{B} ρy=B10B+60e−H/60−60
R e = 69000 ν D Re=69000\nu D Re=69000νD
ν c r = D T i S t \nu_{cr}=\frac{D}{T_{i}S_{t}} νcr=TiStD
ν H = 2000 μ H w 0 ρ \nu_{H}=\sqrt{\frac{2000\mu_{H}w_{0}}{\rho}} νH=ρ2000μHw0
H 1 = H ( ν c r 1.2 ν h ) 1 α H_{1}=H(\frac{\nu_{cr}}{1.2\nu_{h}})^{\frac{1}{\alpha}} H1=H(1.2νhνcr)α1
H 2 = H ( 1.3 ν c r ν H ) 1 α H_{2}=H(\frac{1.3\nu_{cr}}{\nu_{H}})^{\frac{1}{\alpha}} H2=H(νH1.3νcr)α1
λ j = λ j ( H 1 / H ) − λ j ( H 2 / H ) \lambda_{j}=\lambda_{j}(H_{1}/H)-\lambda_{j}(H_{2}/H) λj=λj(H1/H)−λj(H2/H)
w l k = λ j ν c r 2 ϕ j ( z ) 12800 ζ j w_{lk}=\lambda_{j} \nu^{2}_{cr} \frac{\phi_{j}(z)}{12800\zeta_{j}} wlk=λjνcr212800ζjϕj(z)
w = ( 0.6 w k ) 2 + w l k 2 w=\sqrt{(0.6w_{k})^2+w^{2}_{lk}} w=(0.6wk)2+wlk2
