图解AQWA和OpenFAST建立漂浮式风力机全耦合模型程序
导读：8月26日20时，我们将邀请宁波大学副教授杨阳老师，通过仿真秀平台与风电专业相关的研究生、高年级本科生、从事风电行业工作人员和OpenFAST用户开启一堂《漂浮式风力机和风-流集成能源装备仿真技术》线上讲座，敬请关注。
在正式开播前，我们分享一篇杨阳老师原创的迄今为止最满意的工作：FAST2AQWA（F2A）的开发历程给大家，希望能够对学习工程师和理工科学子有所帮助。该文首发杨阳老师的浮草之阳《在AQWA中建立漂浮式风力机全耦合模型-图文详解及程序分享》，略有改动，如果不当，欢迎批评指正。
以下是正文：
一、写在前面
之所以对这项工作最满意，一是解决了AQWA用户的一大技术难点——无法在AQWA中考虑漂浮式风力机的气动伺服弹性与平台运动之间的耦合作用。据我所知，很多人都想实现这项技术，但一直没有人做出来。二是给很多使用AQWA的研究生们提供了一种新的研究工具。自己的工作能给别人带来帮助是最欣慰的事情，也是莫大的肯定
首先说一下，F2A已经在GitHub开源，源代码和算例说明以及基本的用户手册也可以在GitHub平台下载。先简单介绍一下FAST和AQWA这两个软件。
二、风力机一体化仿真分析FAST
FAST是美国可再生能源实验室（NREL）开发的水平轴风力机气动-水动-伺服-弹性仿真计算软件，其发展历程至今有20多年。
初代版本FAST v1.01由Wilson和Walker等人在1996年发布；
1999年，结合20世纪60年代就已经开发的AeroDyn模块，形成了FAST的前身FAST_AD。
2000年，FAST灵魂人物Jason Jonkman正式加入NREL，开始与Marshall Buhl共同维护FAST。
2007年，Jonkman在其博士论文中开发了水动力计算和准静态系泊模块，主导推出了FAST v6，奠定了如今的FAST的核心基础。
2010年，发布了FAST v7。
2011年，NREL风能研究团队决定基于FAST v7版本开展模块化开发，即根据风力机不同部件的结构和载荷特性，开发相对应的模块，将FAST分为AeroDyn（气动）、ElastDyn（弹性）、ServoDyn（伺服）、SubDyn（基础）、HydroDyn（水动）和MoorDyn（系泊）等几大模块。 2016年，发布FAST v8.16之后，在GitHub上启动了OpenFAST的开发与开源交流，更为高效的开发新功能并修复其中的bug。
OpenFAST目前的版本具有很强大的风力机一体化仿真功能，对于超长叶片的大变形问题，采用自由涡尾迹和几何精确梁方法计算其气动弹性响应，在气动和结构方面均具有较高的精度和计算效率。伺服控制方面，可以采用TMD和TCLD进行结构被动控制。
OpenFASTv3.0支持漂浮式平台结构建模，通过SubDyn考虑平台结构响应。HydroDyn中也支持多个浮体的水动力耦合。因为OpenFAST具有太多模块，代码采用很多复杂的派生结构，一时之间难以理解。因此，我选择了功能同样齐全的FAST v7作为F2A的开发基础。在F2A中，FAST的作用相当于就是求解一个陆上风力机的动力学响应，针对大部分情况，FAST v7依然适用。
三、浮体水动力分析AQWA
AQWA是上世纪70年代由Atkins公司开发的海洋工程结构水动力分析软件，可以用于船舶及浮式平台水动力计算和系泊动力学分析，基于势流理论计算静水刚度、附加质量、辐射阻尼和波浪力，采用刚体运动假设在时域中计算平台运动响应。后面ANSYS Inc收购了AQWA的运营方Century Dynamic，将其融入了ANSYS软件包。由于ANSYS在国内相关行业的统治地位，AQWA受众较广，是大部分研究生用于浮体水动力分析的第一选择。
关于AQWA更多功能介绍和使用方法，可以参考海工领域著名工程师、仿真秀专栏作者高巍主编的相关内容。
由于AQWA的应用对象主要是船舶及油气平台，一般没有复杂的气动载荷。对于漂浮式风力机的水动力分析，目前的版本不能准确计算在平台大幅运动影响下的风轮非定常气动力以及叶片/塔架等柔性结构的弹性变形。但AQWA提供了一个用户自定义载荷接口，可通过FORTRAN语言编写的动态链接库（user_force64.dll）或者Python编写的外部载荷服务器计算平台受到的外部载荷。
四、F2A开发原因和基本思路
（1）开发原因
之所以想着去开发F2A，就是为了解决一种多浮体平台的气液伺服弹问题。在英国做博士后的时候，参与的欧盟资助项目是关于10 MW漂浮式风电技术的应用和实现（ARCWIND）。在ARCWIND项目中，提出了三种漂浮式风电平台模型，分别是西班牙Esteyco公司的TELWIND模型、Saitec公司的SATH模型和葡萄牙里斯本大学Carlos Guedes Soares教授团队提出的一种四脚张力腿模型CENTEC-TLP。
这几种模型都在爱尔兰University of Cork进行了波浪水池实验，我的任务其实是根据合作伙伴西班牙可再生能源实验室提供的动态响应数据，提出并验证结构健康监测系统设计方法和故障诊断模型算法等。但是因为种种原因，只能自己想办法计算动态响应数据。所以，针对这几种模型，需要采用气动-水动-伺服-弹性耦合分析工具计算风-浪-流联合作用下风力机各个部件的结构动力学响应。
但TELWIND模型比较特殊，是由上浮体和下浮体组成的多体平台模型，中间由高强度筋腱连接。目前主流软件（比如HAWC2、FAST）并不支持这一类型的平台建模。而AQWA可以用于多个浮体分析，筋腱可以采用系泊建模，所以就想到将FAST的气动-伺服-弹性仿真功能集成到AQWA的user_force64.dll中，结合AQWA的水动力计算和系泊动力学仿真，如此，就可以在AQWA中计算漂浮式风力机的全耦合动态响应了。
（2）基本思路
F2A本质上是AQWA用于外部载荷计算的DLL，在这个DLL中，融合了FAST气动-伺服-弹性仿真，并与AQWA进行数据交换。对于任一漂浮式风力机，风轮、机舱和塔架等结构的载荷和运动响应都在DLL中求解，而平台运动响应以及水动载荷、系泊恢复力则在AQWA里面直接求解。
基本思路是采用AQWA中计算得到的平台运动响应覆盖FAST中对应数值，以计算各结构部件动力响应和相应载荷，再将得到的载荷传递至AQWA，以求解平台的六自由度响应。
具体实现时，需要注意以下几点：
（1）载荷及运动响应坐标变换FAST/DLL计算得到的塔基载荷，经过坐标变换后传递至AQWA中，作用于平台的重心，并用于计算平台的运动响应。在下一时间步，将平台运动响应进行一定的坐标变换，将重心处的响应变换成平台参考点（一般为坐标原点，具体在FAST中设置）的运动响应，传递至FAST/DLL中，用于气动力和结构动力响应计算。平台的速度和加速度都会影响塔架、机舱和风轮的运动状态，因此这一步的坐标变换至关重要。
在AQWA中，采用欧拉角方法进行坐标变化。所以平台重心到平台参考点的运动响应变换采用如下公式进行。需要注意的是，加速度是AQWA的解，当前时间步的加速度是未知的，而FAST/DLL又需要这一数据，所以根据速度数据采用简单的向上差分进行计算，即当前时间步的速度减去上一时间步的速度，再除以时间步长。
在FAST/DLL中完成一个时间步的计算后，将会得到塔架的载荷。因为平台与塔基距离最近，因此，采用塔基载荷作为载荷变换的基础。塔基载荷的参考坐标系是坐落于平台与塔架连接点的局部坐标系，垂向力Fz的方向为平台重心指向塔基的方向。且时刻与平台运动保持一致。因为载荷变换的目的是将塔基载荷转变为作用于平台重心CM处的载荷，参考系由局部参考系转为全局参考系。具体方法如下：
（2）平台运动赋值
在代码实现时，另一个特别需要注意的地方就是FAST程序中对平台运动赋值。由于FAST在进行塔架、机舱和风轮的运动学计算时，只有打开了平台自由度，才会考虑平台运动响应对RNA（Rotor-Nacelle Assembly）和塔架的影响。但是平台的自由度是在AQWA中计算并传递到FAST中的，在FAST程序中并不计算平台的自由度。所以，在对平台运动赋值时，还需要修改FAST中求解equation of motion的相关代码。FAST中，equation of motion的线性方程组是通过Gauss消元法求解的，需要在输入矩阵中，剔除与平台六个自由度相对应的行和列。
针对FAST源代码的修改部分主要包括：
1）将FAST主程序的类型由“Program”改为“Subroutine”，供user_force程序调用，以进行数据交换。
2）修改“TimeMarch”子程序，约束时间步进的次数。因为AQWA似乎是采用Newmark-Beta方法进行时间推进，预测步与校正步时间步长均为dt，且均需要塔基载荷输入，因此，在预测和校正阶段均需要调用user_force64.dll。如果每次调用FAST程序计算的时间长度也为dt，那么FAST中的动态数据就会紊乱。为了数值可行性，每次调用FAST时只推进0.5dt，这样预测校正两个阶段共使得FAST中的数据推进了dt，与AQWA保持一致。一般地，AQWA中时间步长设为0.1s，FAST中时间步长设为0.005s，意味着预测和校正步均需调用FAST程序执行10次循环。这里通过修改“TimeMarch”程序实现。经过测试，AQWA中的时间步并不影响最终结果。
3）修改“Gauss”子程序，从而忽略平台自由度求解。
4）修改“Solver”子程序，在其中对平台运动响应进行赋值。
5）在“CalcOuts”程序中获取塔基载荷，用于载荷变换和传递。
6）“SAVE”所有用于运动学和动力学计算的变量，因为此时FAST程序为一个外部DLL，所有未保存的数据均会被视为首次使用，会采用初始化数值进行计算。
关于FAST的其他修改部分可以在GitHub上下载的源代码查看，所有修改前一行均有以“! Yang”开头的注释。
五、F2A运行与验证
AQWA用户应该都已经明白了，F2A本质上就是AQWA用于外部载荷计算的DLL。因此，使用F2A计算漂浮式风力机动力学响应，需要做一定的前处理，并在AQWA中执行仿真。去年12月13号我做了一个公益讲座，但是参与的人不是很多，只有90人左右。当时详细介绍了如何开发F2A和如何使用F2A用于仿真，相关资料和录屏我下一期发布（再水一期）。
我个人习惯使用cmd窗口执行AQWA仿真命令，如下：
注：去年11月份的时候，采用欧拉角方法代替了FAST中的近似小角度变换矩阵 以湍流风 静水工况和OC4半潜式平台为例，比较了F2A与OpenFAST的计算结果。叶尖位移、风轮转速和风轮推力基本上完全一致，说明FAST中的气动-伺服-弹性计算功能完美的融合进了F2A中。
平台平动位移也非常吻合，垂荡方向的差距稍大一些。OpenFAST中用到的平台频域结果与F2A中完全一致，所以需要对OpenFAST自带的OC4算例进行一点修改。主要是平台的排水体积，这一数值对浮力影响较大，如果HydroDyn中定义的排水体积与AQWA中不同，那么F2A与OpenFAST计算得到的heave的平均值会有所不同。
平台转动位移除了roll方向差别比较大之外，平台较为重要的pitch和yaw差别都很小。一定程度上讲，F2A还是很成功的。
六、F2A开发历程
实际上，在着手F2A开发前，还曾将AeroDyn直接用于AQWA，也是通过user_force实现的，构建了一个伺服控制部分，忽略了叶片和塔架的弹性效应。但是耦合效果不好，所以后来在疫情期间才开始着手F2A的开发工作。
针对FAST源代码的修改，之前断断续续思考过小半年，但是一直没有动手。直到2020年3月10号才正式开始，3月15号调试完成，可以得到很准确的结果了。这6天内，每天调试15个小时以上，中间有十几次以为成功了但是实际上失败了，最终在15号修改了载荷坐标变换程序后，实现了耦合。
其实，这中间还有一个小插曲，因为学术进展不顺利，偷偷玩王者荣耀麻痹自己，3月9号被老婆发现了，在老婆把房子拆了之前赶紧认错收心开始做学术，所以才有了F2A。在此代表F2A的用户感谢我美丽睿智又温柔贤惠的妻子——Jeyla女士。
P.S.：上一次因为偷偷玩游戏被她发现，为了认错逼着自己修改FAST代码，那也是第一次修改FAST代码。后面在FAST基础上开发了SAF，才有了我的博士论文。至此，Jeyla正式锁定“我的一生挚爱”这个任务艰巨的荣誉称号。
六、写在最后
从3月至6月，在F2A强力的尾迹干涉下，写了1篇会议和4篇期刊论文（2xRENEW，ECM和OE），目前都发表了。ECM是我目前发表的影响因子最高的期刊，研究了漂浮式风-流集成能源系统的耦合效应，是我看到的唯一考虑风力机气动弹性与平台和流能装置水动力全耦合的数值研究。不过学术高光已过，从8月份开始就陷入了找工作的烦恼旋涡中，至今还在漂泊。。。
基于OpenFAST漂浮式风力机耦合仿真，实现漂浮式风力机叶片/塔架/系泊等结构疲劳计算及分析
本教程支持试看、回看和同步资料下载，提供答疑专栏90天服务和VIP群订阅用户交流，共享学习资料。另外平台开具电子发票。并且后续根据订阅用户问题不定期加餐直播。
直播《漂浮式风力机和风-流集成能源装备仿真技术》，以下是课程安排，敬请期待！
宁波大学副教授杨阳：漂浮式风力机和风-流集成能源装备仿真技术（8月26日）
最后，想对跟我曾经一样如今正在学术中挣扎的研究生萌新们说一句：没有困难是时间克服不了的，只有投入精力才能解放压力。未来，让我们一起加油吧！
（完）
