基于概率密度演化的風(fēng)機(jī)基礎(chǔ)疲勞可靠度計(jì)算

趙儉斌，王凱威，王一達(dá)，付兵?

（1.沈陽(yáng)建筑大學(xué) 土木工程學(xué)院，遼寧 沈陽(yáng) 110168；2.上海杰地建筑設(shè)計(jì)有限公司，上海 200085）

在過去20 年，風(fēng)電持續(xù)高速發(fā)展[1-3]，到2018 年底，全球風(fēng)電累計(jì)裝機(jī)容量已突破600 GW，新增裝機(jī)容量共53.9 GW.風(fēng)機(jī)受到長(zhǎng)期的動(dòng)力荷載作用，疲勞問題突出，需要準(zhǔn)確預(yù)測(cè)風(fēng)機(jī)結(jié)構(gòu)的疲勞損傷，尤其是在塔底動(dòng)應(yīng)力很高的地方.

影響結(jié)構(gòu)疲勞損傷的因素甚多，其中荷載是最主要的方面，風(fēng)荷載本質(zhì)上具有不可忽略的隨機(jī)性，因此，采用可靠度的分析方法分析結(jié)構(gòu)疲勞損傷是一種自然的選擇.

目前，可靠度的求解方法包括：Monte Carlo[4，5]、一次二階矩[5]、響應(yīng)面法[5]等.基于一次二階矩理論的可靠度分析方法的主要目標(biāo)在于尋求隨機(jī)結(jié)構(gòu)響應(yīng)的二階矩統(tǒng)計(jì)量.在獲得二階矩統(tǒng)計(jì)量之后，通過假定結(jié)構(gòu)響應(yīng)服從正態(tài)分布計(jì)算結(jié)構(gòu)的使用可靠度.響應(yīng)面法的基本思想是將功能函數(shù)的輸入和輸出變量表示為標(biāo)準(zhǔn)正態(tài)分布變量的多項(xiàng)式，而各變量的系數(shù)通過配點(diǎn)法（collocation points）確定，最后由功能方程得到可靠指標(biāo)和失效概率.Monte Carlo 法雖然適用性較強(qiáng)，但以過高的計(jì)算量為代價(jià).近年來陳建兵和李杰[7-9]發(fā)展了一類隨機(jī)結(jié)構(gòu)反應(yīng)的概率密度演化方法，利用這種方法，可以精確定量地給出結(jié)構(gòu)反應(yīng)的演化概率密度曲線族，由此，可方便地根據(jù)指定的位移反應(yīng)限值，直接計(jì)算給出隨機(jī)結(jié)構(gòu)在隨機(jī)荷載作用下的可靠度.

對(duì)于風(fēng)機(jī)這種新的結(jié)構(gòu)形式，其受到的風(fēng)荷載突出，隨機(jī)風(fēng)引起的疲勞損傷問題也更突出，但是有關(guān)風(fēng)機(jī)結(jié)構(gòu)的疲勞可靠度分析成果還很少.

本文基于概率密度演化思想，構(gòu)造一個(gè)虛擬隨機(jī)過程，使得隨機(jī)結(jié)構(gòu)時(shí)域內(nèi)的疲勞損傷為該虛擬隨機(jī)過程的截口隨機(jī)變量.進(jìn)而，建立概率密度演化方程并求解出隨機(jī)結(jié)構(gòu)疲勞損傷的概率密度，在安全域內(nèi)積分給出結(jié)構(gòu)的疲勞可靠度.

1 風(fēng)荷載的數(shù)值模擬

設(shè)x，y，z 為空間中的一個(gè)點(diǎn)，其中z 是離地面的高度，x 是橫向風(fēng)向，y 是順風(fēng)向.在實(shí)踐中，為了簡(jiǎn)化概念，風(fēng)速波動(dòng)可以在平面y=0 中表征[10].當(dāng)只考慮豎向相關(guān)性時(shí)，風(fēng)速場(chǎng)可以寫成：

式中：符號(hào)θ 表示相應(yīng)數(shù)量的隨機(jī)性，V（θ；z，t）是高度z 處的風(fēng)速（m/s），是高度z 處的平均風(fēng)速（m/s），是風(fēng)速的波動(dòng)分量（m/s），假設(shè)為正常的零均值穩(wěn)態(tài)隨機(jī)過程.平均風(fēng)速可用冪函數(shù)表示為[11-12]：

式中：·表示導(dǎo)數(shù)，χk+1為諧波能量調(diào)整系數(shù)，φj，k+1是標(biāo)準(zhǔn)特征向量Фj的k+1 個(gè)元素，λj和Фj分別是相關(guān)矩陣R 的特征值和特征向量[13]RФj=λjФj.R=[ρij]（N×N）.

無量綱的隨機(jī)變量組{ξ1（θ），…，ξr（θ）}滿足下列關(guān)系:

式中：δij為Kronecker-delta 符號(hào).

式（3）中：

只考慮其中λzj和fzj（z）是特征值和相關(guān)函數(shù)RUz（z1，z2）的相應(yīng)本征函數(shù)，可由Fredholm 積分方程得到

式中：δij為Kronecker-delta 符號(hào).本文中rz=5.

2 疲勞可靠度分析的概率密度演化方法

在對(duì)風(fēng)速進(jìn)行正交展開后，可將展開結(jié)果代入塔身動(dòng)力反應(yīng)控制方程，求解控制方程進(jìn)而可以得到對(duì)應(yīng)不同θ 的動(dòng)力反應(yīng)（如塔身某點(diǎn)的速度），將對(duì)應(yīng)于不同θ 的動(dòng)力反應(yīng)代入概率密度演化方程，將求解結(jié)果對(duì)θ 積分，可得所需要的動(dòng)力反應(yīng)隨時(shí)間變化的概率密度.

風(fēng)機(jī)結(jié)構(gòu)的動(dòng)力反應(yīng)控制方程[14]為：

式中：M，C，K 分別為n × n 階質(zhì)量矩陣、阻尼矩陣與剛度矩陣，，X 為n 維加速度、速度與位移向量[13]，Q 為n 維激勵(lì)即風(fēng)荷載向量，θ 為反映激勵(lì)隨機(jī)性的隨機(jī)矩陣，本文中此隨機(jī)矩陣的行向量個(gè)數(shù)為15 個(gè)，行向量的維數(shù)遵循數(shù)論選點(diǎn)法[15，16].其概率密度函數(shù)pΘ（θ）已設(shè)定.

式中：μD為拉力系數(shù)，假定為1.2；ρ 為空氣密度；Bj為塔身單元面積.Vj為風(fēng)速（見式（13），與θ 有關(guān)）.

式中：W（θ）為用推力系數(shù)法[17，18]求得的由旋轉(zhuǎn)葉片傳遞給塔身的風(fēng)荷載推力，W（θ）=，r 為葉片長(zhǎng)度，VT為塔頂風(fēng)速，Thr為推力系數(shù).

對(duì)于疲勞損傷D，就結(jié)構(gòu)動(dòng)力學(xué)問題而言，它必連續(xù)依賴于隨機(jī)參數(shù)θ，構(gòu)造以τ 為虛擬時(shí)間參數(shù)的虛擬隨機(jī)過程Zl

顯然，D（θ，T）為Zl在τ=1 時(shí)的截口隨機(jī)變量，亦即：

（Zl，Θ）的聯(lián)合概率密度函數(shù)的概率密度演化方程為：

其初始條件為：

求解偏微分方程初值問題（18）、（19）即可給出聯(lián)合概率密度函數(shù)，進(jìn)而給出Zl（τ）的概率密度函數(shù)

由式（17）可知，

由式（21）可方便地求解隨機(jī)結(jié)構(gòu)的可靠度，亦即：

3 疲勞累積損傷理論

損傷是指在循環(huán)荷載作用下材料的損壞程度，一般用一個(gè)無量綱參數(shù)D 來表示它，當(dāng)D=0 時(shí)，說明材料完好無損，當(dāng)D ＞1 時(shí)，表示材料已經(jīng)達(dá)到它的疲勞壽命.

在隨機(jī)荷載作用下，結(jié)構(gòu)疲勞損傷分析采用疲勞累積損傷理論.目前普遍采用的理論有Palmgren[20]-Miner[21]線性疲勞損傷準(zhǔn)則.

本文主要分析風(fēng)荷載作用下風(fēng)機(jī)塔身底部混凝土基礎(chǔ)的疲勞可靠度，采用P-M 準(zhǔn)則對(duì)應(yīng)的S-N 曲線[22]為：

式中：Smax為風(fēng)荷載作用下的應(yīng)力范圍，單位為MPa，Nf為對(duì)應(yīng)Smax的導(dǎo)致材料發(fā)生疲勞破壞的循環(huán)次數(shù).疲勞損傷D 的定義為

式中：nk為第k 級(jí)應(yīng)力幅值下的實(shí)際循環(huán)次數(shù)，Nfk為第k 級(jí)應(yīng)力幅值下達(dá)到疲勞破壞時(shí)的允許循環(huán)次數(shù)，由S-N 曲線查得.k 為計(jì)算疲勞損傷時(shí)所涉及到的所有工況所對(duì)應(yīng)的應(yīng)力幅值總數(shù).

4 實(shí)例計(jì)算

本文基于某風(fēng)電場(chǎng)3 MW 風(fēng)力發(fā)電機(jī)進(jìn)行建模分析.該風(fēng)力發(fā)電機(jī)輪轂高度90 m，風(fēng)輪直徑100.8 m，額定風(fēng)速11.9 m/s，設(shè)計(jì)壽命為20 年.基礎(chǔ)混凝土采用圓形臺(tái)柱式擴(kuò)展基礎(chǔ)，底板直徑21.5 m，高3.9 m，埋深3.5 m.塔筒材料為Q345 鋼材，基礎(chǔ)環(huán)采用Q345 鋼材，基礎(chǔ)混凝土采用C35 混凝土，底部采用完全約束，采用ABAQUS 有限元軟件建立模型，選用實(shí)體單元.表1 為各段塔筒的幾何參數(shù)，采用的ABAQUS 中混凝土損傷本構(gòu)模型如圖1 所示.

表1 塔身風(fēng)荷載計(jì)算參數(shù)Tab.1 Calculation parameter of wind load of tower body

圖1 混凝土損傷本構(gòu)模型Fig.1 Concrete damage constitutive model

圖1 中，fcm為屈服強(qiáng)度，εc1為屈服應(yīng)變，Ec為彈性模量，dc為損傷因子，為非線性應(yīng)變，為塑形應(yīng)變，為彈性應(yīng)變，σc為壓應(yīng)力，εc為壓應(yīng)變，dc=取值為0.35～0.7，應(yīng)力應(yīng)變關(guān)系由混凝土規(guī)范[23]提供，將相關(guān)參數(shù)輸入軟件中計(jì)算可得所需數(shù)據(jù).

4.1 風(fēng)荷載的計(jì)算

平均風(fēng)速的選取考慮到輪轂處的工作情況，根據(jù)平均風(fēng)速的指數(shù)模型計(jì)算可選取相應(yīng)10 m 處平均風(fēng)速（標(biāo)準(zhǔn)平均風(fēng)速）12 m/s.

在計(jì)算方程（18）時(shí)，令θ=θq，θq=（θ1，q，θ2，q，…，θs，q）（s=15，q=1，2，…，Nsel），將對(duì)應(yīng)θq的D（θq，T）代入方程中，可得聯(lián)合概率密度函數(shù)，進(jìn)而對(duì)方程（20）積分并結(jié)合式（21）和（22），最終可得D值的疲勞可靠度.關(guān)于Nsel的選取遵循數(shù)論選點(diǎn)法，在Matlab 中實(shí)現(xiàn)選取步驟，對(duì)n=2 422 957 的隨機(jī)列向量采用數(shù)論選點(diǎn)法[10，15，16]進(jìn)行篩選，得到Nsel=182 個(gè)隨機(jī)列向量組，依次編號(hào)1、2、3…，方便后續(xù)整理計(jì)算.

在這182 個(gè)隨機(jī)列向量組的基礎(chǔ)上，根據(jù)算法在Matlab 中進(jìn)行編程，獨(dú)立地生成標(biāo)準(zhǔn)平均風(fēng)速為12 m/s 時(shí)的182 種風(fēng)速時(shí)間歷程，每一個(gè)隨機(jī)列向量對(duì)應(yīng)生成一個(gè)風(fēng)速時(shí)程，將脈動(dòng)風(fēng)速時(shí)程繼承隨機(jī)列向量的編號(hào)，并將脈動(dòng)風(fēng)時(shí)程與平均風(fēng)時(shí)程合并，可得到輪轂處的總風(fēng)速時(shí)程，如圖2.

圖2 標(biāo)準(zhǔn)平均風(fēng)速為12 m/s 輪轂處總風(fēng)速Fig.2 Total wind speed at hub with standard average wind speed of 12 m/s

為了驗(yàn)證風(fēng)速模擬數(shù)值的準(zhǔn)確性，在平均風(fēng)速為10.27 m/s 時(shí)進(jìn)行實(shí)測(cè)，實(shí)測(cè)風(fēng)速的采樣頻率為1/7 Hz，選擇的實(shí)測(cè)風(fēng)速為風(fēng)向穩(wěn)定且基本與應(yīng)變測(cè)點(diǎn)一致的時(shí)間段，實(shí)測(cè)結(jié)果如圖3.采用文中所述方法展開風(fēng)速為10.27 m/s 時(shí)結(jié)果如圖4.

圖3 平均風(fēng)10.27 m/s 的實(shí)測(cè)風(fēng)速Fig.3 Measured wind speed with an average wind of 10.27 m/s

由對(duì)比可知，文中風(fēng)速展開方法與實(shí)測(cè)值趨勢(shì)大體一致，可以由此確定風(fēng)速模擬取值的正確性.

圖4 平均風(fēng)速為10.27 m/s 模擬風(fēng)速Fig.4 Simulated wind speed with an average wind speed of 10.27 m/s

4.2 風(fēng)機(jī)動(dòng)力響應(yīng)的有限元模擬

將風(fēng)荷載時(shí)程加載至風(fēng)機(jī)模型上，可以得到各個(gè)隨機(jī)風(fēng)荷載下的風(fēng)機(jī)動(dòng)力響應(yīng)，本文主要觀察鋼環(huán)與混凝土接觸范圍內(nèi)的動(dòng)力響應(yīng)，取10 m 處平均風(fēng)速為12 m/s.利用公式（14）、（15）計(jì)算塔身處各點(diǎn)的風(fēng)荷載時(shí)程.

將風(fēng)荷載加至有限元模型，并考慮塔身風(fēng)荷載的影響，采用應(yīng)力等值線來表示模型內(nèi)部的應(yīng)力分布情況，可以清晰描述外動(dòng)力響應(yīng)在結(jié)構(gòu)中的分布，從而快速確定模型中的最危險(xiǎn)區(qū)域.疲勞損傷的計(jì)算是在等效應(yīng)力時(shí)程的基礎(chǔ)上計(jì)算，因此提取基礎(chǔ)環(huán)和基礎(chǔ)的等效應(yīng)力云圖.本文提取182 種隨機(jī)風(fēng)荷載時(shí)程的第一和第二個(gè)時(shí)程，編號(hào)為1、2，編號(hào)1、2 動(dòng)力響應(yīng)下基礎(chǔ)環(huán)和基礎(chǔ)等效應(yīng)力云圖如圖5所示.

圖5 基礎(chǔ)鋼環(huán)應(yīng)力云圖Fig.5 Stress of the basic steel ring

由圖5 可看出，基礎(chǔ)環(huán)在與混凝土基礎(chǔ)上部接觸部位出現(xiàn)應(yīng)力集中現(xiàn)象，初步驗(yàn)證了本文工況破壞發(fā)生的位置，進(jìn)一步提取出編號(hào)1 動(dòng)力響應(yīng)混凝土基礎(chǔ)的整體應(yīng)力云圖與剖面應(yīng)力云圖如圖6 所示.

圖6 編號(hào)1 風(fēng)荷載作用下基礎(chǔ)應(yīng)力云圖Fig.6 Stress cloud diagram of the foundation under number 1 wind load

同樣，查看其他編號(hào)動(dòng)力響應(yīng)的應(yīng)力云圖與編號(hào)1 的結(jié)果進(jìn)行對(duì)比，可看到應(yīng)力最大的位置相同，區(qū)別只是應(yīng)力大小和時(shí)程的變化.如圖7 為編號(hào)2動(dòng)力響應(yīng)的風(fēng)機(jī)基礎(chǔ)的動(dòng)力響應(yīng)結(jié)果.

圖7 編號(hào)2 風(fēng)荷載作用下基礎(chǔ)應(yīng)力云圖Fig.7 Stress cloud diagram of the foundation under number 2 wind load

由圖6 和圖7 可知，在順風(fēng)向的鋼環(huán)內(nèi)外側(cè)與混凝土基礎(chǔ)表面接觸處，動(dòng)力響應(yīng)產(chǎn)生的混凝土應(yīng)力值最大，此處為疲勞危險(xiǎn)位置，符合實(shí)際工程發(fā)生的破壞位置.

經(jīng)過對(duì)比發(fā)現(xiàn)鋼環(huán)內(nèi)側(cè)混凝土的應(yīng)力比外側(cè)更高，這意味著，當(dāng)工程實(shí)際遇到本文工況基礎(chǔ)環(huán)外側(cè)混凝土疲勞破壞時(shí)，其實(shí)鋼環(huán)內(nèi)側(cè)混凝土接觸部分更有可能已經(jīng)破壞.提取危險(xiǎn)點(diǎn)的應(yīng)力響應(yīng)，所有編號(hào)荷載作用下的危險(xiǎn)點(diǎn)應(yīng)力響應(yīng)如圖8 所示，編號(hào)1 對(duì)應(yīng)的壓應(yīng)力時(shí)程如圖9 所示，編號(hào)1 的雨流統(tǒng)計(jì)數(shù)據(jù)如圖10 所示.由于風(fēng)荷載是θ 的函數(shù)，故應(yīng)力時(shí)程平均應(yīng)力和最大幅值隨θ=182 而變化，由圖8可見，最大應(yīng)力幅值可達(dá)10 MPa.

圖8 危險(xiǎn)點(diǎn)動(dòng)力響應(yīng)結(jié)果Fig.8 Dynamic response results at dangerous points

圖9 編號(hào)1 風(fēng)速對(duì)應(yīng)的危險(xiǎn)點(diǎn)應(yīng)力時(shí)程Fig.9 Stress time history of dangerous points under number 1 wind load

圖10 編號(hào)1 應(yīng)力時(shí)程降雨計(jì)數(shù)矩陣Fig.10 Stress time history rainfall counting matrix of number 1

根據(jù)有限元模擬所得的應(yīng)力時(shí)程運(yùn)用雨流計(jì)數(shù)法進(jìn)行求解疲勞損傷，根據(jù)雨流計(jì)數(shù)法的基本思想在Matlab 中進(jìn)行編程，得到雨流計(jì)數(shù)結(jié)果之后可將結(jié)果進(jìn)行等效應(yīng)力修正，結(jié)合選取的S-N 曲線對(duì)所有編號(hào)響應(yīng)情況的應(yīng)力時(shí)程根據(jù)線性疊加理論進(jìn)行疲勞損傷計(jì)算，結(jié)果如圖11 所示.

圖11 疲勞損傷Fig.11 Fatigue damage

圖11 的疲勞損傷值為加載時(shí)長(zhǎng)600 s 時(shí)的過程所造成的疲勞損傷，將求解出的疲勞損傷轉(zhuǎn)化為以秒為單位，將之與所需計(jì)算的加載時(shí)長(zhǎng)相乘即可得到所需年限的疲勞損傷，以此可以分別計(jì)算出使用時(shí)長(zhǎng)10 年～30 年時(shí)危險(xiǎn)點(diǎn)處的疲勞損傷累計(jì)值.

將疲勞損傷定義為關(guān)于τ 的隨機(jī)過程，即虛擬隨機(jī)過程，并進(jìn)一步進(jìn)行離散，采用單邊差分法進(jìn)行計(jì)算方程（18）的數(shù)值解.

采用單邊差分方法時(shí)，其中時(shí)間步數(shù)為100，步長(zhǎng)為0.01 s，疲勞離散步數(shù)為50 步，離散區(qū)間為[0，5]，每一編號(hào)對(duì)應(yīng)的疲勞損傷結(jié)果都作為初始條件進(jìn)行一次計(jì)算，將計(jì)算所得所有編號(hào)的疲勞概率分布離散數(shù)值解相加，可以得到此年限下的概率密度數(shù)值解，不同年限時(shí)，代入不同的疲勞累計(jì)損傷即可得到該年限下的疲勞破壞概率密度數(shù)值解，全部結(jié)果展示如圖12 所示.其中疲勞破壞年限為30 年、20 年和10 年的概率密度離散值如圖13 所示.

由圖11、圖12 可看出，使用時(shí)間越短時(shí)，疲勞損傷累計(jì)值小于1 的概率分布越多，隨著年限增加疲勞概率密度曲線總體呈現(xiàn)向疲勞損傷大于1 的方向偏移的趨勢(shì).

圖12 各年限下疲勞損傷概率密度數(shù)值解Fig.12 Numerical solution of probability density of fatigue damage under various years

圖13 疲勞概率分布離散數(shù)值解之和Fig.13 Sum of the discrete numerical solutions of the fatigue probability distribution

當(dāng)疲勞損傷值在1 以下時(shí)，鋼環(huán)側(cè)混凝土不會(huì)發(fā)生疲勞破壞，換言之，結(jié)構(gòu)處于安全狀態(tài).將疲勞損傷值為1 以下的疲勞概率進(jìn)行數(shù)值積分，即為結(jié)構(gòu)處于安全狀態(tài)的概率，即為結(jié)構(gòu)的安全可靠度.可靠度計(jì)算結(jié)果如圖14 所示.

圖14 疲勞可靠度Fig.14 Fatigue reliability

由圖14 可見，隨著時(shí)間的增加，疲勞可靠度不斷降低，隨著年限的增加，疲勞可靠度的降低呈現(xiàn)不斷加快的趨勢(shì).當(dāng)年限到達(dá)30 年時(shí)，疲勞可靠度為76.66%，風(fēng)機(jī)使用時(shí)長(zhǎng)20 年時(shí)，疲勞可靠度為94.67%.從已有文獻(xiàn)[24]來看，風(fēng)機(jī)在使用壽命20 年內(nèi)，基礎(chǔ)環(huán)附近混凝土發(fā)生疲勞破壞是真實(shí)存在的.

5 結(jié)論

本文主要研究了陸上風(fēng)機(jī)在風(fēng)荷載作用下混凝土基礎(chǔ)的疲勞可靠度，具體結(jié)論如下：

1）根據(jù)隨機(jī)動(dòng)力作用的正交展開法和數(shù)論選點(diǎn)法，將風(fēng)荷載展開為分散點(diǎn)集，由此進(jìn)行動(dòng)力荷載的計(jì)算較為合理.這種方法是用概率密度演化方法求解概率密度的前提和基礎(chǔ).

2）提出了一個(gè)基于概率密度演化方法的疲勞概率計(jì)算方法，將雨流計(jì)數(shù)法得到的疲勞損傷D 值代入概率密度演化方程進(jìn)而求解基礎(chǔ)處疲勞損傷的概率密度是一種有效的方法，不但過程簡(jiǎn)單，而且精度較高.采用概率密度演化方法，可精確定量地給出結(jié)構(gòu)反應(yīng)的演化概率密度曲線族，由此，可以方便地得出結(jié)構(gòu)在隨機(jī)荷載作用下的各種可靠度.給定閥值為1，通過對(duì)概率密度在小于閥值的范圍內(nèi)進(jìn)行積分可得對(duì)應(yīng)某風(fēng)速的疲勞可靠度.

3）由表2 可以看出，本文所研究的混凝土基礎(chǔ)的疲勞可靠度隨著時(shí)間降低，當(dāng)年限到達(dá)30 年時(shí)，疲勞可靠度只有大約76.66%，這意味著，風(fēng)機(jī)基礎(chǔ)有23.34%的概率發(fā)生疲勞破壞，這是一個(gè)非常危險(xiǎn)的數(shù)值，結(jié)合本文實(shí)際工況實(shí)例中風(fēng)力發(fā)電機(jī)的設(shè)計(jì)壽命為20 年，風(fēng)機(jī)使用時(shí)間長(zhǎng)20 年時(shí)的疲勞可靠度有94.67%，意味著疲勞破壞概率為5.33%，也印證了本文實(shí)際工況的發(fā)生并不是意外情況，尤其是在風(fēng)機(jī)場(chǎng)的風(fēng)機(jī)數(shù)量基數(shù)大的情況下，發(fā)生本文工況所示的疲勞破壞的概率還是不容小覷的.