基于影響矩陣的風電機組螺栓疲勞壽命分析
——Excel VBA開發

張俊

（東方電氣風電有限公司，四川 德陽，618000）

基于影響矩陣的風電機組螺栓疲勞壽命分析
——Excel VBA開發

張俊

（東方電氣風電有限公司，四川 德陽，618000）

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

風電機組，輪轂與主軸連接螺栓，疲勞，影響矩陣，Excel VBA，MatrixVB，雨流計數算法

0　引言

通常情況下，風電機組螺栓應力對外載呈現出非線性關系，即應力不隨外載線性變化，而現有的疲勞分析軟件 （如FE-Safe、nCode等）都是在應力隨外載線性變化的假設下，通過外載的時間序列乘以單位外載作用于結構的有限元應力結果得出結構的應力時間序列，然后根據累計疲勞損傷準則 （Miner準則）計算結構的疲勞壽命。因此現有的疲勞分析軟件對風電機組螺栓疲勞壽命的計算并不適用。

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

1　風電機組螺栓疲勞壽命分析步驟

風電機組螺栓疲勞壽命分析步驟如下：

（1）從各個有限元工況的分析結果中，針對每個螺栓，提取合成彎矩大小和方向對螺栓危險截面 （靠近螺紋旋合部位的截面）應力的3個影響矩陣 （螺栓采用梁單元模擬，可提取x向的軸向應力及y、z向的彎曲應力）；

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

（3）依據3個應力時間序列，采用三角插值公式計算每個螺栓危險截面圓上每隔30°的點 （12個）的應力時間序列；

（4）針對每個螺栓，對以上12個點的應力時間序列依次進行雨流計數，記錄各個循環的應力范圍，根據GL標準［1］計算螺栓的疲勞等級，根據Eurocode 3-1-9［2］選取SN曲線，依次統計各個點的年損傷值，將12個點中的最大年損傷值定為該螺栓的年損傷值。

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

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

1.1提取影響矩陣

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

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

表1　合成彎矩的大小和方向對單個螺栓危險截面應力的二維影響矩陣

圖2　Excel VBA二次開發程序界面（導入影響矩陣到Excel中）

1.2生成螺栓危險截面應力時間序列

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

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

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

1.3生成螺栓危險截面圓上各點的應力時間序列

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

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

程序在計算過程中 （如圖4所示），將自動計算這些點的應力時間序列并將其保存在內存中。

圖5　螺栓危險截面圓上點的定義（β=0°，30°，...，330°），螺栓連接坐標系

1.4計算螺栓年損傷值

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

2　MatrixVB插件介紹

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

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

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

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

Set Rng=Range（"A1：A5"）

′調用 MatrixVB中的函數 nonzeros，將變量Rng引用的對象作為該函數的參數

result=nonzeros（Rng）

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

圖6　各個指定螺栓運行20年的應力儲備系數 （SRF）

若要將result引用的Matrix對象中的數據輸入到Excel當前工作表B1到B3單元格里，則必須使用Matrix對象中的Simple方法將Matrix對象轉換為VBA中的數組，如下的VBA代碼即可實現這一操作：

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

3　雨流計數算法介紹

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

圖7　時間歷程曲線示意圖

本文提到的二次開發程序中使用的雨流計數算法步驟如下：

（1）根據原始的應力或應變時間序列提取波峰波谷序列；

（2）為了整個計數過程中不出現殘余的半循環，將波峰波谷序列循環移位，使序列中絕對值最大的點位于序列的首位，如圖7所示；

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

圖8　雨流計數算法流程圖

圖9　雨流計數過程示意圖

如下VBA函數的功能對應以上雨流計數算法步驟的 （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　總結

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

［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-），男，工學碩士，2007年3月畢業于華中科技大學機械工程學院機電系，現在東方電氣風電機有限公司從事結構分析工作。