OpenFAST前处理之BModes——塔架模态分析
本文摘要（由AI生成）：
文章主要介绍了BModes软件，它是一种用于计算塔架结构动力学响应的高效模态计算软件。BModes的前身是Modes，由俄勒冈州立大学开发，后来由NREL的Marshall L. Buhl博士负责维护和改进。BModes可以计算塔架的模态振型，并生成相应的关键模态系数计算MATLAB程序。BModes在计算塔架的模态振型时，需要考虑塔架的边界条件、塔架结构属性文件、塔架材料密度、弹性模量、剪切模量、截面直径和厚度等属性。BModes还可以用于计算漂浮式风力机的塔架模态，需要考虑附加质量、静水刚度和系泊等效刚度等参数。最后，文章还介绍了如何通过多元线性回归获取模态振型拟合系数。
这几周宁波大学入职手续终于办好了，杨老师的职业生涯正式开启
上周给英国ASRANet公司培训OpenFAST，再一次熟悉了OpenFAST的基本功能和理论，特别是BModes的使用。所以，今天给大家介绍一下OpenFAST的关键前处理器——BModes。（内附关键模态系数计算MATLAB程序）
OpenFAST之所以是目前计算效率最高的风电仿真软件，其原因之一就是采用假设模态法计算塔架的结构动力学响应。塔架振动仅考虑给定的四阶模态（一阶/二阶前后及一阶/二阶侧向），从而将塔架的自由度数量降低至个位数，极大增加了结构侧计算效率。
这种情况下，需要用户提前计算塔架的模态振型。这时，我们需要用到同样高效的模态计算软件BModes。
一、BModes简介及下载
BModes的前身Modes是上世纪俄勒冈州立大学为分析直升机叶片模态开发的有限元软件。2000年后，由NREL的Marshall L. Buhl博士负责维护/改进，逐渐增加了考虑风电叶片扭角、叶尖刹车装置等功能，并进一步扩充以用于塔架模态分析，并更名为BModes。
2006年，针对大名鼎鼎的OC3项目Monopile柔性基础和Hywind Spar浮式风力机的塔架模态计算，增加了弹性边界功能，从而可以考虑土壤柔性和浮式平台对塔架模态特性的影响。不过由于负责开发这部分功能的博士后提前离职以及其他因素，因此没有发布官方版本的源代码。
之前只在JasonJonkman的NREL个人公共文件夹共享，如今随着NREL官方网站维护，共享链接已失效。为方便大家使用，我已将程序和两个算例（Monopile&Spar）上传至Github供大家下载、使用和交流。其中Monopile的算例是针对FASTv7的，其中桩与塔架统一视为Tower，这一点与OpenFAST有较大不同。
https://github.com/yang7857854/F2A/tree/main/BModes。
二、BModes输入文件介绍和运行
1. 输入文件介绍
BModes输入文件后缀为“.bmi”。以OC3Hywind.bmi为例，其输入文件可以分为八部分，每一部分以空行隔开，不能删除。
（1）常规参数设置
这部分主要需要注意4个参数：
beam_type：对于塔架需要设为2；
radius：塔架高度，即塔顶至静水面距离；
hub_rad：轮毂半径或塔架刚体部分的长度，一般应设为0；
hub_conn：底部边界条件。1：刚性基础悬臂梁；2：全柔性基础；3：除轴向和扭转外的其他自由度。这一边界条件一般选择2或3，如果选3，则忽略平台垂荡和艏摇，但是依然会计算塔架的扭转模态。
（2）塔顶质量矩阵
此处数值需针对自己的模型进行设置
（3）塔架结构属性文件
此处需定义塔架结构属性文件的路径，属性文件的格式如下：
文件第一行为注释，第二行为塔架节点数量，第三行为空行，第四、第五行为属性名及单位（实际为无功能的注释行），第六行至末尾为输入属性，行数需与节点数量保持一致。
属性数据一共有13列。第一列为归一化塔架高度（sec_loc），然后依次是几何扭角（str_tw）、惯性扭角（tw_iner）、线性密度（mass_den）、面外惯性矩（flp_iner）、面内惯性矩（edge_iner）、面外弯曲刚度（flp_stff）、面内弯曲刚度（edge_stff）、扭转刚度（tor_stff）、轴向刚度（axial_stff）、质心偏移（cg_offst）、剪切重心偏移（sc_offst）和扭转中心偏移（tc_offst）。对于塔架，示例文件中为0的变量均需设为0。
假设塔架材料密度为，弹性模量为E，剪切模量为G，截面直径和厚度分别为D和t，则对应的各个属性计算公式如下：
（4）缩放系数
无需更改，一般为1。
（5）输出坐标定义
模态振型输出坐标定义，默认的即可。
（6）平台系统定义
这部分参数都非常关键。
tow_support设为1；
draft为塔基至静水面的距离，且向下为正，塔架的真实长度为radius draft。OC3 Spar模型的塔基高于静水面10m，所以此处draft设置为-10。一般来说，海上风电基础的顶部均位于静水面以上，塔架的draft应该都是负值。
对于漂浮式风力机，cm_pform、mass_pform和i_matrix_pform根据实际情况设置即可。
ref_msl为0。hydro_M、hydro_K和mooring_K分别为附加质量、静水刚度和系泊等效刚度，参考点为（0,0,0）。其中系泊等效刚度矩阵需要采用MAP 计算得到。或者通过free-decay分析获取平台的固有周期，结合无穷频率时的附加质量以及平台质量，获取包含系泊刚度和静水刚度在内的等效刚度矩阵K，将K赋值于hydro_K，将mooring_K设为0。需要特别注意，平台惯性矩的参考点需为（0,0,0），需要对质量矩阵进行一定变换，变换方式如下：
MSWL = [(TransMat)T] Mcm[TransMat]-1
（x_cm, y_cm, z_cm）为平台重心.
对于固定式基础的海上风力机，cm_pform、mass_pform和i_matrix_pform这三个值应设为0，平台的质量和惯性矩包含在下面即将定义的“hydro_M”中。这是因为对于OpenFAST，塔架以下的部分称为“platform”。对于固定式基础的风力机，platform即是指monopile，tripod或jacket等类型的substructure。
使用OpenFAST仿真时，这一结构需要在SubDyn中建模，通过运行SubDyn可以获取platform相对于TP（连接件，也就是塔基）处的质量和刚度矩阵。将SubDyn的sum文件（xxxx.SD.sum）里面的质量和刚度矩阵分别赋值于“hydro_M”和“hydro_K”，将“mooring_K”设为0矩阵，“ref_msl”值与“draft”相同，即平台等效质量和刚度矩阵的参考点为塔基位置。
（7）其他部分
其余部分均不重要，此处不做介绍，使用默认的数值即可。
2. 运行BModes
NREL所有基于FORTRAN开发的软件与OpenFAST一样，都得在命令提示符窗口运行，BModes也没能例外。运行BModes分三步：①打开cmd窗口；②更改当前运行路径；③执行。
（1）打开cmd窗口
快捷键<Win R>，或者“开始”-“运行”，输入cmd，打开cmd窗口。
（2）更改当前运行路径
通过“cd”命令更改当前运行路径至BModes根目录。cd与目标路径之间需要用空格隔开。如果盘符不是C盘，还需要更改当前盘符，以跳至目标路径。
（3）执行
BModes执行方式与OpenFAST相同，需要在输入程序名的同时，加上文件名全称。BModes的输入文件后缀为“.bmi”。文件名与程序名中间同样用空格隔开。显示“all parameters checked: o.k.”则表示输入文件结构正确，同时运行界面还会输出各阶模态振型结果。需要特别注意的是，对于Spar模型，当我们设置边界条件为“free-free”时，平台的六个自由度的固有频率也包括在结果中了。mode 1~2对应sway和surge，mode 3为heave，mode 4~5对应roll和pitch，mode 6为yaw。mode 7及以后的模态才是塔架的模态。
3. 结果查看
BModes运行结果文件类型为.out，通过记事本打开即可查看输出的模态振型。通过上面的分析，我们获知塔架的模态从mode 7开始。通过判断s-s disp和f-a disp的数值，来判断当前模态是塔架的侧向（side-side）还是前后（fore-aft）模态。
很明显，mode7即为一阶side-side模态，mode8为一阶fore-aft模态。
但是二阶side-side模态与一阶扭转（twist）模态固有频率较为接近，可以通过比较twist与s-s disp的绝对值大小来进行判断。这一结果中，mode9的绝对值远大于s-sdisp，说明mode9主模态振型为扭转，所以判断mode9为一阶扭转模态。相应的，mode10为二阶side-side模态，mode 11为二阶fore-aft模态。
因为OpenFAST中为默认给定的塔架模态振型为归一化的振型，即塔顶位置模态位移为1.0，塔基处位移为0。我们通过对模态振型归一化后，可以看到两个方向的一阶模态基本一致，二阶模态有较为明显的差别，风轮在前后方向的偏心导致前后模态位移更大，最大值出现在75%位置附近。
可能有朋友会有点疑惑，模态归一化处理是否正确，因为求解出来的模态位移在塔基位置是不为0的，而归一化后为0。这主要是因为此时模态位移的参考系是惯性参考系，而OpenFAST中在计算时，塔架参考系为平台局部参考系，局部参考系的坐标原点就是塔基位置，意味着塔基位移一直为0，因此，塔基位移归0化处理是没有问题的。需要注意的是，OpenFAST输出的塔架运动和载荷等结果也是基于平台局部参考系的。
三、模态振型多项式系数
得到模态振型后，还需要通过多元线性回归获取模态振型拟合系数，如下：
其中，x为归一化塔架高度，y为归一化模态振型。a1~a5为需要获取的模态系数。
NREL提供了一个EXCEL文件（ModeShapePolyFitting.xls）用于获取这5个系数。该文件可以在Github上下载：
ttps://github.com/old-NWTC/FAST/tree/master/Utilities
你需要将BModes得到的各阶模态振型数据，复 制到“Input”表单的“x,y”处，并修改右侧的“bottom slope”数值（灰色填充部分），即BModes结果中对应的Slope值，或者采用“Estimated slope atbottom”处的值。然后在“6thOrder Polynomial”表单中查看归一化的五个系数。
一般来说，Improved Direction和Projection方法的结果比较接近，也更为准确。可根据右侧的拟合结果选择合适的5个模态系数，用于替换ElastDyn模块中塔架属性文件中的模态系数。
用户只需要重复以上步骤5次，就可以获取4个模态所需的20个系数。
个人觉得这有点麻烦，复 制粘贴来来去去好几回，容易出错。所以，我写了一个通用的MATLAB程序，通过给定BModes结果文件和各阶模态对应的序号，调用“ReadBModOut”、“getModeShapeCoef”和“writeCoeffs”即可直接得到OpenFAST所需的20个系数。各个子程序分享在GitHub平台：
https://github.com/yang7857854/F2A/tree/main/BModes
四、柔性基础效应
据Jason Jonkman在NWTC论坛上“坦白”，所有的算例均没有考虑基础柔性对塔架模态的影响。所以默认算例都是基于刚性基础。
通过BModes重新计算了塔架的模态系数，发现对于Monopile模型，刚性和柔性基础结果均几乎完全一致。说明考虑30m长monopile对塔架模态的影响与否，不影响塔架动力响应的计算结果。这主要是因为该模型的monopile仅有30m长，且泥面处是刚性假设，基础对塔架提供的刚度足够大（平动刚度为9次方，转动刚度为12次方）。所以对于塔架的模态分析，弹性边界假设与刚性边界假设的结果较为接近。
对于Spar模型也得到了类似的结果，尽管平台对模态振型影响较大，但是由于漂浮式风力机在平动方向的刚度主要由系泊系统提供，而悬链线式的系泊提供的平动刚度较小，在风载荷作用下，塔架-平台系统运动主要为平台的平动方向，塔架的弯曲模态本身贡献就很小，因平台导致的模态差异引起的响应差别就更小了。所以，考虑平台刚度对塔架模态的影响与否，对于塔架动力学响应的计算似乎并不重要。
五、结语
尽管最终发现海上风电基础对于塔架模态的影响，不会进一步影响塔架的动力学响应结果。但是上述流程可以增进大家对于OpenFAST中塔架结构动力学处理方法的理解。特别是在自己更换了塔架模型后，如何获取OpenFAST所需的20个模态系数，不再是一个麻烦的事情。
