基于影响矩阵的风电机组螺栓疲劳寿命分析——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月毕业于华中科技大学机械工程学院机电系，现在东方电气风电机有限公司从事结构分析工作。