当前位置:高等教育资讯网  >  中国高校课件下载中心  >  大学文库  >  浏览文档

全风向来流非高斯风场风机疲劳寿命可靠性分析

资源类别:文库,文档格式:PDF,文档页数:10,文件大小:10.51MB,团购合买
在Hermite矩模型基础上,根据Kaimal谱生成某典型风机结构正常风速条件下,三种不同概率特性风场(高斯、非高斯硬化和软化),在考虑来流风向和平均风速联合概率密度条件下,以塔架基础连接处为例,对风机进行疲劳寿命可靠性分析.由叶片的气动模型和多体动力,计算出风机的动力响应,并对响应的时域和频域特性进行分析.基于线性损伤累积理论和Paris公式,对来流全风向条件下的裂纹形成寿命和裂纹扩展寿命进行了详细讨论.结果表明,裂纹形成寿命对风荷载的非高斯性较为敏感,而裂纹扩展寿命对风荷载的非高斯性并不敏感,需要考虑风荷载的非高斯性对风机结构疲劳损伤的影响.此外,在考虑全风向来流条件下,疲劳裂纹形成和扩展阶段的失效位置相同,均在主导风向上.
点击下载完整版文档(PDF)

工程科学学报,第39卷,第9期:1453-1462,2017年9月 Chinese Journal of Engineering,Vol.39,No.9:1453-1462,September 2017 D0l:10.13374/j.issn2095-9389.2017.09.020:htp:/journals..usth.edu.cn 全风向来流非高斯风场风机疲劳寿命可靠性分析 双 妙12),宋波12) 1)北京科技大学土木与资源工程学院,北京1000832)北京科技大学强震区轨道交通工程抗震研究北京市国际科技合作基地,北京 100083 ☒通信作者,E-mail:b20140020@x.usth.edu.cm 摘要在Hermite矩模型基础上,根据Kaimal谱生成某典型风机结构正常风速条件下,三种不同概率特性风场(高斯、非高 斯硬化和软化),在考虑来流风向和平均风速联合概率密度条件下,以塔架基础连接处为例,对风机进行疲劳寿命可靠性分 析.由叶片的气动模型和多体动力,计算出风机的动力响应,并对响应的时域和频域特性进行分析.基于线性损伤累积理论 和Pis公式,对来流全风向条件下的裂纹形成寿命和裂纹扩展寿命进行了详细讨论.结果表明,裂纹形成寿命对风荷载的非 高斯性较为敏感,而裂纹扩展寿命对风荷载的非高斯性并不敏感,需要考虑风荷载的非高斯性对风机结构疲劳损伤的影响 此外,在考虑全风向来流条件下,疲劳裂纹形成和扩展阶段的失效位置相同,均在主导风向上 关键词全风向;非高斯性:雨流计数法:裂纹形成;裂纹扩展 分类号TU311.3 Reliability analysis of the fatigue life of wind turbines under a non-Gaussian wind field with a full-direction inflow SHUANG Miao)SONG Bo'2) 1)School of Civil and Resource Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)Beijing International Cooperation Base for Science and Technology-Aseismic Research of the Rail Transit Engineering in the Strong Motion Area,Uni- versity of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:b20140020@xs.ustb.edu.en ABSTRACT Using the Hermite moment model,three types of wind fields with different probability characteristics,namely the Gaussian,hardening non-Gaussian,and softening non-Gaussian processes,were generated via the Kaimal spectrum for a typical wind turbine under operational conditions.Reliability analysis of the fatigue life of the wind turbine was performed by taking tower base con- nections as an example with consideration of the joint probability density distribution of the wind direction and mean wind speed.The dynamic response was calculated by an aerodynamic model of the blade and multi-body dynamics,and the time-and frequency-domain characteristics of the response were analyzed.Using linear damage accumulation theory and the Paris equation,the fatigue crack initia- tion life and crack growth life of the turbine were discussed in detail.Fatigue estimation shows that the crack initiation life of the tur- bine is more sensitive to non-Gaussian winds,whereas its crack propagation life is less sensitive to the non-Gaussian characteristics of wind load.The influence of the non-Gaussian characteristics of wind load on the fatigue damage of the wind turbine should be consid- ered.In terms of the full-direction inflow,the failure positions of the crack initiation and propagation stages are identical and in the dominant wind direction. KEY WORDS full-direction inflow;non-Gaussian characteristic;rain-flow counting method;crack initiation;crack propagation 收稿日期:2017-01-11 基金项目:国家自然科学基金资助项目(51178045):北京科技大学与台北科技大学专题联合研究计划资助项目(TW201601):高校2016年度 引智项目资助项目(110000201420160126)

工程科学学报,第 39 卷,第 9 期:1453鄄鄄1462,2017 年 9 月 Chinese Journal of Engineering, Vol. 39, No. 9: 1453鄄鄄1462, September 2017 DOI: 10. 13374 / j. issn2095鄄鄄9389. 2017. 09. 020; http: / / journals. ustb. edu. cn 全风向来流非高斯风场风机疲劳寿命可靠性分析 双 妙1,2) 苣 , 宋 波1,2) 1) 北京科技大学土木与资源工程学院, 北京 100083 2) 北京科技大学强震区轨道交通工程抗震研究北京市国际科技合作基地, 北京 100083 苣通信作者, E鄄mail: b20140020@ xs. ustb. edu. cn 收稿日期: 2017鄄鄄01鄄鄄11 基金项目: 国家自然科学基金资助项目(51178045);北京科技大学与台北科技大学专题联合研究计划资助项目(TW201601);高校 2016 年度 引智项目资助项目(110000201420160126) 摘 要 在 Hermite 矩模型基础上,根据 Kaimal 谱生成某典型风机结构正常风速条件下,三种不同概率特性风场(高斯、非高 斯硬化和软化),在考虑来流风向和平均风速联合概率密度条件下,以塔架基础连接处为例,对风机进行疲劳寿命可靠性分 析. 由叶片的气动模型和多体动力,计算出风机的动力响应,并对响应的时域和频域特性进行分析. 基于线性损伤累积理论 和 Paris 公式,对来流全风向条件下的裂纹形成寿命和裂纹扩展寿命进行了详细讨论. 结果表明,裂纹形成寿命对风荷载的非 高斯性较为敏感,而裂纹扩展寿命对风荷载的非高斯性并不敏感,需要考虑风荷载的非高斯性对风机结构疲劳损伤的影响. 此外,在考虑全风向来流条件下,疲劳裂纹形成和扩展阶段的失效位置相同,均在主导风向上. 关键词 全风向; 非高斯性; 雨流计数法; 裂纹形成; 裂纹扩展 分类号 TU311郾 3 Reliability analysis of the fatigue life of wind turbines under a non鄄Gaussian wind field with a full鄄direction inflow SHUANG Miao 1,2) 苣 , SONG Bo 1,2) 1) School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China 2) Beijing International Cooperation Base for Science and Technology鄄Aseismic Research of the Rail Transit Engineering in the Strong Motion Area, Uni鄄 versity of Science and Technology Beijing, Beijing 100083, China 苣Corresponding author, E鄄mail: b20140020@ xs. ustb. edu. cn ABSTRACT Using the Hermite moment model, three types of wind fields with different probability characteristics, namely the Gaussian, hardening non鄄Gaussian, and softening non鄄Gaussian processes, were generated via the Kaimal spectrum for a typical wind turbine under operational conditions. Reliability analysis of the fatigue life of the wind turbine was performed by taking tower base con鄄 nections as an example with consideration of the joint probability density distribution of the wind direction and mean wind speed. The dynamic response was calculated by an aerodynamic model of the blade and multi鄄body dynamics, and the time鄄 and frequency鄄domain characteristics of the response were analyzed. Using linear damage accumulation theory and the Paris equation, the fatigue crack initia鄄 tion life and crack growth life of the turbine were discussed in detail. Fatigue estimation shows that the crack initiation life of the tur鄄 bine is more sensitive to non鄄Gaussian winds, whereas its crack propagation life is less sensitive to the non鄄Gaussian characteristics of wind load. The influence of the non鄄Gaussian characteristics of wind load on the fatigue damage of the wind turbine should be consid鄄 ered. In terms of the full鄄direction inflow, the failure positions of the crack initiation and propagation stages are identical and in the dominant wind direction. KEY WORDS full鄄direction inflow; non鄄Gaussian characteristic; rain鄄flow counting method; crack initiation; crack propagation

·1454· 工程科学学报,第39卷,第9期 风机结构自振频率与高度在60~160m之间的高 轮教高度90m叶片直径126m的5MW三叶片陆上风 层建筑相似,并且风机在设计中不考虑舒适度,因此风 机为例)].该风机采用变速变桨控制系统,额定风速 振响应比一般结构更为显著.同时,由于叶片与来流 为11.4ms1,额定转速12.1rmin.对结构在平稳 之间相互作用且机舱结构接近于刚体质点,动力响应 高斯、硬化和软化非高斯3种不同概率特性风荷载作 按传统有限元计算会产生较大误差.常用方法可分为 用下的动力响应和疲劳可靠性进行分析 多体动力-)和基于叶片与来流的气动模型两大类. 算例采用FAST(fatigue,aerodynamics,structures 目前的风机设计规范中,风场通常简化为平稳高 and turbulence)建模,风机叶片划分为I7个单元,塔筒 斯过程.对于平坦均匀地形的边界层风场,高斯假设 划分为10个单元(如图1),塔底与地面刚接.叶片和 是恰当的.然而,实测数据表明,复杂地形的脉动风往 塔体的材料参数和气动参数可参见文献[15]. 往显示出明显的非高斯特性.同时,非高斯荷载作用 10 20 下的建筑结构和海岸工程具有较大的极值响应和疲劳 损伤.因此需要根据风机结构的动力和控制特性,对 280 310 非高斯荷载作用下结构的广义动力响应和疲劳可靠性 截面2 进行分析.此外,由于来流风向随时间变化,疲劳估计 590 520 中需根据工程所处地区气象统计资料,考虑全风向来 流对疲劳损伤累积的影响. 对于风机等细长结构的风致疲劳问题,Do等[5-] 931 940 95 面1 964 根据Monte Carlo模拟和多体动力分析了塔架基础连 3.87m 09 接处的疲劳寿命,并提出基于性能的设计方法.Da- 270 wood等[)在疲劳损伤累积准则和雨流计数法基础上 5 mm /6=25mm 考虑材料的不确定性,对典型细长高耸结构的风致疲 截面1 截面2 180P 截面1 劳进行研究.Repetto与Solari[s-]根据动力响应的背 图1风机结构和风场示意图 景和共振两部分,采用双峰计数法[]研究了高斯荷载 Fig.1 Sketch of the wind turbine and wind field 作用下高耸结构的风致疲劳. 现有的风荷载非高斯性对风机结构影响的研究多 1.2平均风速和来流风向 集中于极值响应分析-],对疲劳损伤的研究相对较 风机结构任意高度处的风荷载均可分解为平均和 少.由破坏实例可知,风机失效多为叶片根部和塔架 脉动两部分.对于平均风荷载,可根据指数律由参考 基础连接处.目前,复合材料的广泛使用提高了叶 高度处平均风速换算得到.平均风速()的概率密度 片强度和疲劳寿命,因此本文将控制截面限定于塔架 函数可按Weibull分布简化,即: 基础连接处.结合结构和风荷载特性可知,塔体底部 ee (1) 截面由面内、面外弯曲应力和自重引起的压应力控制, 式中,c和k分别是尺度和形状参数,可根据工程所处 处于单轴应力状态.同时,控制截面不同位置的弯曲 地的风速统计资料按最大似然估计法(MLM)、统计矩 应力与来流风向相关,因此在可靠性分析中需考虑来 法(MM)和能量密度法(PDM)计算得到.考虑到方法 流风向的影响. 的简便性及结果的准确性,本文按能量密度法对两参 基于以上分析,本文在Hermite矩模型基础上,通 数Weibull分布进行估计,即uo]: 过非高斯穿越过程模拟了多维相关平稳高斯、非高斯 3.69 硬化和软化三种不同概率特性的风场,详细研究了三 k=1+ 种风场作用下风机结构广义动力响应的时域和频域特 (Ei (2) 性.结合气象资料,根据来流风向将控制截面划分成 不同控制点,采用线性疲劳损伤累积理论,对全风向来 r小+) 流不同概率特性风场作用下不同控制点的疲劳裂纹形 式中,T(·)是Gamma函数;E是能量系数,定义为: 成寿命进行讨论:并采用Pais公式,计算出疲劳裂纹 扩展寿命 E=7 (3) 1风机与来流特性 如图2所示为NREL135m气象塔88m高度处 超声风速仪测得的2014年全年10min平均风速(即轮 1.1风机结构 毂高度处的平均风速山)的玫瑰图.图3所示为3种 本文以NREL(national renewable energy laboratory) 不同方法Weibull分布参数估计与实测数据的比较

工程科学学报,第 39 卷,第 9 期 风机结构自振频率与高度在 60 ~ 160 m 之间的高 层建筑相似,并且风机在设计中不考虑舒适度,因此风 振响应比一般结构更为显著. 同时,由于叶片与来流 之间相互作用且机舱结构接近于刚体质点,动力响应 按传统有限元计算会产生较大误差. 常用方法可分为 多体动力[1鄄鄄2]和基于叶片与来流的气动模型[3鄄鄄4]两大类. 目前的风机设计规范中,风场通常简化为平稳高 斯过程. 对于平坦均匀地形的边界层风场,高斯假设 是恰当的. 然而,实测数据表明,复杂地形的脉动风往 往显示出明显的非高斯特性. 同时,非高斯荷载作用 下的建筑结构和海岸工程具有较大的极值响应和疲劳 损伤. 因此需要根据风机结构的动力和控制特性,对 非高斯荷载作用下结构的广义动力响应和疲劳可靠性 进行分析. 此外,由于来流风向随时间变化,疲劳估计 中需根据工程所处地区气象统计资料,考虑全风向来 流对疲劳损伤累积的影响. 对于风机等细长结构的风致疲劳问题,Do 等[5鄄鄄6] 根据 Monte Carlo 模拟和多体动力分析了塔架基础连 接处的疲劳寿命,并提出基于性能的设计方法. Da鄄 wood 等[7]在疲劳损伤累积准则和雨流计数法基础上 考虑材料的不确定性,对典型细长高耸结构的风致疲 劳进行研究. Repetto 与 Solari [8鄄鄄9] 根据动力响应的背 景和共振两部分,采用双峰计数法[10]研究了高斯荷载 作用下高耸结构的风致疲劳. 现有的风荷载非高斯性对风机结构影响的研究多 集中于极值响应分析[11鄄鄄13] ,对疲劳损伤的研究相对较 少. 由破坏实例可知,风机失效多为叶片根部和塔架 基础连接处[14] . 目前,复合材料的广泛使用提高了叶 片强度和疲劳寿命,因此本文将控制截面限定于塔架 基础连接处. 结合结构和风荷载特性可知,塔体底部 截面由面内、面外弯曲应力和自重引起的压应力控制, 处于单轴应力状态. 同时,控制截面不同位置的弯曲 应力与来流风向相关,因此在可靠性分析中需考虑来 流风向的影响. 基于以上分析,本文在 Hermite 矩模型基础上,通 过非高斯穿越过程模拟了多维相关平稳高斯、非高斯 硬化和软化三种不同概率特性的风场,详细研究了三 种风场作用下风机结构广义动力响应的时域和频域特 性. 结合气象资料,根据来流风向将控制截面划分成 不同控制点,采用线性疲劳损伤累积理论,对全风向来 流不同概率特性风场作用下不同控制点的疲劳裂纹形 成寿命进行讨论;并采用 Paris 公式,计算出疲劳裂纹 扩展寿命. 1 风机与来流特性 1郾 1 风机结构 本文以 NREL(national renewable energy laboratory) 轮毂高度 90 m 叶片直径 126 m 的 5 MW 三叶片陆上风 机为例[15] . 该风机采用变速变桨控制系统,额定风速 为 11郾 4 m·s - 1 ,额定转速 12郾 1 r·min - 1 . 对结构在平稳 高斯、硬化和软化非高斯 3 种不同概率特性风荷载作 用下的动力响应和疲劳可靠性进行分析. 算例采用 FAST ( fatigue, aerodynamics, structures and turbulence)建模,风机叶片划分为 17 个单元,塔筒 划分为 10 个单元(如图 1),塔底与地面刚接. 叶片和 塔体的材料参数和气动参数可参见文献[15]. 图 1 风机结构和风场示意图 Fig. 1 Sketch of the wind turbine and wind field 1郾 2 平均风速和来流风向 风机结构任意高度处的风荷载均可分解为平均和 脉动两部分. 对于平均风荷载,可根据指数律由参考 高度处平均风速换算得到. 平均风速( u)的概率密度 函数可按 Weibull 分布简化,即: f(u) = k ( c u ) c k - 1 ·e - (u / c) k . (1) 式中,c 和 k 分别是尺度和形状参数,可根据工程所处 地的风速统计资料按最大似然估计法(MLM)、统计矩 法(MM)和能量密度法(PDM)计算得到. 考虑到方法 的简便性及结果的准确性,本文按能量密度法对两参 数 Weibull 分布进行估计,即[16] : k = 1 + 3郾 69 (Epf) 2 ; c = u 祝 ( 1 + 1 ) k ì î í ï ïï ï ïï . (2) 式中,祝(·)是 Gamma 函数;Epf是能量系数,定义为: Epf = u 3 u 3 . (3) 如图 2 所示为 NREL 135 m 气象塔[17] 88 m 高度处 超声风速仪测得的 2014 年全年 10 min 平均风速(即轮 毂高度处的平均风速 uhub )的玫瑰图. 图 3 所示为 3 种 不同方法 Weibull 分布参数估计与实测数据的比较, ·1454·

双妙等:全风向来流非高斯风场风机疲劳寿命可靠性分析 ·1455· ☐25ms3)和硬化过程(y,<3).任意标准非高斯过程 湍流积分尺度,分别取L。=8.1A,L.=2.74.和L.= X(t)均可由单调函数变换成标准高斯过程U(t),即:

双 妙等: 全风向来流非高斯风场风机疲劳寿命可靠性分析 图 2 风速玫瑰图 Fig. 2 Wind rose diagram 图中柱状图是实测数据,曲线分别是按 PDM、MM 和 MLM 方法计算得到的两参数 Weibull 分布的概率密度 函数(PDF). 图 3 三种方法 Weibull 分布参数估计与实测数据的比较 Fig. 3 Comparison of parameters estimated by the Weibull distribu鄄 tion using three methods and the measured data 由图 3 可知,3 种方法的计算结果相近,但能量密 度法不需通过迭代即可快速估计出 Weibull 分布的形 状和尺度参数,且结果精度高,具有明显的优势. 1郾 3 脉动风速 根据脉动风速谱,采用 Monte Carlo 方法模拟如图 1 所示风场,风场中心位于轮毂位置. 本文脉动风速谱 采用 IEC 61400鄄鄄 1 ( International Electronic Technical Commission)规范中 Kaimal 谱模型(式(4)),但本文方 法并不限于 Kaimal 谱,对其他风速谱同样成立. Sk(f) = 4滓 2 kLk / uhub (1 + 6fLk / uhub ) 5 / 3 . (4) 式中:f 表示角频率;滓k是脉动风速的标准差( k 取 u、v 和 w,分别表示顺风向、横风向和垂直向),取 滓u = 0郾 16(0郾 75uhub + 5郾 6),滓v = 0郾 8滓u和 滓w = 0郾 5滓u ;Lk 是 湍流积分尺度,分别取 Lu = 8郾 1撰u ,Lv = 2郾 7撰u和 Lw = 0郾 66撰u ,湍流尺度系数 撰u = 0郾 7min(30, Hhub ),Hhub是 轮毂高度. 本文根据 IEC 规范,仅考虑空间中任意 i,j 两点顺 风向脉动风荷载的相关性,表示为: cohij,u (f) = e - au (f·r/ uhub) 2 + (bu r) 2 . (5) 式中,r 表示生成风速网格点间的距离;au 、bu表示衰减 和偏移系数,按 IEC 规范取值,即 au = 12 和 bu = 3郾 527 伊 10 - 4 . 风机正常风速条件下按 IEC 规范中的设计荷载条 件(DLC)1郾 1 和正常湍流模型(NTM)考虑载荷. 2 风场模拟 2郾 1 高斯风场模拟 根据谐波叠加法,生成平稳高斯脉动风速第 m 个 时程在 t 时刻的风速,表示为: uG,m (t) = 移 N/ 2 -1 j = -N/ 2 Fm (棕j)·e i棕j t . (6) 式中,Fm (棕j)表示第 m 个待生成脉动风速时程的第 j 个傅里叶谱,即: Fm (棕j) = 移 m r = 1 驻棕Hmr(棕j)·e i渍rj . (7) 式中:棕j = j驻棕,驻棕 是圆频率间隔;渍rj是第 r 个待生成 时程相位谱 渍r的第 j 个随机相位角;Hmr(棕j)是傅里叶 幅值谱,由脉动风速谱 Cholesky 分解得到,即: S(棕) = H(棕)H *T (棕). (8) 式中,S(棕) 是目标谱矩阵,H * ( 棕) 为 H( 棕) 的共轭 矩阵. 2郾 2 非高斯风场模拟 非高斯过程根据峰度系数( 酌4 ) 的不同分为软化 过程(酌4 >3)和硬化过程(酌4 < 3). 任意标准非高斯过程 X(t)均可由单调函数变换成标准高斯过程 U(t),即: ·1455·

·1456· 工程科学学报,第39卷,第9期 x=g(u)=F:[Φ(u)]. (9) 45m a 式中:u和x分别是U(t)和X()的随机变量;中(·)表 95%C 示标准高斯过程的概率分布函数(CDF);F表示标 30 准非高斯过程CDF的逆函数;g(·)表示非高斯穿越过 程,可根据Hermite矩模型us]或CDF映射等方法 95%CL 得到. 100 200 300 400 500 600 对于非高斯风荷载模拟,首先需对日标谱S(ω) 45 标准化,即: (b) 95%CI o)8 j,k=1,2,…,m. (10) 二 式中,Oc是第j个时程的标准差. 再根据式(6)生成标准高斯过程的风速时程 95%C ((t))后,通过非高斯转化过程(式(9))生成标准 100 200 300400 500 600 t切s 45 非高斯过程的风速时程(()): (c 95%CI+ u8(t)=g[u2(t)] (11) 为了提高非高斯过程的准确性,可对生成时程的 功率谱密度进行迭代,即20-2]: S(ω)= TS.c(o)1 Is2()s () (12) 100 200 300400500600 式中,B是收敛速度和收敛精度的优化系数,通常取 图4轮毂高度处风速时程.(a)硬化过程:(b)高斯过程:(c)软 B=0.3] 化过程 2.3非高斯风场正常风速模拟 Fig.4 Wind speed time history at the hub height:(a)hardening 在风机运行状态,将切入风速(3ms·)到切出风 process;(b)Gaussian process:(c)softening process 速(25ms1)划分为12个单元,单元宽度为2ms; 在风机停机状态,考虑低于切入风速的1ms和2m· 系统对结构动力响应有较大影响.因此本文采用基于 s风速,共14种不同平均风速的硬化、软化和高斯过 耦合气动模型的FAST程序计算结构的动力响应.为 程风场.本文按Hermite矩模型生成非高斯风场,软化 考虑来流风向影响,按图2所示将来流划分为16个风 过程的偏斜和峰度分别取0和4.5,硬化过程取0和 向角. 2.0.每种情况按相同的随机种子模拟硬化、软化和高 风机结构动力响应分析中,采用叶素动量理论考 斯风场各50条,每条时程的长度为630s,时间间隔为 虑叶片尾流模型,按Beddoes-Leishman模型考虑叶片 0.05$.为了消除启动瞬间的影响,从30s后对结构动 的气动效应.由FAST程序计算出塔体底部控制截面 力响应进行分析 的内力后,根据材料属性和几何尺寸计算出控制截面 图4所示为轮毂高度处平均风速25m·s时,按 不同位置的应力.由于面内弯曲应力远大于面外弯曲 式(6)和式(11)生成的硬化、高斯和软化风荷载.图 应力,因此本文仅对控制截面面内弯曲正应力(S,)的 中虚线表示均值加减1.96倍标准差(对于高斯分布即 动力响应进行分析 95%置信度区间C).图中小图比较了目标谱与生成 图5所示为风机结构在三种不同概率特性风荷载 荷载的模拟谱,右图实线为与风荷载均值和方差相同 作用下,控制截面面内最大弯曲正应力的响应时程 的高斯分布概率密度函数(PDF),柱状图为模拟风荷 可以看出,风机结构应力响应的非高斯性较风荷载的 载PDF.由图可知,三种不同风场功率谱密度与目标 非高斯性大大减弱. 谱均吻合较好,软化过程和硬化过程均表现出明显的 图6所示为不同风速区间内不同概率特性风场各 非高斯性.在相同目标谱、平均风速和湍流度条件下, 50条计算得到的应力响应的均值、标准差、偏态和峰 高斯过程的极值大小和数量介于硬化和软化过程 度的平均值.从图中可知在高斯和非高斯风荷载作用 之间. 下,均值、标准差和偏态系数基本相同.然而,由于高 阶矩的不稳定性,3种荷载作用下动力响应的峰度系 3风机结构的动力响应分析 数有所不同,但整体趋势仍保持一致 在风荷载作用下,叶片与来流、叶片与塔体间均相 图7所示为轮毂高度处平均风速25m·s条件 互耦合,动力响应分析具有高度的非线性.同时,控制 下,3种不同概率特性风荷载作用下,控制截面弯曲正

工程科学学报,第 39 卷,第 9 期 x = g(u) = F - 1 X [椎(u)]. (9) 式中:u 和 x 分别是 U(t)和 X(t)的随机变量;椎(·)表 示标准高斯过程的概率分布函数(CDF);F - 1 X 表示标 准非高斯过程 CDF 的逆函数;g(·)表示非高斯穿越过 程,可根据 Hermite 矩模型[18] 或 CDF 映射[19] 等方法 得到. 对于非高斯风荷载模拟,首先需对目标谱 S T NG (棕) 标准化,即: S 0 NG,jk(棕) = ST,NG,jk(棕) 滓NG,j滓NG,k ,j,k = 1,2,…,m. (10) 式中,滓NG,j是第 j 个时程的标准差. 再根据式 ( 6 ) 生 成 标 准 高 斯 过 程 的 风 速 时 程 (u 0 G (t))后,通过非高斯转化过程(式(9)) 生成标准 非高斯过程的风速时程(u 0 NG (t)): u 0 NG (t) = g[u 0 G (t)]. (11) 为了提高非高斯过程的准确性,可对生成时程的 功率谱密度进行迭代,即[20鄄鄄21] : S (i + 1) G (棕) = [ ST,NG (棕) S (i) NG (棕 ] ) 茁 S (i) G (棕). (12) 式中,茁 是收敛速度和收敛精度的优化系数,通常取 茁 = 0郾 3 [22] . 2郾 3 非高斯风场正常风速模拟 在风机运行状态,将切入风速(3 m·s - 1 )到切出风 速(25 m·s - 1 )划分为 12 个单元,单元宽度为 2 m·s - 1 ; 在风机停机状态,考虑低于切入风速的 1 m·s - 1和 2 m· s - 1风速,共 14 种不同平均风速的硬化、软化和高斯过 程风场. 本文按 Hermite 矩模型生成非高斯风场,软化 过程的偏斜和峰度分别取 0 和 4郾 5,硬化过程取 0 和 2郾 0. 每种情况按相同的随机种子模拟硬化、软化和高 斯风场各 50 条,每条时程的长度为 630 s,时间间隔为 0郾 05 s. 为了消除启动瞬间的影响,从 30 s 后对结构动 力响应进行分析. 图 4 所示为轮毂高度处平均风速 25 m·s - 1时,按 式(6)和式(11)生成的硬化、高斯和软化风荷载. 图 中虚线表示均值加减 1郾 96 倍标准差(对于高斯分布即 95% 置信度区间 CI ± ). 图中小图比较了目标谱与生成 荷载的模拟谱,右图实线为与风荷载均值和方差相同 的高斯分布概率密度函数( PDF),柱状图为模拟风荷 载 PDF. 由图可知,三种不同风场功率谱密度与目标 谱均吻合较好,软化过程和硬化过程均表现出明显的 非高斯性. 在相同目标谱、平均风速和湍流度条件下, 高斯过程的极值大小和数量介于硬化和软化过程 之间. 3 风机结构的动力响应分析 在风荷载作用下,叶片与来流、叶片与塔体间均相 互耦合,动力响应分析具有高度的非线性. 同时,控制 图 4 轮毂高度处风速时程. (a)硬化过程;( b)高斯过程;( c)软 化过程 Fig. 4 Wind speed time history at the hub height: ( a) hardening process; (b) Gaussian process; (c) softening process 系统对结构动力响应有较大影响. 因此本文采用基于 耦合气动模型的 FAST 程序计算结构的动力响应. 为 考虑来流风向影响,按图 2 所示将来流划分为 16 个风 向角. 风机结构动力响应分析中,采用叶素动量理论考 虑叶片尾流模型,按 Beddoes鄄鄄 Leishman 模型考虑叶片 的气动效应. 由 FAST 程序计算出塔体底部控制截面 的内力后,根据材料属性和几何尺寸计算出控制截面 不同位置的应力. 由于面内弯曲应力远大于面外弯曲 应力,因此本文仅对控制截面面内弯曲正应力( Sb )的 动力响应进行分析. 图 5 所示为风机结构在三种不同概率特性风荷载 作用下,控制截面面内最大弯曲正应力的响应时程. 可以看出,风机结构应力响应的非高斯性较风荷载的 非高斯性大大减弱. 图 6 所示为不同风速区间内不同概率特性风场各 50 条计算得到的应力响应的均值、标准差、偏态和峰 度的平均值. 从图中可知在高斯和非高斯风荷载作用 下,均值、标准差和偏态系数基本相同. 然而,由于高 阶矩的不稳定性,3 种荷载作用下动力响应的峰度系 数有所不同,但整体趋势仍保持一致. 图 7 所示为轮毂高度处平均风速 25 m·s - 1 条件 下,3 种不同概率特性风荷载作用下,控制截面弯曲正 ·1456·

双妙等:全风向来流非高斯风场风机疲劳寿命可靠性分析 ·1457· 70 60 60 (a 95%CI+ ( 。硬化过程 口软化过程 一高斯过程 002010 均值 40 0 0 0-0 95%CI -10 100 200 300 400500 600 20 70 标准差 (b) 60 95%CI+ 0000000 50 0 10 15 20 25 2 m·s) 0 。硬化过程 0 95%C 口软化过程 -10 一高斯过程 100 200 300 400 500 600 tis 70 峰度 60 (c) 95%C+ 0888 50 40 30 偏态 20 10 0 95%CI -1 0 100 200 300 400 500 600 10 15 20 25 tis 丛/m·s) 图5风荷载作用下面内弯曲应力时程.(a)硬化过程:(b)高斯 图63种不同概率特性风荷载作用下控制截面面内弯曲正应力 过程;(c)软化过程 的前四阶统计矩.(a)均值与标准差:(b)偏态和峰度 Fig.5 Time history of the in-plane bending stress of the critical sec- Fig.6 First four statistical moments of the in-plane bending stress of tion under a wind load:(a)hardening process;(b)Gaussian the critical section under wind loads with different probability charac- process;(c)softening process teristics:(a)mean and standard deviation;(b)skewness and kurto- sis 应力的功率谱密度(PSD).由于动力响应时程的时间 102 间隔为0.O5s,根据Nyquist抽样定理,PSD的有效频率 硬化过程 1P 高斯过程 范围为[0,10Hz].由图可知,三种不同概率特性风荷 10°F 2P 软化过程 载作用下弯曲应力的PSD基本相同. 4P 图8是3种不同概率特性风场作用下,由控制截 102 1P:0.327Hz 面面内弯曲正应力的功率谱密度计算得到的带宽系数 's 10E 2P:0.578Hz α、a,和q(式(13)).根据随机振动理论,当g≈0或 3P:1.302Hz 4P:2.915Hz α心,≈1时,应力过程是窄带过程.由图8可知,3种应 10 力过程均表现出明显的宽带特性,且带宽系数基本 10-2 10 10 10 相同. f/Hz 入 图73种不同概率特性风荷载作用下控制截面面内弯曲应力的 a=() 功率谱密度 入2 Fig.7 PSDs of the in-plane bending stress of the critical section un- a2= (oA4)迈, (13) der wind loads with different probability characteristics 9=/1- 4疲劳可靠性分析 入。A2 如图9所示为硬化、高斯和软化风荷载作用下,不 结构的使用寿命通常由疲劳裂纹形成寿命和裂纹 同风速区间内由50个应力时程得到的响应最大值. 扩展寿命两部分组成.本文在考虑全风向来流条件 结合图6和图9可知,在均值和方差近似相等的条件 下,采用线性损伤累积理论对裂纹形成寿命进行计算, 下,软化、高斯和硬化过程的最大值依次减小,且软化 并考虑由结构和荷载两方面不确定性引起的疲劳损伤 过程最大值的离散程度最大 不确定性[2a].同时,采用Paris公式计算疲劳裂纹扩展

双 妙等: 全风向来流非高斯风场风机疲劳寿命可靠性分析 图 5 风荷载作用下面内弯曲应力时程. ( a)硬化过程;( b)高斯 过程;(c)软化过程 Fig. 5 Time history of the in鄄plane bending stress of the critical sec鄄 tion under a wind load: ( a ) hardening process; ( b ) Gaussian process; (c) softening process 应力的功率谱密度(PSD). 由于动力响应时程的时间 间隔为 0郾 05 s,根据 Nyquist 抽样定理,PSD 的有效频率 范围为[0, 10 Hz]. 由图可知,三种不同概率特性风荷 载作用下弯曲应力的 PSD 基本相同. 图 8 是 3 种不同概率特性风场作用下,由控制截 面面内弯曲正应力的功率谱密度计算得到的带宽系数 琢1 、琢2和 q(式(13)). 根据随机振动理论,当 q抑0 或 琢2抑1 时,应力过程是窄带过程. 由图 8 可知,3 种应 力过程均表现出明显的宽带特性,且带宽系数基本 相同. 琢1 = 姿1 (姿0姿2 ) 1 / 2 , 琢2 = 姿2 (姿0姿4 ) 1 / 2 , q = 1 - 姿 2 1 姿0姿2 ì î í ï ï ï ï ï ï ï ï . (13) 如图 9 所示为硬化、高斯和软化风荷载作用下,不 同风速区间内由 50 个应力时程得到的响应最大值. 结合图 6 和图 9 可知,在均值和方差近似相等的条件 下,软化、高斯和硬化过程的最大值依次减小,且软化 过程最大值的离散程度最大. 图 6 3 种不同概率特性风荷载作用下控制截面面内弯曲正应力 的前四阶统计矩. (a)均值与标准差; (b)偏态和峰度 Fig. 6 First four statistical moments of the in鄄plane bending stress of the critical section under wind loads with different probability charac鄄 teristics: (a)mean and standard deviation; ( b) skewness and kurto鄄 sis 图 7 3 种不同概率特性风荷载作用下控制截面面内弯曲应力的 功率谱密度 Fig. 7 PSDs of the in鄄plane bending stress of the critical section un鄄 der wind loads with different probability characteristics 4 疲劳可靠性分析 结构的使用寿命通常由疲劳裂纹形成寿命和裂纹 扩展寿命两部分组成. 本文在考虑全风向来流条件 下,采用线性损伤累积理论对裂纹形成寿命进行计算, 并考虑由结构和荷载两方面不确定性引起的疲劳损伤 不确定性[23] . 同时,采用 Paris 公式计算疲劳裂纹扩展 ·1457·

·1458· 工程科学学报,第39卷,第9期 1.0 00000g 0-000 200o 0.8 。硬化过程 150H 口软化过程 高斯过程 ow 0.4 目日 50 000000-0000000p 10 20 30 0 10 20 30 t/m-s-) m·9 200西 图83种不同概率特性风荷载作用下控制截面面内弯曲应力的 150 频域特性 Fig.8 Frequency-domain characteristics of the in-plane bending stress of the critical section under wind loads with different probability characteristics 50 寿命. 本文将结构的不确定性限定于S-N曲线参数,即 10 20 30 u/(m-s-) 不考虑施工误差等因素:对于外荷载的不确定性,仅考 200 虑由不同应力样本所引起的损伤不确定性,并将疲劳 损伤的概率分布近似为高斯分布. 150 4.1考虑荷载不确定性的裂纹形成寿命 在风荷载作用下,风机结构的面内弯曲应力远大 8 0 于面外弯曲应力和切应力,且面内弯曲应力沿垂直方 50 向(z)与自重引起的压应力方向相同,呈现出典型的单 轴应力状态.根据雨流计数法,对全风向来流作用下z 10 30 30 轴方向正应力S,进行统计分析,得到应力范围S.和相 (m-s-) 应幅值的频率∫后,将零均值的等效应力幅表示为: 图9风荷载作用下控制截面面内弯曲应力的最大值分布.(a) =(∑) 硬化过程:(b)高斯过程:(c)软化过程 (14) Fig.9 Distribution of the maximum responses of the in-plane ben- 为考虑应力均值影响,按Goodman准则将零均值 ding stress of the critical section under wind loads with different prob- 应力过程转化为应力均值为$的应力过程,表示为: ability characteristics:(a)hardening process;(b)Gaussian process; (c)softening process s=-受)广 (15) 式中,S是材料的强度极限. 根据中心极限定理,通常认为疲劳损伤服从高斯 根据S-N曲线,在等效应力幅作用下,风机控制 分布2四.因此,风机结构的疲劳寿命期望E[t]和方 截面裂纹形成时的应力循环次数表示为: 差2[,]分别表示为: N=AS (16) 式中,A和m是S-N曲线参数.按美国公路桥梁设计 (18) 规范(AASHTO LRFD bridge design specifications)中E 12(+1.25) 类曲线取值,即m=3.0和A=36.1×10°MPa3.为考 虑疲劳极限,应力阀值取31MPa,即S<15.5MPa时 式中,4,和σ。是时间,内控制点处的疲劳损伤期望和 不引起疲劳损伤. 标准差,分别表示为: 由线性损伤累积准则可知,时间【内控制截面处 =ua(u.0)f(u.0)dud0, 的疲劳损伤可表示为: (19) Ds。u o。=Jo(u,o)u,o)duda. (17) 式中,和σ:分别表示平均风速和来流风向一定条件 式中,。表示应力过程的均值穿越率 下控制点处的疲劳损伤期望和标准差,∫(“,)是平均

工程科学学报,第 39 卷,第 9 期 图 8 3 种不同概率特性风荷载作用下控制截面面内弯曲应力的 频域特性 Fig. 8 Frequency鄄domain characteristics of the in鄄plane bending stress of the critical section under wind loads with different probability characteristics 寿命. 本文将结构的不确定性限定于 S鄄鄄N 曲线参数,即 不考虑施工误差等因素;对于外荷载的不确定性,仅考 虑由不同应力样本所引起的损伤不确定性,并将疲劳 损伤的概率分布近似为高斯分布. 4郾 1 考虑荷载不确定性的裂纹形成寿命 在风荷载作用下,风机结构的面内弯曲应力远大 于面外弯曲应力和切应力,且面内弯曲应力沿垂直方 向(z)与自重引起的压应力方向相同,呈现出典型的单 轴应力状态. 根据雨流计数法,对全风向来流作用下 z 轴方向正应力 Sz进行统计分析,得到应力范围 Sri和相 应幅值的频率 f i后,将零均值的等效应力幅表示为: S 0 reff = ( 移i f iS 3 ri ) 1 / 3 . (14) 为考虑应力均值影响,按 Goodman 准则将零均值 应力过程转化为应力均值为 Sm的应力过程,表示为: Sreff = S 0 reff·( 1 - Sm S ) u - 1 . (15) 式中,Su是材料的强度极限. 根据 S鄄鄄N 曲线,在等效应力幅作用下,风机控制 截面裂纹形成时的应力循环次数表示为: Nf = AS - m reff . (16) 式中,A 和 m 是 S鄄鄄N 曲线参数. 按美国公路桥梁设计 规范(AASHTO LRFD bridge design specifications) 中 E 类曲线取值,即 m = 3郾 0 和 A = 36郾 1 伊 10 10 MPa 3 . 为考 虑疲劳极限,应力阀值取 31 MPa,即 Sreff < 15郾 5 MPa 时 不引起疲劳损伤. 由线性损伤累积准则可知,时间 t 内控制截面处 的疲劳损伤可表示为: D = v + 0·t Nf . (17) 式中,v + 0 表示应力过程的均值穿越率. 图 9 风荷载作用下控制截面面内弯曲应力的最大值分布. ( a) 硬化过程; (b)高斯过程; (c)软化过程 Fig. 9 Distribution of the maximum responses of the in鄄plane ben鄄 ding stress of the critical section under wind loads with different prob鄄 ability characteristics: (a)hardening process; (b) Gaussian process; (c) softening process 根据中心极限定理,通常认为疲劳损伤服从高斯 分布[24] . 因此,风机结构的疲劳寿命期望 E[ t f ] 和方 差 滓 2 [t f]分别表示为[25] : E[t f] = 1 滋 ( D 1 + 滓 2 D 2滋 ) D t, 滓 2 [t f] = 滓 2 D 滋 3 ( D 1 + 1郾 25 滓 2 D 滋 ) D t 2 ì î í ï ï ï ï . (18) 式中,滋D和 滓D是时间 t 内控制点处的疲劳损伤期望和 标准差,分别表示为: 滋D = 蓦滋d (u,兹)f(u,兹)dud兹, 滓D = 蓦滓d (u,兹)f(u,兹)dud兹 ì î í ïï ïï . (19) 式中,滋d和 滓d分别表示平均风速和来流风向一定条件 下控制点处的疲劳损伤期望和标准差,f( u, 兹)是平均 ·1458·

双妙等:全风向来流非高斯风场风机疲劳寿命可靠性分析 ·1459· 风速和来流风向的联合概率密度函数. 弯曲应力为压应力,则该点的均值应力为压应力, 图10是在考虑平均风速概率密度函数条件 由于Goodman准则不对均值应力是压应力的应力 下,来流风向为0°时,控制截面16个控制点,由不 过程修正,因此在应力范围相同情况下该点的疲劳 同风速区间内50条应力时程分别计算得到的疲劳 损伤小于0°点.此外,由于面外弯曲应力方向在 损伤的均值.由于0°点处的面内弯曲应力为拉应 [270°,360°]范围内与面内弯曲应力方向相同,因 力,则该点的均值应力为拉应力,根据Goodman准 此均值应力大于[0°,90]范围内相对应各点,疲 则可知该点的疲劳损伤最大.而180°点处的面内 劳损伤也相应增加. 37.5l410925 硬化过程 (a) 一高斯过程 14 ·软化过程 315.0 45.09 b) 12 ■·硬化过程 ●一高斯过程 292.59 67.5° 10 ▲一软化过程 4 /D 270.0° 90.0 6 247.5 112.5° 风 225.0° 135.0 045 90135180225270315360 202.50 180.0° 157.59 180° 角度) 图10来流风向为0°时,控制截面不同控制点(a)和其疲劳损伤(b) Fig.10 Critical points (a)and fatigue damage (b)at different points of the critical section with a wind direction of 图11是在考虑平均风速和来流风向联合概率密 引起的裂纹形成寿命的期望.结合图2可知,工程所 度函数条件下,由式(18)和式(19)计算得到的3种不 处地区主导风向为WNW,因此在292.5°处的疲劳裂 同概率特性风荷载作用下控制截面各点处由正应力S, 纹形成寿命最小 0单位: ……硬化过程 a 337.5120 22.50 一高斯过程 120 100 315.09 ·软化过程 b 0 45.0° I00 C60 292.5 40 67.5 90 30 60 270.0° 90.0° ■·硬化过程 一高斯过程 247.59 ▲·软化过程 112.5 0 225.0° 135.0 004590135180225270315360 角度() 202.50 157.5 180.0° 图11考虑平均风速和来流风向时,控制截面不同控制点处(a)和其疲劳寿命(b) Fig.11 Critical points (a)and fatigue life (b)at different points of the critical section considering the mean wind speed and wind direction 表1为三种不同概率特性风荷载作用下,292.5° 表1考虑全风向来流条件下,292.5°处的疲劳裂纹形成寿命的期 处的疲劳裂纹形成寿命的期望和标准差.由表可知, 望和标准差 与高斯风荷载相比较,硬化风荷载作用下控制点的裂 Table 1 Mean and standard deviations of the fatigue crack initiation life at a critical point of 292.5 considering a full-direction inflow 纹形成寿命增加约5a,软化风荷载作用下则减小约5 a.因此,需要考虑风荷载的非高斯性对风机结构裂纹 荷载类型 E[t]/a o[t]/a 形成寿命的影响.分析可知,对于地形复杂地区,由于 硬化过程 70.74 6.4×10-3 风场呈现非高斯性,高斯假设会过高的估计风机疲劳 高斯过程 65.42 5.8×10-3 寿命,考虑软化非高斯影响更为合理.此外,由样本引 软化过程 60.15 4.1×10-3 起的寿命标准差远小于寿命的期望,可不考虑样本不 确定性对裂纹形成寿命的影响. 4.2考虑结构参数不确定性的裂纹形成寿命 S-N曲线参数具有很强的离散性2],设计规范中

双 妙等: 全风向来流非高斯风场风机疲劳寿命可靠性分析 风速和来流风向的联合概率密度函数. 图 10 是 在 考 虑 平 均 风 速 概 率 密 度 函 数 条 件 下,来流风向为 0毅时,控制截面 16 个控制点,由不 同风速区间内 50 条应力时程分别计算得到的疲劳 损伤的均值. 由于 0毅 点处的面内弯曲应力为拉应 力,则该点的均值应力为拉应力,根据 Goodman 准 则可知该点的疲劳损伤最大. 而 180毅 点处的面内 弯曲应力为压应力,则该点的均值应力为压应力, 由于 Goodman 准则不对均值应力是压应力的应力 过程修正,因此在应力范围相同情况下该点的疲劳 损伤小于 0毅 点. 此外,由于面外弯曲应力方向在 [270毅, 360毅] 范围内与面内弯曲应力方向相同,因 此均值应力大于[ 0毅, 90毅] 范围内相对应各点,疲 劳损伤也相应增加. 图 10 来流风向为 0毅时,控制截面不同控制点(a)和其疲劳损伤(b) Fig. 10 Critical points (a) and fatigue damage (b) at different points of the critical section with a wind direction of 0毅 图 11 是在考虑平均风速和来流风向联合概率密 度函数条件下,由式(18)和式(19)计算得到的 3 种不 同概率特性风荷载作用下控制截面各点处由正应力 Sz 引起的裂纹形成寿命的期望. 结合图 2 可知,工程所 处地区主导风向为 WNW,因此在 292郾 5毅处的疲劳裂 纹形成寿命最小. 图 11 考虑平均风速和来流风向时,控制截面不同控制点处(a)和其疲劳寿命(b) Fig. 11 Critical points (a) and fatigue life (b) at different points of the critical section considering the mean wind speed and wind direction 表 1 为三种不同概率特性风荷载作用下,292郾 5毅 处的疲劳裂纹形成寿命的期望和标准差. 由表可知, 与高斯风荷载相比较,硬化风荷载作用下控制点的裂 纹形成寿命增加约 5 a,软化风荷载作用下则减小约 5 a. 因此,需要考虑风荷载的非高斯性对风机结构裂纹 形成寿命的影响. 分析可知,对于地形复杂地区,由于 风场呈现非高斯性,高斯假设会过高的估计风机疲劳 寿命,考虑软化非高斯影响更为合理. 此外,由样本引 起的寿命标准差远小于寿命的期望,可不考虑样本不 确定性对裂纹形成寿命的影响. 表 1 考虑全风向来流条件下,292郾 5毅处的疲劳裂纹形成寿命的期 望和标准差 Table 1 Mean and standard deviations of the fatigue crack initiation life at a critical point of 292郾 5毅 considering a full鄄direction inflow 荷载类型 E[t f] / a 滓[t f] / a 硬化过程 70郾 74 6郾 4 伊 10 - 3 高斯过程 65郾 42 5郾 8 伊 10 - 3 软化过程 60郾 15 4郾 1 伊 10 - 3 4郾 2 考虑结构参数不确定性的裂纹形成寿命 S鄄鄄N 曲线参数具有很强的离散性[26] ,设计规范中 ·1459·

·1460· 工程科学学报,第39卷,第9期 通常选取存活率为95%以上的S-N曲线进行设计 P=1-(B). (21) 同时,对于已带裂纹构件或重要结构,Sam等建议将 式中,中(·)是标准高斯分布的概率分布函数:B是可 现有规范参数降低50%取值,严重低估裂纹形成寿命 靠性指数,表示为: 现有研究表明,S-N曲线参数A可近似认为服从 B=(a+)-3In Sm-In N (22) 对数高斯分布2].文献[7]中给出了按对数高斯分布 √+ 统计得到的参数A的均值和变异系数.本文按AASH- 其中,入、入45和分别表示为: T0规范中E类曲线取值,即4=65.9×10°MPa3, 84=0.26. 4sln,-号入=ln4~9 在Miner准则中,认为疲劳损伤D累积到界限△ 54=√1n(1+6),54=√1n(1+8).(23) 时,控制截面产生裂纹,工程中通常取△=l.Chung 图12和图13表示考虑来流风向和平均风速联合 等[2]指出需要考虑累积界限△的不确定性,并认为△ 概率密度条件下,使用年限分别为50a和100a时,控 也服从对数高斯分布,取44=1和6=0.3.由于考虑 制截面不同控制点裂纹形成阶段的失效概率.与图11 了累积界限的不确定性,控制截面裂纹形成时的应力 所示规律相同,292.5°处的50a和100a的失效概率最 循环次数表示为: 大.由图12(b)和图13(b)可知,硬化和软化非高斯性 N,=4.4 对裂纹形成阶段失效概率均有影响,且与考虑荷载不 (20) 确定性时类似,对于地形复杂地区,考虑软化非高斯影 裂纹形成的失效概率可表示为: 响更为合理 (a) 337.5 12%)22.5° 硬化过程 10 高斯过程 12 315.0 .8 45.0 软化过程 b 6 292.5° 4 67.5 8 270.0 190.0 ■硬化过程 一高斯过程 247.5 112.5 2 ▲·软化过程 0 225.0 135.0° 0 4590135180225270315360 角度/) 202.5° 157.5 180.0° 图12使用年限为50a时不同控制点裂纹形成阶段的失效概率 Fig.12 Critical points (a)and failure probabilities (b)of different critical points with service lives of 50 a in the crack initiation stage …硬化过程 (a) 337.525(%) 22.50 一高斯过程 2 一··软化过程 b 315.09 45.0° 1 20 292.5 10 67.5 15 源★ 10 ◆ 270.0 90.0° ■·硬化过程 5 ●一高斯过程 247.5 112.5° ▲·软化过程 06 4590135180225270315360 225.0° 135.0 角度/() 202.50 1575 180.0° 图13使用年限为100a时不同控制点(a)及其裂纹形成阶段的失效概率(b) Fig.13 Critical points (a)and failure probabilities (b)of different critical points with service lives of 100 a in the crack initiation stage 4.3裂纹扩展寿命 中,结构仍有疲劳寿命[).但现有设计规范不考虑疲 从控制截面裂纹形成到裂纹扩展至临界长度过程 劳裂纹扩展寿命,因此需根据断裂力学理论计算裂纹

工程科学学报,第 39 卷,第 9 期 通常选取存活率为 95% 以上的 S鄄鄄 N 曲线进行设计. 同时,对于已带裂纹构件或重要结构,Stam 等 [27]建议将 现有规范参数降低50%取值,严重低估裂纹形成寿命. 现有研究表明,S鄄鄄N 曲线参数 A 可近似认为服从 对数高斯分布[28] . 文献[7]中给出了按对数高斯分布 统计得到的参数 A 的均值和变异系数. 本文按 AASH鄄 TO 规范中 E 类曲线取值,即 滋A = 65郾 9 伊 10 10 MPa 3 , 啄A = 0郾 26. 在 Miner 准则中,认为疲劳损伤 D 累积到界限 驻 时,控制截面产生裂纹,工程中通常取 驻 = 1. Chung 等[28]指出需要考虑累积界限 驻 的不确定性,并认为 驻 也服从对数高斯分布,取 滋驻 = 1 和 啄驻 = 0郾 3. 由于考虑 了累积界限的不确定性,控制截面裂纹形成时的应力 循环次数表示为: Nf = A·驻 S 3 reff . (20) 裂纹形成的失效概率可表示为: Pf = 1 - 椎(茁). (21) 式中,椎(·)是标准高斯分布的概率分布函数;茁 是可 靠性指数,表示为: 茁 = (姿驻 + 姿A ) - 3ln Sreff - ln N 灼 2 驻 + 灼 2 A . (22) 其中,姿驻 、姿A 、灼驻和 灼A分别表示为: 姿A = ln 滋A - 灼 2 A 2 ,姿驻 = ln 滋驻 - 灼 2 驻 2 , 灼A = ln (1 + 啄 2 A ),灼驻 = ln (1 + 啄 2 驻 ). (23) 图 12 和图 13 表示考虑来流风向和平均风速联合 概率密度条件下,使用年限分别为 50 a 和 100 a 时,控 制截面不同控制点裂纹形成阶段的失效概率. 与图 11 所示规律相同,292郾 5毅处的 50 a 和 100 a 的失效概率最 大. 由图 12(b)和图 13(b)可知,硬化和软化非高斯性 对裂纹形成阶段失效概率均有影响,且与考虑荷载不 确定性时类似,对于地形复杂地区,考虑软化非高斯影 响更为合理. 图 12 使用年限为 50 a 时不同控制点裂纹形成阶段的失效概率 Fig. 12 Critical points (a) and failure probabilities (b) of different critical points with service lives of 50 a in the crack initiation stage 图 13 使用年限为 100 a 时不同控制点(a)及其裂纹形成阶段的失效概率(b) Fig. 13 Critical points (a) and failure probabilities (b) of different critical points with service lives of 100 a in the crack initiation stage 4郾 3 裂纹扩展寿命 从控制截面裂纹形成到裂纹扩展至临界长度过程 中,结构仍有疲劳寿命[29] . 但现有设计规范不考虑疲 劳裂纹扩展寿命,因此需根据断裂力学理论计算裂纹 ·1460·

双妙等:全风向来流非高斯风场风机疲劳寿命可靠性分析 ·1461· 扩展寿命. 通常根据裂纹扩展速率将裂纹扩展分为失稳扩展 a=a。= 1. 0.707Ke) (26) 和亚临界扩展两种.图14所示为裂纹扩展速率与应 式中,K是应力强度因子临界值,由夏比冲击实验 力强度因子范围△K的变化曲线,当△K小于其临界值 (CVN)确定;S是塔体控制截面处的最大拉应力, Kc时,裂纹扩展速率da/dW可由Paris公式确定: S.={SCF[Sm/2-P/A]+o,},SCF表示应力集 =C(K) da (24) 中系数,取1.5,P/A表示结构自重引起的压应力,0, 和σ分别表示钢材的屈服强度和极限强度,本文选取 式中,C和m是材料参数;△K是裂纹尖端的等效交 的ASTM A36型钢材,材料参数见图14. 变应力强度因子范围 根据平均风速的概率密度函数和应力循环次数 10- (式(25)),可得到裂纹扩展寿命,表示为: ASTM A36钢 105 0,=250 MPa ,=400 MPa t= (27) D.(u.0)f(u.0)dud0 10-6 Kic=100 MPa.m2 式中,D(u,0)表示平均风速和来流风向一定条件下 107E 的断裂疲劳损伤,即: da/dN=C△Km 108 m=3 C-1.65x10-m-MPa3.m号 0侧是 (28) 10 10P 10 102 103 图15是初始裂纹2a:=28时,考虑全风向来流和 △K/MPa·ma 平均风速联合概率密度条件下,控制截面不同点处的 图14疲劳裂纹扩展速率曲线 疲劳裂纹扩展寿命.由于裂纹扩展寿命仅受拉应力影 Fig.14 Fatigue crack growth rates 响,且工程所处位置的主导风向为WNW,因此292.5° 位置处的裂纹扩展寿命最小.结合图12和图13可 根据Paris公式,疲劳裂纹扩展的应力循环次数 知,裂纹形成和扩展阶段的失效位置相同,这是由于裂 N可表示为: 纹扩展是受拉应力控制,Goodman准则仅对均值拉应 a2-a2 Ni= (25) 力进行修正,因此结构最易发生疲劳失效的位置在主 C(F△S√m)".(1-m/2) 导风向上.此外,由图15(b)可知,裂纹扩展寿命对风 式中:a,和a,分别表示最终裂纹尺寸和初始裂纹长度, 荷载的非高斯性并不敏感,这主要是由于在计算临界 通常认为初始裂纹贯穿结构壁厚(8),即2a,=2δ;AS 裂纹尺寸时(式(26)),最大拉应力S.中材料的屈服 是有效应力范围,△S=S:F是关于裂纹几何形状和 强度远大于等效应力幅S,从而使由荷载非高斯性 尺寸的量纲一参数,当裂纹长度与构件宽度比值α≤ 引起的最大拉应力差异不明显.因此,在裂纹扩展阶 0.4时,取F=1.120 段可忽略非高斯性影响 根据英国标准委员会颁布的《金属结构中缺陷验 收评定方法导则》(BS7910)中1A等级,临界裂纹尺 5结论 寸a.(结构破坏时,ar=ae)表示为: (1)软化非高斯、高斯和硬化非高斯3种不同概 a 337.5 g单位25 硬化过程 高斯过程 315.0° 一··软化过程 45.0° ■硬化过程 ●一高斯过程 ▲·软化过程 292.52 67.5 270.0° 190.0° 247.5 112.5o 225.0° 135.0° 00459013518025270315360 角度) 202.50 180.0° 157.5 图15初始裂纹2a:=28条件下控制截面不同控制点(a)及其裂纹扩展寿命(b) Fig.15 Critical points (a)and crack propagation lives (b)at different points of the critical section with an initial crack of 2=2

双 妙等: 全风向来流非高斯风场风机疲劳寿命可靠性分析 扩展寿命. 通常根据裂纹扩展速率将裂纹扩展分为失稳扩展 和亚临界扩展两种. 图 14 所示为裂纹扩展速率与应 力强度因子范围 驻K 的变化曲线,当 驻K 小于其临界值 KIC时,裂纹扩展速率 da / dN 可由 Paris 公式确定: da dN = C (驻Kreff) m . (24) 式中,C 和 m 是材料参数;驻Kreff是裂纹尖端的等效交 变应力强度因子范围. 图 14 疲劳裂纹扩展速率曲线 Fig. 14 Fatigue crack growth rates 根据 Paris 公式,疲劳裂纹扩展的应力循环次数 Nif可表示为: Nif = a 1 - m/ 2 f - a 1 - m/ 2 i C (F驻S 仔) m·(1 - m/ 2) . (25) 图 15 初始裂纹 2ai = 2啄 条件下控制截面不同控制点(a)及其裂纹扩展寿命(b) Fig. 15 Critical points (a) and crack propagation lives (b) at different points of the critical section with an initial crack of 2ai = 2啄 式中:af和 ai分别表示最终裂纹尺寸和初始裂纹长度, 通常认为初始裂纹贯穿结构壁厚( 啄),即 2ai = 2啄;驻S 是有效应力范围,驻S = Sreff;F 是关于裂纹几何形状和 尺寸的量纲一参数,当裂纹长度与构件宽度比值 琢臆 0郾 4 时,取 F = 1郾 12 [30] . 根据英国标准委员会颁布的《金属结构中缺陷验 收评定方法导则》 (BS 7910) 中 1A 等级,临界裂纹尺 寸 ac(结构破坏时,af = ac)表示为: af = ac = 1 仔 ·( 0郾 707KIC St, ) max 2 . (26) 式中,KIC 是应力强度因子临界值,由夏比冲击实验 (CVN)确定;St,max是塔体控制截面处的最大拉应力, St,max = { SCF[ Sreff / 2 - P/ A] + 滓y },SCF 表示应力集 中系数,取 1郾 5,P/ A 表示结构自重引起的压应力,滓y 和 滓u分别表示钢材的屈服强度和极限强度,本文选取 的 ASTM A36 型钢材,材料参数见图 14. 根据平均风速的概率密度函数和应力循环次数 (式(25)),可得到裂纹扩展寿命,表示为: t = 1 蓦 ^Di(u,兹)f(u,兹)dud兹 . (27) 式中, ^Di(u,兹) 表示平均风速和来流风向一定条件下 的断裂疲劳损伤,即: ^Di(u,兹) = v + 0 Nif ·t. (28) 图 15 是初始裂纹 2ai = 2啄 时,考虑全风向来流和 平均风速联合概率密度条件下,控制截面不同点处的 疲劳裂纹扩展寿命. 由于裂纹扩展寿命仅受拉应力影 响,且工程所处位置的主导风向为 WNW,因此 292郾 5毅 位置处的裂纹扩展寿命最小. 结合图 12 和图 13 可 知,裂纹形成和扩展阶段的失效位置相同,这是由于裂 纹扩展是受拉应力控制,Goodman 准则仅对均值拉应 力进行修正,因此结构最易发生疲劳失效的位置在主 导风向上. 此外,由图 15(b)可知,裂纹扩展寿命对风 荷载的非高斯性并不敏感,这主要是由于在计算临界 裂纹尺寸时(式(26)),最大拉应力 St,max中材料的屈服 强度远大于等效应力幅 Sreff,从而使由荷载非高斯性 引起的最大拉应力差异不明显. 因此,在裂纹扩展阶 段可忽略非高斯性影响. 5 结论 (1)软化非高斯、高斯和硬化非高斯 3 种不同概 ·1461·

·1462· 工程科学学报,第39卷,第9期 率特性风荷载作用下风机结构动力响应的时域和频域 collapsed large wind turbine tower.Eng Fail Anal,2011,18 特性基本相同,功率谱密度也近似相等.但动力响应 (1):295 的最大值依次减小,且响应的非高斯性较风荷载的非 [15]Jonkman J,Butterfield S,Musial W,et al.Definition of a 5MW 高斯性减弱. Reference Wind Turbine for Offshore System Development.Techni- cal Report,NREL/TP-500-38060,National Renewable Energy (2)裂纹形成阶段对风荷载的非高斯性较为敏 Laboratory,2009 感,因此风机重要部件需考虑荷载的非高斯性影响:裂 [16]Akdag S A,Dinler A.A new method to estimate Weibull param- 纹扩展阶段对风荷载的非高斯性不敏感,可忽略非高 eters for wind energy applications.Energy Coners Manage, 斯性对裂纹扩展寿命的影响. 2009,50(7):1761 (3)在考虑平均风速和来流风向联合概率密度条 [17]Clifton A.135 m Meteorological Towers at the National Wind 件下,控制截面的疲劳失效出现在最大来流风向位置, Technology Center.NREL Report TP-500-55915,2016 且裂纹形成和扩展阶段的失效位置相同. [18]Ding J,Chen X Z.Moment-based translation model for hardening non-Gaussian response processes.J Eng Mech,2016,142(2): 06015006-】 参考文献 [19]Grigoriu M.Spectralrepresentation for a class of non-Gaussian processes.J Eng Mech,2004,130(5):541 [1]Saravia C M,Machado S P,Cortinez V H.A composite beam fi- [20]Shields M D,Deodatis G.A simple and efficient methodology to nite element for multibody dynamics:application to large wind tur- approximate a general non-Gaussian stationary stochastic vector bine modeling.Eng Struct,2013,56:1164 process by a translation process with applications in wind velocity [2]Chen X B,Li J,Chen J Y.Wind-induced response analysis of a simulation.Probabilist Eng Mech,2013,31:19 wind turbine tower including the blade-tower coupling effect. [21]Shields M D,Deodatis G,Bocchini P.A simple and efficient Zhejiang Unin-SCIENCE A,2009,10(11):1573 methodology to approximate a general non-Gaussian stationary [3]Jonkman J M.Buhl Jr M L.FAST User's Guide.Technical Report stochastic process by a translation process.Probabilist Eng Mech, No.NREL/EL-500-38230,2005 2011,26(4):511 [4]Larsen T J,Hansen A M.How 2 HAWC2 the User's Manual.Risg [22]Bocchini P.Probabilistic Approaches in Ciril Engineering:Genera- National Laboratory,2007 tion of Random Fields and Structural Identification with Genetic [5]Do T Q,van de Lindt J W,Mahmoud H.Fatigue life fragilities Algorithms[Dissertation ]Bologna:Universita di Bologna, and performance-based design of wind turbine tower base connec- 2008 tions..J Struct Eng,2014,141(7):04014183-1 [23]Bengtsson A,Rychlik I.Uncertainty in fatigue life prediction of [6]Do T Q,Mahmoud H,van de Lindt J W.Fatigue life of wind tur- structures subject to Gaussian loads.Probabilist Eng Mech, bine tower bases throughout Colorado.I Perform Constr Fac, 2009,24(2):224 2015,29(4):04014109-1 [24]Low Y M.Variance of the fatigue damage due to a Gaussian nar- [7]Dawood M,Goyal R,Dhonde H,et al.Fatiguelife assessment of rowband process.Struct Saf,2012,34(1):381 cracked high-Mast Illumination Poles.J Perform Constr Fac, [25]Ding J,Chen X Z,Zuo D L,et al.Fatigue life assessment of 2014,28(2):311 traffic-signal support structures from an analytical approach and [8]Repetto M P,Solari G.Closed-form prediction of the alongwind- long-term vibration monitoring data.Struct Eng,2016,142 induced fatigue of structures.J Struct Eng,2012,138(9):1149 (6):04016017-1 [9]Repetto M P.Solari G.Closed form solution of the alongwind-in- [26]Hanaki S,Yamashita M,Uchida H,et al.On stochastic evalua- duced fatigue damage to structures.Eng Struct,2009.31(10): tion of S-N data based on fatigue strength distribution.Int Fa- 2414 tigue,2010,32(3):605 [10]Repetto M P.Cycle counting methods for bi-modal stationary [27]Stam A,Richman N.Pool C,et al.Fatigue Life of Steel Base Gaussian processes.Probabilist Eng Mech,2005,20(3):229 Plate to Pole Connections for Traffic Structures.Austin:Center [11]Gong K M,Chen X Z.Influence of non-Gaussian wind charac- for Transportation Research,University of Texas at Austin,2011 teristics on wind turbine extreme response.Eng Struct,2014, [28]Chung H,Manuel L,Frank K H.Optimal Inspection of Fracture- 59:727 Critical Steel Trapesoidal Girders.Austin:Center for Transporta- [12]Gong K M,Ding J,Chen X Z.Estimation of long-term extreme tion Research,University of Texas at Austin,2003 response of operational and parked wind turbines:validation and [29]Mahmoud H N.Dexter R J.Propagation rate of large cracks in some new insights.Eng Struct,2014,81:135 stiffened panels under tension loading.Mar Struct,2005,18 [13]Ding J.Gong K M,Chen X Z.Comparison of statistical extrapo- (3):265 lation methods for the evaluation of long-term extreme response of [30]Dowling N E.Mechanical Beharior of Materials:Engineering wind turbine.Eng Struct,2013,57:100 Methods for Deformation,Fracture,and Fatigue.4th Ed.Bos- [14]Chou JS,Tu W T.Failure analysis and risk management of a ton:Pearson,2012

工程科学学报,第 39 卷,第 9 期 率特性风荷载作用下风机结构动力响应的时域和频域 特性基本相同,功率谱密度也近似相等. 但动力响应 的最大值依次减小,且响应的非高斯性较风荷载的非 高斯性减弱. (2)裂纹形成阶段对风荷载的非高斯性较为敏 感,因此风机重要部件需考虑荷载的非高斯性影响;裂 纹扩展阶段对风荷载的非高斯性不敏感,可忽略非高 斯性对裂纹扩展寿命的影响. (3)在考虑平均风速和来流风向联合概率密度条 件下,控制截面的疲劳失效出现在最大来流风向位置, 且裂纹形成和扩展阶段的失效位置相同. 参 考 文 献 [1] Saravia C M, Machado S P, Cort侏nez V H. A composite beam fi鄄 nite element for multibody dynamics: application to large wind tur鄄 bine modeling. Eng Struct, 2013, 56: 1164 [2] Chen X B, Li J, Chen J Y. Wind鄄induced response analysis of a wind turbine tower including the blade鄄tower coupling effect. J Zhejiang Univ鄄SCIENCE A, 2009, 10(11): 1573 [3] Jonkman J M, Buhl Jr M L. FAST User蒺s Guide. Technical Report No. NREL / EL鄄鄄500鄄鄄38230, 2005 [4] Larsen T J, Hansen A M. How 2 HAWC2 the User蒺s Manual. Ris覬 National Laboratory, 2007 [5] Do T Q, van de Lindt J W, Mahmoud H. Fatigue life fragilities and performance鄄based design of wind turbine tower base connec鄄 tions. J Struct Eng, 2014, 141(7): 04014183鄄1 [6] Do T Q, Mahmoud H, van de Lindt J W. Fatigue life of wind tur鄄 bine tower bases throughout Colorado. J Perform Constr Fac, 2015, 29(4): 04014109鄄1 [7] Dawood M, Goyal R, Dhonde H, et al. Fatiguelife assessment of cracked high鄄Mast Illumination Poles. J Perform Constr Fac, 2014, 28(2): 311 [8] Repetto M P, Solari G. Closed鄄form prediction of the alongwind鄄 induced fatigue of structures. J Struct Eng, 2012, 138(9): 1149 [9] Repetto M P, Solari G. Closed form solution of the alongwind鄄in鄄 duced fatigue damage to structures. Eng Struct, 2009, 31 (10): 2414 [10] Repetto M P. Cycle counting methods for bi鄄modal stationary Gaussian processes. Probabilist Eng Mech, 2005, 20(3): 229 [11] Gong K M, Chen X Z. Influence of non鄄Gaussian wind charac鄄 teristics on wind turbine extreme response. Eng Struct, 2014, 59: 727 [12] Gong K M, Ding J, Chen X Z. Estimation of long鄄term extreme response of operational and parked wind turbines: validation and some new insights. Eng Struct, 2014, 81: 135 [13] Ding J, Gong K M, Chen X Z. Comparison of statistical extrapo鄄 lation methods for the evaluation of long鄄term extreme response of wind turbine. Eng Struct, 2013, 57: 100 [14] Chou J S, Tu W T. Failure analysis and risk management of a collapsed large wind turbine tower. Eng Fail Anal, 2011, 18 (1): 295 [15] Jonkman J, Butterfield S, Musial W, et al. Definition of a 5MW Reference Wind Turbine for Offshore System Development. Techni鄄 cal Report, NREL / TP鄄鄄500鄄鄄38060, National Renewable Energy Laboratory,2009 [16] Akdag姚 S A, Dinler A. A new method to estimate Weibull param鄄 eters for wind energy applications. Energy Convers Manage, 2009, 50(7): 1761 [17] Clifton A. 135 m Meteorological Towers at the National Wind Technology Center. NREL Report TP鄄鄄500鄄鄄55915, 2016 [18] Ding J,Chen X Z. Moment鄄based translation model for hardening non鄄Gaussian response processes. J Eng Mech, 2016, 142(2): 06015006鄄1 [19] Grigoriu M. Spectralrepresentation for a class of non鄄Gaussian processes. J Eng Mech, 2004, 130(5): 541 [20] Shields M D, Deodatis G. A simple and efficient methodology to approximate a general non鄄Gaussian stationary stochastic vector process by a translation process with applications in wind velocity simulation. Probabilist Eng Mech, 2013, 31: 19 [21] Shields M D, Deodatis G, Bocchini P. A simple and efficient methodology to approximate a general non鄄Gaussian stationary stochastic process by a translation process. Probabilist Eng Mech, 2011, 26(4): 511 [22] Bocchini P. Probabilistic Approaches in Civil Engineering: Genera鄄 tion of Random Fields and Structural Identification with Genetic Algorithms [ Dissertation ]. Bologna: Universit伽 di Bologna, 2008 [23] Bengtsson A, Rychlik I. Uncertainty in fatigue life prediction of structures subject to Gaussian loads. Probabilist Eng Mech, 2009, 24(2): 224 [24] Low Y M. Variance of the fatigue damage due to a Gaussian nar鄄 rowband process. Struct Saf, 2012, 34(1): 381 [25] Ding J, Chen X Z, Zuo D L, et al. Fatigue life assessment of traffic鄄signal support structures from an analytical approach and long鄄term vibration monitoring data. J Struct Eng, 2016, 142 (6): 04016017鄄1 [26] Hanaki S, Yamashita M, Uchida H, et al. On stochastic evalua鄄 tion of S鄄鄄N data based on fatigue strength distribution. Int J Fa鄄 tigue, 2010, 32(3): 605 [27] Stam A, Richman N, Pool C, et al. Fatigue Life of Steel Base Plate to Pole Connections for Traffic Structures. Austin: Center for Transportation Research, University of Texas at Austin, 2011 [28] Chung H, Manuel L, Frank K H. Optimal Inspection of Fracture鄄 Critical Steel Trapezoidal Girders. Austin: Center for Transporta鄄 tion Research, University of Texas at Austin, 2003 [29] Mahmoud H N, Dexter R J. Propagation rate of large cracks in stiffened panels under tension loading. Mar Struct, 2005, 18 (3): 265 [30] Dowling N E. Mechanical Behavior of Materials: Engineering Methods for Deformation, Fracture, and Fatigue. 4th Ed. Bos鄄 ton: Pearson, 2012 ·1462·

点击下载完整版文档(PDF)VIP每日下载上限内不扣除下载券和下载次数;
按次数下载不扣除下载券;
24小时内重复下载只扣除一次;
顺序:VIP每日次数-->可用次数-->下载券;
已到末页,全文结束
相关文档

关于我们|帮助中心|下载说明|相关软件|意见反馈|联系我们

Copyright © 2008-现在 cucdc.com 高等教育资讯网 版权所有