CN103955555B - 一种用于防风抗振的高耸塔器疲劳寿命设计方法 - Google Patents

Abstract

本发明属于高耸塔器的结构设计与安全评定领域，具体涉及一种用于防风抗振的高耸塔器疲劳寿命设计方法。其具体包括以下步骤：设定参考风速，制定离散化的风速风向联合概率分布区间；设计危险截面部位的多种结构形式，并求解其应力集中因子；模拟单位时间内高耸塔器在各参考风速下危险截面不同方向上的名义应力时程响应；计算高耸塔器在危险截面部位不同结构形式和缺陷评定等级下的风致疲劳寿命；判定其寿命是否满足要求，若满足则确定设计方案，否则需重新设计高耸塔器整体结构，并再次进行上述过程。本发明提供的方法可为沿海多风地区高耸塔器的合理设计、安全运行提供借鉴。

Description

一种用于防风抗振的高耸塔器疲劳寿命设计方法

技术领域

本发明属于高耸塔器的结构设计与安全评定领域，具体涉及一种用于防风抗振的高耸塔器疲劳寿命设计方法。

技术背景

塔器是化工、石油化工和炼油等生产中最重要的单元操作设备之一，它的性能好坏，对于整个装置的生产能力、三废处理和环境保护等各个方面都有较大影响。风载荷是塔器较为敏感的载荷之一，它能引起塔器沿着风向的顺风向振动以及垂直于风向的横风向振动，而过于频繁地顺风向或横风向振动又可能会导致其结构疲劳开裂甚至疲劳失效。对于大型直立高耸塔器而言，如脱甲烷塔、精馏塔、焦炭塔，虽然在传统设计中已经考虑了风载荷的影响，但是对沿海等地区多风环境下高耸结构可能会发生风致疲劳的潜在危害性考虑不周。当一些高耸塔器在实际使用中因结构设计不合理产生了局部高应力集中区，或者在制造时本身存在先天缺陷，或者在腐蚀性介质的长期作用下产生的各类新生缺陷，这些都会使高耸塔器的承载能力大幅度降低，使得高耸塔器在长期的风致振动作用下产生疲劳裂纹，最终造成塔器开裂、易燃易爆或有毒有害气体泄漏，给企业的生产带来了安全隐患，近年来高耸塔器由此发生的事故也对相关企业造成了较大的经济损失。

然而国内外的现行标准规范，如美国土木工程师协会（ASCE）标准、我国建筑行业的GB50135《高耸结构设计规范》以及压力容器设计标准JB/T4710《钢制塔式容器》，虽然都对该类结构的风载荷及风振响应给出了明确的计算方法，但是这些方法都将风载荷视为短期载荷，尚未考虑多风环境下高耸塔器的疲劳问题，且这种将顺风向载荷等效为静载荷以及仅考虑横风向临界共振响应的传统方法，并不适合对高耸塔器进行疲劳寿命设计。同时，高耸塔器所在地区的风速风向都有一定的分布规律，这无疑会影响高耸塔器的风致疲劳寿命结果，而随着塔器高度的增加与高强材料的使用，其风敏感性问题也更加突出。

从目前国内外有关研究来看，针对高耸塔器疲劳寿命设计方面的研究几乎还未涉足，此方面的研究多集中在诸如桅杆、电视塔、高楼等结构的风致疲劳寿命分析领域。此外，高耸塔器作为一种大型的石油化工设备，不仅外部有操作平台、扶梯等特殊结构，而且内部还有大量的工艺反应介质，因而其风致疲劳特性及设计要求也与上述结构存在较大的差别。为确保多风地区高耸塔器的运行安全，如何为该类结构提供适合的疲劳寿命校核及设计方法，现已成为一个迫切的需求。

发明内容

针对现有技术的不足，本发明提供了一种用于防风抗振的高耸塔器疲劳寿命设计方法，该方法可为沿海地区高耸塔器的合理设计、安全运行提供参考。

为实现上述目的，本发明采用了以下技术方案：一种用于防风抗振的高耸塔器疲劳寿命设计方法，包括以下步骤：

1、一种用于防风抗振的高耸塔器疲劳寿命设计方法，其包括如下步骤：

S1、以离地面10m高度处的平均风速作为参考风速，设定参考风速为V_r，r＝1，2，…，13；最小参考风速V₁取4m/s，最大参考风速V₁₃取28m/s，中间相邻的参考风速取值间隔为2m/s。

根据高耸塔器所在的地区，选择风速风向联合分布函数，所述风速风向联合分布函数的表达式如下：

$P(u<U_r,{\mathrm{\&theta;}}_j)=f_j\&times;\{1-\exp [-{(U_r/a_j)}^{r_j}\left]\right\},j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,16$

式中，u为风速风向联合分布函数的风速自变量，单位为m/s；

U_r为风速风向联合分布函数的风速值，单位为m/s；

θ_j为0°到360°之间平均分布的16个风向角，且按顺时针依次增大，θ_j＝（j－1)×22.5°，各个风向角θ_j代表不同的风向：如θ₁＝0°代表正北（N）风向，θ₅＝90°代表正东（E）方向，θ₉＝180°代表正南（S）风向，θ₁₁＝225°代表西南（SW）风向；

f_j代表风向角θ_j的频度参数，a_j、r_j分别是风向角θ_j上风速分布的尺度参数、形状参数，f_j、a_j、r_j采用下表中的数据：

风速风向联合分布函数中f_j、a_j、r_j的默认参数值

j	风向	fj	aj	rj		j	风向	fj	aj	rj
											1	N	0.123	4.724	1.592		9	S	0.050	4.374	1.267
2	NNE	0.084	3.222	1.160		10	SSW	0.050	5.140	1.927
											3	NE	0.063	2.630	1.073		11	SW	0.035	4.215	1.448
4	ENE	0.031	2.057	0.822		12	WSW	0.016	3.868	1.377
											5	E	0.076	3.975	1.033		13	W	0.006	3.810	1.315
6	ESE	0.168	3.154	1.009		14	WNW	0.003	3.622	1.596
											7	SE	0.147	2.351	0.908		15	NW	0.023	3.597	1.133
8	SSE	0.055	3.923	1.380		16	NNW	0.069	4.347	1.593

将风速风向联合分布函数离散成与各参考风速V_r相关的多个连续区间，使得风速风向联合分布函数的每个区间[U_r，U_r+1)的中值为V_r，r＝1，2，…，13；参考风速V_r在第j个风向角θ_j上的概率Pr_j为

P_rj＝P(u＜U_r+1,θ_j)-P(u＜U_r,θ_j)

S2、利用ANSYS软件，建立高耸塔器壳单元有限元模型，在塔顶施加弯矩进行静力学求解，高耸塔器轴向应力最大值所在的水平截面为高耸塔器的危险截面；

S3、按照无损检测及容许缺陷要求，设计若干种达到无损检测及容许缺陷要求的危险截面部位结构形式，并分别建立危险截面部位的局部实体单元ANSYS有限元模型，求解危险截面部位在上述若干种结构形式下的应力集中因子；

S4、模拟单位时间内高耸塔器在各参考风速V_r下的危险截面部位的不同方向上的名义应力时程响应，包括如下步骤：

1）、模拟风载荷

（1）将高耸塔器沿竖直方向划分为n段塔节，划分标准参考高耸塔器的设计文件，确定每段塔节的长度l_h、外直径D_h、等效直径D_eh以及每段塔节的中心点距地面高度z_h，其中h＝1，2，3，…，n；各段塔节的中心点即风载荷的模拟点，l_h、D_h、D_eh和z_h的单位为m；

（2）计算各参考风速V_r下各模拟点的平均风速V_rh，其表达式如下：

$V_{\mathrm{rh}}=V_r{\left(\frac{z_h}{10}\right)}^{0.16},h=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n$

（3）计算各参考风速V_r下的顺风向风载荷时程样本

首先计算参考风速V_r下的顺风向脉动风速谱矩阵中的上标a代指顺风向，其表达式如下

式中，矩阵元素 $S_{r{h}_1{h}_2}^a\left(\mathrm{\&omega;}\right)=V_r^2\frac{0.0043X_r^2}{\mathrm{\&pi;\&omega;}{(1+X_r^2)}^{4/3}}\exp [-\frac{28\mathrm{\&pi;\&omega;}|z_{h_1}-z_{h_2}|}{V_{{\mathrm{rh}}_1}+V_{{\mathrm{rh}}_2}}],h_1/h_2=1,$ $2,\&CenterDot;\&CenterDot;\&CenterDot;,n;$

n为高耸塔器的塔节数；

X_r为脉动风速谱的无量纲莫宁坐标，

ω为圆频率，单位为rad/s；

为第h₁、h₂段塔节的中心点距地面高度，即第h₁、h₂个模拟点的高度，单位为m；

为参考风速V_r下第h₁、h₂个模拟点的平均风速，单位为m/s；

则各参考风速V_r下的顺风向脉动风速时程样本由谐波叠加法获得：

式中，h＝1,2,…,n；

Δω＝(ω_up－ω_down)/N；

N为频率等分数，N≥1024；

ω_up为截止频率，ω_down为起始频率，单位为rad/s；

为均匀分布于0～2π之间的随机相位角，单位为rad/s；

ω_l＝ω_down＋(l－1/2)Δω，l＝1，2，…，N；

为的下三角矩阵第h行、m列元素，单位为m·s^1/2/rad^1/2；

时间t的增量Δt必须满足：Δt≤π/ω_up，t＝Δt，2Δt，3Δt，…，T，单位为s；

T为单位时间，T为增量Δt的整数倍，即每次模拟样本的采样时间，单位为s；

于是，各参考风速V_r下的顺风向风载荷时程样本可写成如下形式：

$P_{\mathrm{rh}}^a\left(t\right)=0.5{\mathrm{\&mu;}}_s\mathrm{\&rho;}{D}_{\mathrm{eh}}{l}_h{[V_{\mathrm{rh}}+v_{\mathrm{rh}}^a\left(t\right)]}^2,h=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n$

式中，为高耸塔器第h段塔节的顺风向风载荷时程样本，单位为牛顿（N）；

μ_s为体型系数，μ_s＝0.7；

ρ为空气密度，取ρ＝1.25Kg/m³;

l_h为高耸塔器第h段塔节的长度，单位为m；

D_eh为高耸塔器第h段塔节的等效直径，单位为m；

（4）计算各参考风速V_r下的横风向风载荷时程样本

首先计算参考风速V_r下的横风向风载荷谱矩阵中的c代指横风向，其表达式如下：

式中，矩阵元素 $S_{{\mathrm{rh}}_1{h}_2}^c\left(\mathrm{\&omega;}\right)=S_r^{c1}(z_{h_1},z_{h_2},\mathrm{\&omega;})+S_r^{c2}(z_{h_1},z_{h_2},\mathrm{\&omega;}),h_1/h_2=\mathrm{1,2,3},\&CenterDot;\&CenterDot;\&CenterDot;,$ $n;$

n为高耸塔器的塔节数；

ω为圆频率，单位为rad/s；

为第h₁、h₂个模拟点的高度，单位为m；

为参考风速V_r下的旋涡脱落力互谱密度，单位为N²·s/rad；

为参考风速V_r下的来流紊流力互谱密度，单位为N²·s/rad；

由以下二式获得：

$S_r^{c1}(z_{h_1},z_{h_2},\mathrm{\&omega;})=\sqrt{S_r^{c1}(z_{h_1},\mathrm{\&omega;})S_r^{c1}(z_{h_2},\mathrm{\&omega;})}\cos \left(\frac{2}{3}R\right)\exp (-\frac{1}{3}{R}^2),h_1,h_2=\mathrm{1,2,3},\&CenterDot;\&CenterDot;\&CenterDot;,n;$

$S_r^{c1}(z_h,\mathrm{\&omega;})={\left[\frac{1}{2}\mathrm{\&rho;}{D}_h{V}_{\mathrm{rh}}^2\right]}^2\frac{{\mathrm{\&sigma;}}_{\mathrm{CL}}^2}{2\mathrm{\&pi;}\sqrt{\mathrm{\&pi;}}B{\mathrm{\&omega;}}_{\mathrm{sh}}}\exp \{-{[(1-\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_{\mathrm{sh}}})/B]}^2\},h=\mathrm{1,2,3},\&CenterDot;\&CenterDot;\&CenterDot;,n;$

式中，为参考风速V_r下的旋涡脱落力自谱密度，单位为N²·s/rad；；

R为距离参数， $R=2|z_{h_1}-z_{h_2}|/(D_{h_1}+D_{h_2});$

ρ为空气密度，取ρ＝1.25Kg/m³;

σ_CL为升力系数的方差，取0.14；

B是谱曲线带宽经验参数，取0.18；

ω_sh为高耸塔器h段塔节的旋涡脱落圆频率，ω_sh＝V_rh/(10πD_h)，单位为rad/s；；

分别为第h₁、h₂段塔节的外直径，单位为m；

由以下两式获得：

$S_r^{c2}(z_{h_1},z_{h_2},\mathrm{\&omega;})=\sqrt{S_r^{c2}(z_{h_1},\mathrm{\&omega;})S_r^{c2}(z_{h_2},\mathrm{\&omega;})}\exp \left(-18.7\frac{\mathrm{\&omega;\&pi;}|z_{h_1}-z_{h_2}|}{V_{{\mathrm{rh}}_1}+V_{{\mathrm{rh}}_2}}\right),h_1,h_2=\mathrm{1,2,3},\&CenterDot;\&CenterDot;\&CenterDot;,n;$

$S_r^{c2}(z_h,\mathrm{\&omega;})={\left[\frac{1}{2}\mathrm{\&rho;}{C}_{\mathrm{rh}}{D_{h}V}_{\mathrm{rh}}\right]}^2\frac{15{\mathrm{\&mu;}}_{*r}^2{X}_{\mathrm{rh}}^2}{2\mathrm{\&pi;\&omega;}{(1+9.5X_{\mathrm{rh}}^c)}^{5/3}},h=\mathrm{1,2,3},\&CenterDot;\&CenterDot;\&CenterDot;,n;$

式中，为来流紊流谱的无量纲莫宁坐标,

为来流紊流力自谱密度，单位为N²·s/rad；

C_rh为阻力系数, $C_{\mathrm{rh}}=\left\{\begin{array}{cc}1.2& {\mathrm{Re}}_{\mathrm{rh}}\&le;3\&times;{10}^5\\ 0.2& {\mathrm{Re}}_{\mathrm{rh}}=4\&times;{10}^5\\ 0.6& {\mathrm{Re}}_{\mathrm{rh}}={10}^7\end{array}\right.,$ 中间值按线性插值计算；

Re_rh为雷诺数，Re_rh＝69000V_rhD_h；

u_*r为摩擦速度，u_*r＝0.064V_r，单位为m/s；

则各参考风速V_r下的横风向单位长度的风载荷时程样本由谐波叠加法公式获得：

式中，h＝1,2,…,n；

Δω＝(ω_up－ω_down)/N；

N为频率等分数，N≥1024；

ω_up为截止频率，ω_down为起始频率，单位为rad/s；

为均匀分布于0～2π之间的随机相位角，单位为rad/s；

ω_l＝ω_down＋(l－1/2)Δω，l＝1，2，…，N；

为的下三角矩阵第h行、m列元素，单位为m·s^1/2/rad^1/2；

T为单位时间，即每次模拟样本的采样时间，单位为s；

时间t的增量Δt必须满足：Δt≤π/ω_up，t＝Δt，2Δt，3Δt，…，T，T为增量Δt的整数倍，单位为s；

于是，各参考风速Vr下各点横风向风载荷时程样本为：

$P_{\mathrm{rh}}^c\left(t\right)=l_h{f}_{\mathrm{rh}}^c\left(t\right),h=\mathrm{1,2,3},\&CenterDot;\&CenterDot;\&CenterDot;,n;$

式中，为高耸塔器第h段塔节的横风向风载荷时程样本，单位为牛顿（N）；

l_h为高耸塔器第h段塔节长度，单位为m；

2）、模拟名义应力时程响应

利用ANSYS有限元软件建立高耸塔器梁单元有限元模型，分别加载各参考风速V_r下的顺风向风载荷时程样和横风向风载荷时程样本进行FULL瞬态动力学模块求解；利用ANSYS软件APDL编程，提取各参考风速下危险截面的顺风向名义应力时程响应和横风向名义应力时程响应

求解高耸塔器顺风向风振响应时，高耸塔器的阻尼比为ζ^a＝0.01；

求解高耸塔器横风向风振响应时，参考风速V_r下高耸塔器的阻尼比为

式中，为求解横风向风振响应时高耸塔器的阻尼比；

M^*为高耸塔器的广义质量，

M_h为高耸塔器第h段塔节的质量，由设计文件给出，单位为Kg；

l_h为高耸塔器第h段塔节的长度，单位为m；

为高耸塔器第h段塔节振型函数值，计算公式如下：

H为高耸塔器总高度，单位为m;

A_rh为高耸塔器第h段塔节的气动力系数；

$A_{\mathrm{rh}}=\left\{\begin{array}{cc}0& V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}}<0.85\\ a_h(3.5V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}}-2.95)& 0.85\&le;V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}}<1.0\\ 0.55a_h& 1.0\&le;V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}}<1.1\\ a_h(2.75-2V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}})& \mathrm{1,1}\&le;V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}}<1.3\\ a_h(0.46-0.25V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}})& \mathrm{1,3}\&le;V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}}<1.84\\ 0& 1.84\&le;V_{\mathrm{rh}}/V_{\mathrm{rh}}^{\mathrm{cr}}\end{array}\right.$

为第h段塔节的临界风速，单位为m/s；

n_s为高耸塔器的一阶自振频率，可由ANSYS软件模态求解得出；

a_h为第h段塔节的气动力参数，a_h＝a_1ha_2ha_3h，h＝1，2，3，…，n；

a_1h、a_2h、a_3h分别由以下三式确定：

$a_{1h}=\left\{\begin{array}{cc}1.0& {\mathrm{Re}}_{\mathrm{rh}}<{10}^4\\ 1.8& {10}^4\&le;{\mathrm{Re}}_{\mathrm{rh}}={10}^5\\ 1.0& {{10}^5\&le;\mathrm{Re}}_{\mathrm{rh}}\end{array}\right.$

$a_{2h}=\left\{\begin{array}{cc}1.8& V_r\&le;12m/s\\ 0.9& V_r>12m/s\end{array}\right.$

$a_{3h}=\left\{\begin{array}{cc}1.0& z_h/D_h>12.5\\ 1.0-0.04[12.5-z_h/D_h& z_h/D_h<12.5\end{array}\right.$

Re_rh为雷诺数，Re_rh＝69000V_rhD_h；

3）、叠加名义应力时程响应

当参考风速为Vr、风向角为θ_j时，高耸塔器的危险截面在方向θk（风向角的角度同时也是危险截面上的方向位置）上单位时间T内的名义应力时程响应σ_rjk(t)为

${\mathrm{\&sigma;}}_{\mathrm{rjk}}\left(t\right)={\mathrm{\&sigma;}}_r^a\left(t\right)\cos ({\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_k)+{\mathrm{\&sigma;}}_r^c\left(t\right)\sin ({\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_k)$

式中，为顺风向名义应力时程响应，单位为MPa；

为横风向名义应力时程响应，单位为MPa；

S5、计算高耸塔器在步骤S3中不同结构形式和缺陷评定等级下的风致疲劳寿命，包括如下步骤：

1）、统计应力循环

（1）将步骤S₄中的叠加后的名义应力时程响应σ_rjk(t)修正为

σ_irjk(t)＝K_i·σ_rjk(t)，i＝1,2,…,p

式中，K_i为步骤S3中所求的第i种危险截面部位结构形式的应力集中因子；

p为S3中设计的危险截面部位结构形式个数；

σ_irjk(t)为修正后的应力时程响应，单位为MPa；

（2）利用雨流法统计σ_irjk(t)对应的应力循环；

2）、计算疲劳损伤

根据不同的焊接形式和焊接缺陷所对应的疲劳曲线计算单位时间T内上述应力循环对应的疲劳损伤D_irjk；

3）、计算疲劳寿命

结合离散后的风速风向联合分布，计算单位时间T内危险截面θ_k方向上引起的疲劳损伤，其表达式为

$D_{\mathrm{ik}}=\underset{r=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{16}{\mathrm{\&Sigma;}}}{D}_{\mathrm{irjk}}{P}_{\mathrm{rj}},k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,16$

取上述疲劳损伤的最大值，记为则高耸塔器在第i种危险截面部位结构形式下的风致疲劳寿命为

在实际工作中，还可以增加以下三步：

S6、根据S5的结果判定是否存在满足设计寿命的评定区域，且该区域应包含制造允许的容许缺陷要求；

S7、若步骤S6判定结果为“否”，则需重新设计高耸塔器整体结构，并返回步骤S2；

S8、若步骤S6判定结果为“是”，则确定结构设计方案。

实际工作中，在步骤S5中计算疲劳损伤D_irjk时，针对不同的焊接形式和焊接缺陷，英国标准BSI PD5500:2009给出了D、E、F、F2、G、W六种评定等级下的疲劳曲线；使用时根据高耸塔器具体的焊接形式与焊接缺陷确定评定等级，再由评定等级下的疲劳曲线计算单位时间T内上述应力循环对应的疲劳损伤D_irjk即可。

所述不同的焊接形式和焊接缺陷所对应的疲劳曲线是反应自然规律的人工绘制的曲线，前述按照英国标准BSI PD5500:2009给出的D、E、F、F2、G、W六种评定等级下的疲劳曲线进行计算是为了提高计算效率，而不是必须按照英国标准BSI PD5500:2009进行计算。

本发明的主要优点在于：

1）、本发明模拟了高耸塔器在危险截面各方向上的风振动态响应过程，这与标准JB/T4710-2005《钢制塔式容器》与GB50135-2006《高耸结构设计规范》只能计算单一的等效风振响应值相比，能够再现高耸塔器的风振响应随时间变化的过程；

2）、本发明计算疲劳损伤时可以单独考虑应力集中因子的影响，使得该方法能够仅在计算一次风振响应的条件下，就能分析危险截面不同结构对风致疲劳寿命带来的影响，这对基于寿命预测的危险截面局部优化设计带来了很大的便利；

3）、危险截面的应力循环是根据雨流法对应力幅和平均应力进行直接统计，这比其他计算疲劳损伤的间接等效计算方法更加精确，此外还可以提供应力幅和平均应力的循环周次分布情况；

4）、该发明考虑了风速风向对高耸塔器风致疲劳寿命的影响，可以根据高耸塔器所在地的气象条件选择适合的风速风向联合分布函数，这使得高耸塔器风致疲劳寿命分析更加科学准确；

5）、本发明可以根据不同的无损检测要求和容许缺陷等级对其结构疲劳寿命进行校核，设计结果的预知性保证了塔器的安全运行。

附图说明

图1为本发明的操作流程图。

图2为风速与风向角示意图。

图3a、图3b为实施例中两个风向的离散风速分布情况。

图4a为实施例高耸塔器的危险截面示意图。

图4b是图4a中的危险截面部位局部放大图。

图5a为实施例中的结构设计形式示意图。

图5b为当焊脚高度ht＝45mm时的危险截面部位的结构示意图。

图5c为当焊脚高度ht＝25mm时的危险截面部位的结构示意图。

图6a、6b均为实施例中焊脚高度为ht=45mm时危险截面部位的实体有限元模型。

图6c、6d为图6a、6b中危险截面部位的局部放大图。

图7为实施例塔节的划分与模拟点示意图。

图8为实施例中的顺风向风载荷时程样本。

图9为实施例中的横风向风载荷时程样本。

图10为实施例中顺风向名义应力时程响应。

图11为实施例中横风向名义应力时程响应。

图12a、12b为实施例中参考风速为28m/s时两个夹角上的名义应力时程响应。

图13为实施例中危险截面某一方向的应力循环统计情况。

具体实施方式

本实施例中的具体计算部分没有给出，在实际计算时请参考现有技术中的计算规则。

本发明步骤S5中的疲劳曲线、焊接形式与焊接缺陷评定等级的确定，可参考英国标准BSI5500-2009《Specification for Unfired fusion welded pressure vessels》与BS7910-2005《Guidance on methods for assessing the acceptability of flaws inmetallic structures》；技术方案中的具体计算过程需通过MATLAB软件编程实现，涉及有限元法求解部分需利用ANSYS有限元软件实现。

以下结合附图和实施例对本方法做进一步说明。

对沿海的某台高耸塔器进行用于抗风防振的疲劳寿命设计，该高耸塔器高约52.5m，筒体内直径2.6m，封头壁厚14mm，上筒体壁厚12mm，下筒体壁厚14mm，材料为Q345R。

步骤S1，首先设定参考风速为V_r＝4，6，8，…，26，28m/s，相邻的参考风速间隔为2m/s；其次，步骤S1中的频度参数f_j、尺度参数a_j、形状参数r_j除选用表格中的默认参数值以外，还可以选用该地区的风速风向联合分布函数的统计数据，本实施例中，由于高耸塔器所在地区尚无风速风向联合分布函数的统计数据，故频度参数f_j、尺度参数a_j、形状参数r_j均选用默认参数值，并将其离散为与参考风速对应的风速分布区间。风向与风向角的对应关系如图2所示，其中NNE和E两个风向上的风速分布离散化结果如图3a、3b所示。

步骤S2，使用ANSYS有限元软件SHELL93壳单元确定危险截面位置，水平截面A-A为该高耸塔器的危险截面，参见图4a、4b所示。

步骤S3，首先设计危险截面部位不同的结构形式。标准JB/T4710-2005规定该处的焊脚高度ht应不小于1.7倍的裙座名义厚度δ_ns（参见图5a），即焊脚高度ht不小于23.8mm，为此设定25mm、30mm、35mm、40mm、45mm五种可行的焊脚高度值，图5b所示的为焊脚高度ht的上限值，图5c所示的为焊脚高度ht的下限值。

其次，利用ANSYS有限元软件SOLID95实体单元，如图6a、6b、6c、6d所示，建立高耸塔器危险截面部位的局部有限元模型，分别计算上述五种焊脚高度值所对应的结构的应力集中因子，其结果如表1所示。

表1有限元法求解应力集中因子

焊脚高度ht/mm	45	40	35	30	25
						应力集中因子Ki	3.393	3.858	4.391	5.009	5.751

步骤S4，分以下几步获得该高耸塔器各风向上的名义应力时程响应。

1）模拟风载荷

（1）将该塔器分为10段塔节，塔节的划分和模拟点的分布见图7，每段塔节的l_h、D_h、D_eh和z_h参数如表2所示。

表2塔节的相关参数值（单位：m）

h	1	2	3	4	5	6	7	8	9	10
											lh	5.65	7.35	5.00	5.00	5.00	5.00	5.00	5.00	5.00	4.55
Dh	3.50	2.70	2.70	2.70	2.70	2.70	2.70	2.70	2.70	2.70
											Deh	5.13	4.33	4.33	4.33	4.33	4.33	4.33	4.33	4.33	4.33
zh	2.83	9.33	15.5	20.5	25.5	30.5	35.5	40.5	45.5	50.3

（2）计算各参考风速下各模拟点的平均风速，以参考风速V_r=16m/s为例，各模拟点的平均风速V_rh如表3所示。

表3各模拟点的平均风速（单位：m）

h	1	2	3	4	5	6	7	8	9	10
											Vrh	13.1	15.8	17.2	17.9	18.6	19.1	19.6	20.0	20.4	20.7

（3）利用MATLAB编程，模拟了该高耸塔器在各参考风速下顺风向风载荷时程样本。谐波叠加法模拟参数如下：起始频率为ω_down=0.001rad/s，截止频率为ω_up=8rad/s,频率等分数为N=2048，单位时间为T=600s，采样时距为t=0.1s。以上参数值适用于各参考风速，图8示出了参考风速V_r=20m/s时，第4模拟点前200s的顺风向风载荷时程样本。

（4）利用MATLAB编程，模拟了该高耸塔器在各参考风速下横风向风载荷时程样本。谐波叠加法的模拟参数如下：起始频率为ω_down=0.001rad/s，截止频率为ω_up=12rad/s,频率等分数为N=2048，单位时间为T=600s，采样时距为t=0.1s。以上参数值适用于各参考风速，但不同参考风速还需给出特定的参数值，如摩擦速度u_*r、旋涡脱落圆频率ω_srh和阻力系数C_rh。以参考风速V_r=24m/s为例，计算所得的摩擦速度为u_*r＝1.536m/s，参数ω_srh和C_rh的值如表4所示，模拟获得的第4点的横风向风载荷时程样本如图9所示。

表4各模拟点的ω_srh、C_rh参数值

2）模拟名义应力时程响应

首先，求解各参考风速下的顺风向风振响应。利用ANSYS软件建立该高耸塔器的梁单元有限元模型，阻尼比设为0.01，分别加载各参考风速下顺风向风载荷时程样本，对模型进行FULL瞬态动力学求解，并提取顺风向名义应力时程响应图10示出了参考风速为28m/s时高耸塔器危险截面的顺风向名义应力时程响应。

其次，求解各参考风速下的横风向风振响应。各参考风速下的高耸塔器横风向风振响应阻尼比的计算结果如表5所示，据此设定各参考风速下的阻尼比，分别加载各参考风速下横风向风载荷时程样本，对模型进行FULL瞬态动力学求解，并提取横风向名义应力时程响应图11示出了参考风速为28 m/s时高耸塔器危险截面的横风向名义应力时程响应。

表5 各参考风速下的横风向阻尼比

3）叠加名义应力时程响应

根据公式 ${\mathrm{\&sigma;}}_{\mathrm{rjk}}\left(t\right)={\mathrm{\&sigma;}}_r^a\left(t\right)\cos ({\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_k)+{\mathrm{\&sigma;}}_r^c\left(t\right)\sin ({\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_k),$ 将各参考风速下的顺风向名义应力时程响应与横风向名义应力时程响应，叠加为不同风向角和危险截面方向上的名义应力时程。图12a、12b示出了参考风速为28 m/s时两种夹角下危险截面的名义应力时程响应。

步骤S5，针对不同的结构设计形式，分别利用雨流法统计了各参考风速下危险截面各方向上的应力循环。参考风速为28 m/s、焊脚高度为45 mm时，危险截面某一夹角下T=600s内的应力循环统计情况如图13所示。最后，计算该高耸塔器在不同结构形式和缺陷评定等级下的风致疲劳寿命，计算结果如表6所示。

表6不同焊脚高度及评定等级下的疲劳寿命（单位：年）

步骤S6，判定是否满足设计条件。参考标准BSI 5500-2009中缺陷等级的评定部分，其规定危险截面处的单面焊焊接接头在无明显缺陷时应定为E级，而F、F2、G、W均代表不同等级的缺陷，W代表最严重的缺陷等级。已知该塔器的设计寿命为20年，则表3中灰色背景的区域均可以满足设计寿命要求，且部分结构容许F、F2级缺陷。故判定结果为“是”，即经初次结构设计后不必执行步骤S7与返回步骤S2。

步骤S8，确定结构设计方案。当焊脚高度ht为30、35mm时，其疲劳寿命在无明显缺陷的E等级下可以满足要求，但是这两种结构形式没有容许的缺陷等级，若该塔器在制造时产生了缺陷或者无损检测时未能发现存在的缺陷，其使用寿命可能达不到设计要求；当焊脚高度分别ht为40、45mm时，其疲劳寿命可以满足要求，并且该塔器在缺陷等级分别为F、F2时，其疲劳寿命仍可以满足设计要求。因此，应将危险截面处的焊接结构设计为焊脚高度ht在40mm以上，使之其风致疲劳寿命不仅可以满足设计要求，而且还有一定的缺陷容许等级，从而保证该塔器的安全运行。

Claims (1)

1.一种用于防风抗振的高耸塔器疲劳寿命设计方法，其包括如下步骤：

S1、以离地面10m高度处的平均风速作为参考风速，设定参考风速为V_r，r＝1，2，…，13；最小参考风速V₁取4m/s，最大参考风速V₁₃取28m/s，中间相邻的参考风速取值间隔为2m/s；

根据高耸塔器所在的地区，选择风速风向联合分布函数，所述风速风向联合分布函数的表达式如下：

$P(u<U_r,{\mathrm{\&theta;}}_j)=f_j\&times;\{1-\exp \&lsqb;-{(U_r/a_j)}^{r_j}\&rsqb;\},j=1,2,......,16$

式中，u为风速风向联合分布函数的风速自变量，单位为m/s；

U_r为风速风向联合分布函数的风速值，单位为m/s；

θ_j为0°到360°之间平均分布的16个风向角，且按顺时针依次增大，

θ_j＝(j－1)×22.5°，各个风向角θ_j代表不同的风向：如θ₁＝0°代表正北(N)风向，θ₅＝90°代表正东(E)方向，θ₉＝180°代表正南(S)风向，θ₁₁＝225°代表西南(SW)风向；

f_j代表风向角θ_j的频度参数，a_j、r_j分别是风向角θ_j上风速分布的尺度参数、形状参数，f_j、a_j、r_j采用下表中的数据：

风速风向联合分布函数f_j、a_j、r_j的默认参数值

j 风向 f_j a_j r_j j 风向 f_j a_j r_j 1 N 0.123 4.724 1.592 9 S 0.050 4.374 1.267 2 NNE 0.084 3.222 1.160 10 SSW 0.050 5.140 1.927 3 NE 0.063 2.630 1.073 11 SW 0.035 4.215 1.448 4 ENE 0.031 2.057 0.822 12 WSW 0.016 3.868 1.377 5 E 0.076 3.975 1.033 13 W 0.006 3.810 1.315 6 ESE 0.168 3.154 1.009 14 WNW 0.003 3.622 1.596 7 SE 0.147 2.351 0.908 15 NW 0.023 3.597 1.133 8 SSE 0.055 3.923 1.380 16 NNW 0.069 4.347 1.593

将风速风向联合分布函数离散成与各参考风速V_r相关的多个连续区间，使得风速风向联合分布函数的每个区间[U_r，U_r+1)的中值为V_r，r＝1，2，…，13；参考风速V_r在第j个风向角θ_j上的概率P_rj为

P_rj＝P(u＜U_r+1,θ_j)-P(u＜U_r,θ_j)

S2、利用ANSYS软件，建立高耸塔器壳单元有限元模型，在塔顶施加弯矩进行静力学求解，高耸塔器轴向应力最大值所在的水平截面为高耸塔器的危险截面；

S3、按照无损检测及容许缺陷要求，设计若干种达到无损检测及容许缺陷要求的危险截面部位结构形式，并分别建立危险截面部位的局部实体单元ANSYS有限元模型，求解危险截面部位在上述若干种结构形式下的应力集中因子；

S4、模拟单位时间内高耸塔器在各参考风速V_r下的危险截面部位的不同方向上的名义应力时程响应，包括如下步骤：

1)、模拟风载荷

(1)将高耸塔器沿竖直方向划分为n段塔节，划分标准参考高耸塔器的设计文件，确定每段塔节的长度l_h、外直径D_h、等效直径D_eh以及每段塔节的中心点距地面高度z_h，其中h＝1，2，3，…，n；各段塔节的中心点即风载荷的模拟点，l_h、D_h、D_eh和z_h的单位为m；

(2)计算各参考风速V_r下各模拟点的平均风速V_rh，其表达式如下：

$V_{rh}=V_r{\left(\frac{z_h}{10}\right)}^{0.16},h=1,2,3,...,n$

(3)计算各参考风速V_r下的顺风向风载荷时程样本

首先计算参考风速V_r下的顺风向脉动风速谱矩阵 中的上标a代指顺风向，其表达式如下

式中，矩阵元素

n为高耸塔器的塔节数；

X_r为脉动风速谱的无量纲莫宁坐标，

ω为圆频率，单位为rad/s；

为第h₁、h₂段塔节的中心点距地面高度，即第h₁、h₂个模拟点的高度，单位为m；

为参考风速V_r下第h₁、h₂个模拟点的平均风速，单位为m/s；

则各参考风速V_r下的顺风向脉动风速时程样本由谐波叠加法获得：

式中，h＝1,2,…,n；

Δω＝(ω_up－ω_down)/N；

N为频率等分数，N≥1024；

ω_up为截止频率，ω_down为起始频率，单位为rad/s；

为均匀分布于0～2π之间的随机相位角，单位为rad/s；

ω_l＝ω_down+(l－1/2)Δω，l＝1，2，…，N；

为的下三角矩阵第h行、m列元素，单位为m·s^1/2/rad^1/2；

时间t的增量Δt必须满足：Δt≤π/ω_up，t＝Δt，2Δt，3Δt，…，T，单位为s；

T为单位时间，T为增量Δt的整数倍，即每次模拟样本的采样时间，单位为s；

于是，各参考风速V_r下的顺风向风载荷时程样本可写成如下形式：

$P_{rh}^a\left(t\right)=0.5{\mathrm{\&mu;}}_s{\mathrm{\&rho;D}}_{eh}{l}_h{\&lsqb;V_{rh}+v_{rh}^a\left(t\right)\&rsqb;}^2,h=1,2,...,n$

式中，为高耸塔器第h段塔节的顺风向风载荷时程样本，单位为牛顿(N)；

μ_s为体型系数，μ_s＝0.7；

ρ为空气密度，取ρ＝1.25Kg/m³；

l_h为高耸塔器第h段塔节的长度，单位为m；

D_eh为高耸塔器第h段塔节的等效直径，单位为m；

(4)计算各参考风速V_r下的横风向风载荷时程样本

首先计算参考风速V_r下的横风向风载荷谱矩阵 中的c代指横风向，其表达式如下：

式中，矩阵元素

n为高耸塔器的塔节数；

ω为圆频率，单位为rad/s；

为第h₁、h₂个模拟点的高度，单位为m；

为参考风速V_r下的旋涡脱落力互谱密度，单位为N²·s/rad；

为参考风速V_r下的来流紊流力互谱密度，单位为N²·s/rad；

由以下二式获得：

$S_r^{c1}(z_{h_1},z_{h_2},\mathrm{\&omega;})=\sqrt{S_r^{c1}(z_{h_1},\mathrm{\&omega;})S_r^{c1}(z_{h_2},\mathrm{\&omega;})}\cos \left(\frac{1}{2}R\right)\exp (-\frac{1}{3}{R}^2),h_1,h_2=1,2,3,...,n;$

$S_r^{c1}(z_h,\mathrm{\&omega;})={\&lsqb;\frac{1}{2}{\mathrm{\&rho;D}}_h{V}_{rh}^2\&rsqb;}^2\frac{{\mathrm{\&sigma;}}_{CL}^2}{2\mathrm{\&pi;}\sqrt{\mathrm{\&pi;}}{\mathrm{B\&omega;}}_{sh}}\exp \{-{\&lsqb;(1-\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_{sh}})/B\&rsqb;}^2\},h=1,2,3,...,n;$

式中，为参考风速V_r下的旋涡脱落力自谱密度，单位为N²·s/rad；

R为距离参数，

ρ为空气密度，取ρ＝1.25Kg/m³；

σ_CL为升力系数的方差，取0.14；

B是谱曲线带宽经验参数，取0.18；

ω_sh为高耸塔器h段塔节的旋涡脱落圆频率，ω_sh＝V_rh/(10πD_h)，单位为rad/s；

分别为第h₁、h₂段塔节的外直径，单位为m；

由以下两式获得：

$S_r^{c1}(z_{h_1},z_{h_2},\mathrm{\&omega;})=\sqrt{S_r^{c1}(z_{h_1},\mathrm{\&omega;})S_r^{c1}(z_{h_2},\mathrm{\&omega;})}\exp (-18.7\frac{\mathrm{\&omega;}\mathrm{\&pi;}|z_{h_1}-z_{h_2}|}{V_{{\mathrm{rh}}_1}+V_{{\mathrm{rh}}_2}}),h_1,h_2=1,2,3,...,n;$

$S_r^{c2}(z_h,\mathrm{\&omega;})={\left(\frac{1}{2}{\mathrm{\&rho;C}}_{rh}{D}_h{V}_{rh}\right)}^2\frac{15u_{*r}^2{X}_{rh}^c}{2\mathrm{\&pi;}\mathrm{\&omega;}{(1+9.5X_{rh}^c)}^{5/3}},h=1,2,3,...,n;$

式中，为来流紊流谱的无量纲莫宁坐标,

为来流紊流力自谱密度，单位为N²·s/rad；

C_rh为阻力系数,中间值按线性插值计算；

Re_rh为雷诺数，Re_rh＝69000V_rhD_h；

u_*r为摩擦速度，u_*r＝0.064V_r，单位为m/s；

则各参考风速V_r下的横风向单位长度的风载荷时程样本由谐波叠加法公式获得：

式中，h＝1,2,…,n；

Δω＝(ω_up－ω_down)/N；

N为频率等分数，N≥1024；

ω_up为截止频率，ω_down为起始频率，单位为rad/s；

为均匀分布于0～2π之间的随机相位角，单位为rad/s；

ω_l＝ω_down+(l－1/2)Δω，l＝1，2，…，N；

为的下三角矩阵第h行、m列元素，单位为m·s^1/2/rad^1/2；

T为单位时间，即每次模拟样本的采样时间，单位为s；

时间t的增量Δt必须满足：Δt≤π/ω_up，t＝Δt，2Δt，3Δt，…，T，T为增量Δt的整数倍，单位为s；

于是，各参考风速V_r下各点横风向风载荷时程样本为：

$P_{rh}^c\left(t\right)=l_h{f}_{rh}^c\left(t\right),h=1,2,3,...,n;$

式中，为高耸塔器第h段塔节的横风向风载荷时程样本，单位为牛顿(N)；

l_h为高耸塔器第h段塔节长度，单位为m；

2)、模拟名义应力时程响应

利用ANSYS有限元软件建立高耸塔器梁单元有限元模型，分别加载各参考风速V_r下的顺风向风载荷时程样和横风向风载荷时程样本进行FULL瞬态动力学模块求解；利用ANSYS软件APDL编程，提取各参考风速下危险截面的顺风向名义应力时程响应和横风向名义应力时程响应

求解高耸塔器顺风向风振响应时，高耸塔器的阻尼比为ζ^a＝0.01；

求解高耸塔器横风向风振响应时，参考风速V_r下高耸塔器的阻尼比为

式中，为求解横风向风振响应时高耸塔器的阻尼比；

M^*为高耸塔器的广义质量，

M_h为高耸塔器第h段塔节的质量，由设计文件给出，单位为Kg；

l_h为高耸塔器第h段塔节的长度，单位为m；

为高耸塔器第h段塔节振型函数值，计算公式如下：

H为高耸塔器总高度，单位为m；

A_rh为高耸塔器第h段塔节的气动力系数；

$A_{rh}=\left\{\begin{array}{cc}0& V_{rh}/V_{rh}^{cr}<0.85\\ a_h(3.5V_{rh}/V_{rh}^{cr}-2.95)& 0.85\&le;V_{rh}/V_{rh}^{cr}<1.0\\ 0.55a_h& 1.0\&le;V_{rh}/V_{rh}^{cr}<1.1\\ a_h(2.75-2V_{rh}/V_{rh}^{cr})& 1.1\&le;V_{rh}/V_{rh}^{cr}<1.3\\ a_h(0.46-0.25V_{rh}/V_{rh}^{cr})& 1.3\&le;V_{rh}/V_{rh}^{cr}<1.84\\ 0& 1.84\&le;V_{rh}/V_{rh}^{cr}\end{array}\right.$

为第h段塔节的临界风速，单位为m/s；

n_s为高耸塔器的一阶自振频率，可由ANSYS软件模态求解得出；

a_h为第h段塔节的气动力参数，a_h＝a_1ha_2ha_3h，h＝1，2，3，…，n；

a_1h、a_2h、a_3h分别由以下三式确定：

$a_{1h}=\left\{\begin{array}{cc}1.0& {\mathrm{Re}}_{rh}<{10}^4\\ 1.8& {10}^4\&le;{\mathrm{Re}}_{rh}<{10}^5\\ 1.0& {10}^5\&le;{\mathrm{Re}}_{rh}\end{array}\right.$

$a_{2h}=\left\{\begin{array}{cc}1.8& V_r\&le;12m/s\\ 0.9& V_r>12m/s\end{array}\right.$

$a_{3h}=\left\{\begin{array}{cc}1.0& z_h/D_h>12.5\\ 1.0-0.04\&lsqb;12.5-z_h/D_h\&rsqb;& z_h/D_h<12.5\end{array}\right.$

Re_rh为雷诺数，Re_rh＝69000V_rhD_h；

3)、叠加名义应力时程响应

当参考风速为V_r、风向角为θ_j时，高耸塔器的危险截面在方向θ_k(风向角的角度同时也是危险截面上的方向位置)上单位时间T内的名义应力时程响应σ_rjk(t)为

${\mathrm{\&sigma;}}_{rjk}\left(t\right)={\mathrm{\&sigma;}}_r^a\left(t\right)cos({\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_k)+{\mathrm{\&sigma;}}_r^c\left(t\right)sin({\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_k)$

式中，为顺风向名义应力时程响应，单位为MPa；

为横风向名义应力时程响应，单位为MPa；

S5、计算高耸塔器在步骤S3中不同结构形式和缺陷评定等级下的风致疲劳寿命，包括如下步骤：

1)、统计应力循环

(1)将步骤S₄中的叠加后的名义应力时程响应σ_rjk(t)修正为

σ_irjk(t)＝K_i·σ_rjk(t)，i＝1,2,…,p

式中，K_i为步骤S3中所求的第i种危险截面部位结构形式的应力集中因子；

p为S3中设计的危险截面部位结构形式个数；

σ_irjk(t)为修正后的应力时程响应，单位为MPa；

(2)利用雨流法统计σ_irjk(t)对应的应力循环；

2)、计算疲劳损伤

根据不同的焊接形式和焊接缺陷所对应的疲劳曲线计算单位时间T内上述应力循环对应的疲劳损伤D_irjk；

3)、计算疲劳寿命

结合离散后的风速风向联合分布，计算单位时间T内危险截面θ_k方向上引起的疲劳损伤，其表达式为

$D_{ik}=\underset{r=1}{\overset{N}{\&Sigma;}}\underset{j=1}{\overset{16}{\&Sigma;}}{D}_{irjk}{P}_{rj},k=1,2,...,16$

取上述疲劳损伤的最大值，记为则高耸塔器在第i种危险截面部位结构形式下的风致疲劳寿命为