張俊

（東方電氣風(fēng)電有限公司，四川 德陽，618000）

基于影響矩陣的風(fēng)電機(jī)組螺栓疲勞壽命分析
——Excel VBA開發(fā)

張俊

（東方電氣風(fēng)電有限公司，四川 德陽，618000）

通常情況下，風(fēng)電機(jī)組螺栓應(yīng)力對外載呈現(xiàn)出非線性關(guān)系，即應(yīng)力不隨外載線性變化。而現(xiàn)有的疲勞分析軟件（如FE-Safe、nCode等）應(yīng)用的前提條件都是應(yīng)力隨外載線性變化，故現(xiàn)有的疲勞分析軟件對風(fēng)電機(jī)組螺栓疲勞壽命的計(jì)算并不適用。文章介紹基于影響矩陣的風(fēng)電機(jī)組螺栓疲勞壽命分析步驟，以輪轂與主軸連接螺栓為例并結(jié)合Excel VBA二次開發(fā)程序詳細(xì)介紹各個步驟的具體實(shí)現(xiàn)方式，最后簡單介紹該二次開發(fā)程序中用到的MatrixVB插件及雨流計(jì)數(shù)算法。

風(fēng)電機(jī)組，輪轂與主軸連接螺栓，疲勞，影響矩陣，Excel VBA，MatrixVB，雨流計(jì)數(shù)算法

0　引言

通常情況下，風(fēng)電機(jī)組螺栓應(yīng)力對外載呈現(xiàn)出非線性關(guān)系，即應(yīng)力不隨外載線性變化，而現(xiàn)有的疲勞分析軟件 （如FE-Safe、nCode等）都是在應(yīng)力隨外載線性變化的假設(shè)下，通過外載的時間序列乘以單位外載作用于結(jié)構(gòu)的有限元應(yīng)力結(jié)果得出結(jié)構(gòu)的應(yīng)力時間序列，然后根據(jù)累計(jì)疲勞損傷準(zhǔn)則 （Miner準(zhǔn)則）計(jì)算結(jié)構(gòu)的疲勞壽命。因此現(xiàn)有的疲勞分析軟件對風(fēng)電機(jī)組螺栓疲勞壽命的計(jì)算并不適用。

目前主要采用曲線擬合法或影響矩陣法計(jì)算風(fēng)電機(jī)組螺栓的疲勞壽命。曲線擬合法通過有限元分析結(jié)果擬合出各個載荷分量的大小對螺栓應(yīng)力的關(guān)系曲線，再結(jié)合各個載荷分量大小的時間序列算出螺栓的各個應(yīng)力時間序列，最后合并各個應(yīng)力時間序列得出螺栓總的應(yīng)力時間序列，進(jìn)而統(tǒng)計(jì)螺栓的疲勞損傷；影響矩陣法通過有限元分析結(jié)果提取合成彎矩的大小和方向?qū)β菟☉?yīng)力（1個軸向，2個彎曲）的3個二維影響矩陣，再將合成彎矩大小和方向的時間序列 （由Bladed后處理生成）在此3個影響矩陣的基礎(chǔ)上分別進(jìn)行二維插值，得出螺栓的3個應(yīng)力時間序列 （1個軸向，2個彎曲），然后通過三角插值公式計(jì)算螺栓危險截面上各個點(diǎn)的總的應(yīng)力時間序列，進(jìn)而統(tǒng)計(jì)螺栓的疲勞損傷。因影響矩陣法同時考慮了載荷的大小和方向?qū)β菟☉?yīng)力的影響，故較曲線擬合法更精確。但影響矩陣法涉及到數(shù)據(jù)的二維插值操作，而Bladed目前版本的后處理中只能實(shí)現(xiàn)數(shù)據(jù)的一維插值操作，故采用此方法分析風(fēng)電機(jī)組螺栓的疲勞壽命時，有必要編寫二次開發(fā)程序以實(shí)現(xiàn)分析數(shù)據(jù)的快速處理。

1　風(fēng)電機(jī)組螺栓疲勞壽命分析步驟

風(fēng)電機(jī)組螺栓疲勞壽命分析步驟如下：

（1）從各個有限元工況的分析結(jié)果中，針對每個螺栓，提取合成彎矩大小和方向?qū)β菟ㄎｋU截面 （靠近螺紋旋合部位的截面）應(yīng)力的3個影響矩陣 （螺栓采用梁單元模擬，可提取x向的軸向應(yīng)力及y、z向的彎曲應(yīng)力）；

（2）將合成彎矩大小和方向的時間序列在3個影響矩陣的基礎(chǔ)上分別進(jìn)行二維插值，得到每個螺栓危險截面上的3個應(yīng)力時間序列；

（3）依據(jù)3個應(yīng)力時間序列，采用三角插值公式計(jì)算每個螺栓危險截面圓上每隔30°的點(diǎn) （12個）的應(yīng)力時間序列；

（4）針對每個螺栓，對以上12個點(diǎn)的應(yīng)力時間序列依次進(jìn)行雨流計(jì)數(shù)，記錄各個循環(huán)的應(yīng)力范圍，根據(jù)GL標(biāo)準(zhǔn)［1］計(jì)算螺栓的疲勞等級，根據(jù)Eurocode 3-1-9［2］選取SN曲線，依次統(tǒng)計(jì)各個點(diǎn)的年損傷值，將12個點(diǎn)中的最大年損傷值定為該螺栓的年損傷值。

下面以輪轂與主軸連接螺栓為例并結(jié)合Excel VBA二次開發(fā)程序詳細(xì)介紹以上各個步驟的具體實(shí)現(xiàn)方式，輪轂與主軸螺栓連接有限元模型如圖1所示。

圖1　輪轂與主軸螺栓連接有限元模型

1.1提取影響矩陣

在有限元模型中，首先分析求解螺栓預(yù)緊力工況，然后將合成彎矩分為12個方向，每個方向進(jìn)行一次分析求解，求解時從第二載荷步開始（第一載荷步為已分析求解的螺栓預(yù)緊力工況）施加合成彎矩的大小，分4個子步進(jìn)行施加，故需求解的有限元工況數(shù)為 1+12×4=49個。通過APDL宏 （即命令流）從各個有限元工況的分析結(jié)果中提取合成彎矩的大小和方向?qū)β菟ㄎｋU截面應(yīng)力的3個影響矩陣 （見表 1）到txt文件，每個螺栓對應(yīng)一個txt文件。

基于Excel VBA的風(fēng)電機(jī)組螺栓疲勞壽命分析二次開發(fā)程序界面如圖 2所示，該程序的運(yùn)行涉及到很多的矩陣運(yùn)算，因此在編寫VBA代碼時，引用了MatrixVB插件。單擊界面上的按鈕RetrieveData，彈出打開文件對話框，選擇保存影響矩陣的所有txt文件，打開后各個螺栓的3個影響矩陣將被導(dǎo)入 Excel的各個對應(yīng)的工作頁（“bolt1”，“bolt2”，...“bolt60”）中。

表1　合成彎矩的大小和方向?qū)蝹€螺栓危險截面應(yīng)力的二維影響矩陣

圖2　Excel VBA二次開發(fā)程序界面（導(dǎo)入影響矩陣到Excel中）

1.2生成螺栓危險截面應(yīng)力時間序列

首先在二次開發(fā)程序界面中單擊如圖3所示的按鈕CopyLTSInOneDir，將目錄D：LTS_M_YZ下的所有子目錄中后綴名為S101的所有載荷時間序列文件拷貝到目錄D：LTS_M_YZall下，方便程序讀取載荷時間序列文件，計(jì)算前，必須保證疲勞工況 （Load Case）、頻次 （Frequency）及指定螺栓編號 （Selected Bolts）和疲勞等級 （DC）均已正確輸入，因?yàn)槌绦蛟谟?jì)算時將用到這些信息。然后單擊按鈕Go，程序開始計(jì)算 （如圖4所示），計(jì)算過程中通過二維插值生成各個指定螺栓危險截面的3個應(yīng)力時間序列（σaxial（t），σbending_1（t）和σbending_2（t）），并將其保存在內(nèi)存中。

圖3　Excel VBA二次開發(fā)程序界面（拷貝所有載荷時間序列文件到一個目錄下）

圖4　Excel VBA二次開發(fā)程序界面（計(jì)算各個指定螺栓的年損傷值）

1.3生成螺栓危險截面圓上各點(diǎn)的應(yīng)力時間序列

輪轂與主軸連接螺栓不僅承受軸向拉伸載荷，還承受彎曲載荷，故有必要對螺栓應(yīng)力截面圓上的多個點(diǎn)進(jìn)行疲勞計(jì)算，如圖5所示，每隔30°的點(diǎn)的應(yīng)力時間序列如下：

其中，β=0°，30°，…，330°。

程序在計(jì)算過程中 （如圖4所示），將自動計(jì)算這些點(diǎn)的應(yīng)力時間序列并將其保存在內(nèi)存中。

圖5　螺栓危險截面圓上點(diǎn)的定義（β=0°，30°，...，330°），螺栓連接坐標(biāo)系

1.4計(jì)算螺栓年損傷值

程序計(jì)算完成后，將自動輸出各個指定螺栓的年損傷值、最大年損傷值點(diǎn)的位置β、壽命和運(yùn)行20年的應(yīng)力儲備系數(shù) （SRF）。若所有指定螺栓的SRF都大于1，則螺栓的疲勞壽命滿足要求（如圖6所示）。

2　MatrixVB插件介紹

MatrixVB是由 MATHWORKS公司提供的COM組件，包含了大量與MATLAB相似的函數(shù)與調(diào)用語法，可以加強(qiáng)VB的數(shù)學(xué)運(yùn)算與圖形顯示功能，在VB程序代碼中可以像使用VB自己的函數(shù)一樣使用MatrixVB的函數(shù)，從而輕松地在Visual Basic中完成矩陣運(yùn)算與圖形繪制及顯示等功能。

MatrixVB插件安裝好后，在VBA編輯器中，單擊菜單上的 “工具”—＞“引用”，然后選擇“MMatrix”，即完成了使用MatrixVB的準(zhǔn)備工作。

調(diào)用MatrixVB中的函數(shù)時，可直接將Excel中的Range對象作為函數(shù)的參數(shù)，如下面的VBA代碼 （Excel當(dāng)前工作表中A1到A5單元格的數(shù)值分別為1，0，2，0，3，代碼中有單引號的行為注釋文本）：

′獲取當(dāng)前工作表中代表A1到A5單元格的Range對象的引用

Set Rng=Range（"A1：A5"）

′調(diào)用 MatrixVB中的函數(shù) nonzeros，將變量Rng引用的對象作為該函數(shù)的參數(shù)

result=nonzeros（Rng）

MatrixVB中函數(shù)的返回值類型一般為Matrix對象，以上的result即為一個引用Matrix對象的變量：

圖6　各個指定螺栓運(yùn)行20年的應(yīng)力儲備系數(shù) （SRF）

若要將result引用的Matrix對象中的數(shù)據(jù)輸入到Excel當(dāng)前工作表B1到B3單元格里，則必須使用Matrix對象中的Simple方法將Matrix對象轉(zhuǎn)換為VBA中的數(shù)組，如下的VBA代碼即可實(shí)現(xiàn)這一操作：

Range（"B1：B3"）.Value=result.Simple

3　雨流計(jì)數(shù)算法介紹

雨流計(jì)數(shù)法又名 “塔頂法”，由Matsuishi和T.Endo提出。雨流計(jì)數(shù)法在疲勞壽命計(jì)算中應(yīng)用非常廣泛，用來精確統(tǒng)計(jì)各個應(yīng)力或應(yīng)變區(qū)域（區(qū)域的大小由劃分的bin數(shù)確定，bin數(shù)越多，區(qū)域越小，統(tǒng)計(jì)結(jié)果越精確）的循環(huán)次數(shù)。把應(yīng)力或應(yīng)變-時間歷程曲線圖 （見圖7）順時針轉(zhuǎn)90°，使時間坐標(biāo)軸豎直向下，曲線猶如一系列屋頂，雨水順著屋頂往下流，故稱為雨流計(jì)數(shù)法。

圖7　時間歷程曲線示意圖

本文提到的二次開發(fā)程序中使用的雨流計(jì)數(shù)算法步驟如下：

（1）根據(jù)原始的應(yīng)力或應(yīng)變時間序列提取波峰波谷序列；

（2）為了整個計(jì)數(shù)過程中不出現(xiàn)殘余的半循環(huán)，將波峰波谷序列循環(huán)移位，使序列中絕對值最大的點(diǎn)位于序列的首位，如圖7所示；

（3）如圖8所示的流程圖中，dSC表示循環(huán)移位后的波峰波谷序列，dBuf表示為進(jìn)行雨流計(jì)數(shù)而定義的緩沖區(qū) （即VBA數(shù)組），dBuf（Index）、dBuf（Index-1）、dBuf（Index-2）分別存放圖9中的A′、B′、C′處的應(yīng)力或應(yīng)變值，dRngArr表示計(jì)數(shù)過程中記錄應(yīng)力或應(yīng)變范圍的動態(tài)數(shù)組，dNPnt表示dSC中的數(shù)據(jù)個數(shù)。

圖8　雨流計(jì)數(shù)算法流程圖

圖9　雨流計(jì)數(shù)過程示意圖

如下VBA函數(shù)的功能對應(yīng)以上雨流計(jì)數(shù)算法步驟的 （2）和 （3）。

Function RainFlowCount（ByVal targeCol As_

Variant）As Matrix

'Convert targeCol to matrix and retrieve the

'abs.max.value and its index as matrix.

mTC=plus（targeCol，0）

mMaxP=mmax（mTC）

mMinP=mmin（mTC）

mMaxPnt=plus（times（ge（mabs（mMaxP），_

mabs（mMinP）），mMaxP），times（lt（mabs（mMaxP），_

mabs（mMinP）），mMinP））

mIndexMP=findstr（mMaxPnt，mTC）

dShift=minus（mIndexMP，1）.Simple

'Circularly shift the max.value in"mTC"to the top.

mSC=RowShiftUp（mTC，dShift）

'Append the abs.max.value"mMaxPnt"to the

'bottom of"mSC".

mSC=vertcat（mSC，mMaxPnt）

'Get number of points in"mSC".

dNPnt=Length（mSC）.Simple

'Use double array dSC（）to store mSC.

dSC=mSC.Simple

′Define buffer dBuf（）for rainflow count.

Dim dBuf（1 To 8 192）As Double

′Define dynamic double array dRngArr（）to store

′rainflow counting results.

Dim dRngArr（）As Double

′RAINFLOW COUNT.

Index=0

k=0

For i=1 To dNPnt

Index=Index+1

dBuf（Index）=dSC（i，1）

Do While Index＞2

If Abs（dBuf（Index-1）-dBuf（Index-2））＜=_

Abs（dBuf（Index）-dBuf（Index-1））Then

dRng=Abs（dBuf（Index-1）-dBuf（Index-2））

Index=Index-2

dBuf（Index）=dBuf（Index+2）

′Record the range.

k=k+1

ReDim Preserve dRngArr（1 To 1，1 To k）

dRngArr（1，k）=dRng

Else

Exit Do

End If

Loop

Next

Set RainFlowCount=Transpose（plus（dRngArr，0））

End Function

4　總結(jié)

本文介紹了基于影響矩陣的風(fēng)電機(jī)組螺栓疲勞壽命分析步驟，以輪轂與主軸連接螺栓為例并結(jié)合Excel VBA二次開發(fā)程序詳細(xì)介紹了各個步驟的具體實(shí)現(xiàn)方式，最后簡單介紹了該二次開發(fā)程序中用到的MatrixVB插件及雨流計(jì)數(shù)算法。本文介紹的影響矩陣法同時考慮了載荷的大小和方向?qū)β菟☉?yīng)力的影響，較曲線擬合法更精確。另外，利用Excel VBA二次開發(fā)的螺栓疲勞壽命計(jì)算程序操作起來比較方便和靈活，具有較高的實(shí)用價值。

［1］Germanischer Lloyd.Guideline for the Certification of Wind Turbines［S］，2010

［2］EN 1993-1-9，Eurocode 3：Design of steel structures-part 1-9：Fatigue，January 2006

［3］VDI2230 Part 1，Systematic Calculation of High Duty Bolted Joints，Joints with One Cylindrical Bolt［S］

［4］MatrixVB，MatrixVB User's Guide［DB］，June 2000

［5］ANSYS，ANSYS Workbench 14.0 Help Documentation［DB］，Mechanical APDL ANSYS Parametric Design Language Guide

［6］Steve Saunders，Jeff Webb，Programming Excel with VBA and.NET［M］，O'Reilly，2006

Influence Matrix-based Fatigue Life Analysis of Wind Turbine Bolts—Excel VBA Development

Zhang Jun

（Dongfang Electric Wind Power Co.，Ltd.，Deyang Sichuan，618000）

The relationship between the stress of wind turbine bolts and the external load is usually nonlinear，in other words，the stress doesn't change linearly with the external load.However，fatigue analysis softwares（e.g.FE-Safe，nCode，etc.）available now can only be used if the stress changes linearly with the external load，thus they're not suitable for the fatigue life analysis of wind turbine bolts.In this article，procedures of the influence Matrix-based fatigue life analysis of wind turbine bolts are introduced.Besides，an Excel VBA program is developed for analysis procedures，and the hub-main shaft bolted connection is taken as an example to describe in detail program operations for various procedures.Finally，the MatrixVB add-in and rainflow counting algorithm used in the Excel VBA program are briefly described.

wind turbine，hub-main shaft bolted connection，fatigue，influence matrix，Excel VBA，MatrixVB，rainflow counting algorithm

TK83

B

1674-9987（2015）02-0035-06

10.13808/j.cnki.issn1674-9987.2015.02.007

張俊 （1983-），男，工學(xué)碩士，2007年3月畢業(yè)于華中科技大學(xué)機(jī)械工程學(xué)院機(jī)電系，現(xiàn)在東方電氣風(fēng)電機(jī)有限公司從事結(jié)構(gòu)分析工作。