柔性风力发电机叶片中流固耦合的计算方法

《Journal of Fluids and Structures》：Computational methods for fluid-structure interaction in flexible wind turbine blades

　　

摘要：为提供准确且计算高效的流固耦合（FSI）研究方法，本研究开发了一种基于弱耦合策略的CFD/CSD耦合模型，并通过全自主开发的求解器实现。在CFD部分，采用基于重叠网格方法的双网格/双求解器策略，以发挥不同数值技术的各自优势并规避其个体局限性。在近体区域，采用非结构混合网格有限体积方法求解RANS方程，以便以相对较低的网格代价准确模拟叶面上的边界层流动。与此同时，在远体区域，采用自适应笛卡尔网格上的五阶WENO格式，以提高尾流和梢涡的分辨率；笛卡尔网格的规则性使得高阶格式能够高效实施，而无需生成全域高阶弯曲网格。在CSD部分，基于几何精确梁理论开发了柔性叶片的多维非线性结构动力学模型。此外，应用有限旋转法来处理CFD网格变形。通过将高阶格式局限于远体区域并仅在需要处采用自适应网格加密，所提出的方法在提高尾流分辨率的同时避免了全局高阶方法带来的过度计算开销。该模型具有精确的气动计算、改善的尾流分辨率以及对几何非线性结构特征建模等优势。

引言：随着风能理论成熟和制造技术进步，风电机组正在显著大型化。截至2023年底，首台16兆瓦风电机组的叶片长度已达123米，其转子扫掠面积相当于约七个足球场。叶片长度的增加使其柔性增强，导致更为显著的流固耦合效应。对大型风电机组叶片流固耦合现象认知不足已成为其设计中的关键瓶颈。因此，开发准确的流固耦合计算方法对于理解这些现象至关重要。

流固耦合计算方法可分为整体式和分区式。整体式耦合方法将流体和结构视为集成系统并同步求解。虽然该方法具有方程形式简便、同步求解和计算精度高优点，但其实施面临重大挑战。流体和结构物理性质的显著差异导致求解复杂度呈指数增长。此外，由于流体和结构不分别求解，该方法在代码开发和模块化方面存在困难，缺乏灵活性。相比之下，分区式耦合方法分别求解流体和结构域，更适合代码开发以及集成现有的计算流体力学（CFD）和计算结构动力学（CSD）求解器（Farhat等，2006年）。分区法通常依次求解流场和结构响应，并在两个域之间交替迭代。该策略已成为主流的FSI求解方法，包含三个关键组成部分：CFD计算、CSD计算以及二者之间的耦合交互。

（1）CFD计算方法

最普遍的方法采用非结构网格有限体积法（FVM）求解雷诺平均Navier-Stokes（RANS）方程。该方法因其算法稳定性和鲁棒性而受到青睐，同时还具备低网格代价、可处理贴合体面网格以及准确表达复杂几何形状等优势。

RANS方程的求解通常采用基于压力的方法或基于密度的方法。在风电机组应用中，流动状态为不可压，通常优先选择基于压力的方法。在EllipSys3D和W2A2K3D等专业软件中，SIMPLE和PISO等算法被广泛使用。然而，为了开发能够模拟可压和不可压流动的通用求解器，采用基于密度的方法结合低速预处理技术是一种有效的替代方案。SU2、OVERFLOW和NSU3D等软件包采用这些技术来模拟低速流动。研究人员已广泛使用这些工具开展研究。例如，Zahle等（Zahle等，2009年）使用EllipSys3D研究了机舱塔架与转子之间的耦合效应；Kirby等（Kirby等，2019年）使用W2A2K3D研究了风电场尺度效应；Yirtici等（Yirtici等，2019年）使用SU2分析了结冰风电机组叶片的气动特性；Ravishankara等（Ravishankara等，2020年）使用SU2研究了叶片前缘粗糙度对风电机组性能的影响。此外，Duque（Duque等，2003年）和Potsdam等（Potsdam和Mavriplis，2009年）分别使用OVERFLOW-D2和NSU3D计算了Phase VI风电机组的气动特性，为后续研究提供了有价值的CFD数据。

非结构网格有限体积法仅能达到二阶空间精度，这对于以涡流为主导流动的保真度模拟来说是不够的。对于此类复杂流动，FVM可以引入大涡模拟（LES）或直接数值模拟（DNS）等方法，而非仅依赖传统湍流模型。然而，LES和DNS因计算效率低而限制了其工程应用的实际性（Slotnick等，2014年）。或者采用高阶空间离散格式具有优势，因为它们能减少数值耗散并提高流场分辨率。不过，提高非结构网格FVM的精度会扩大计算模板，从而导致计算需求显著增加（Zhang等，2013年）。相比之下，笛卡尔网格的规则性使得高阶格式能够更高效地运行，并且便于随流场演化进行网格自适应。然而，粗笛卡尔网格可能难以准确捕捉固体表面的几何特征。此外，实施高阶格式有效地需要在固体表面生成高阶弯曲网格（Whitehouse等，2022年）。总体而言，非结构网格FVM和笛卡尔网格各有优势和劣势，其选择取决于具体的模拟需求和可用的计算资源。

为了发挥各自优势并克服各自的局限性，双网格/双求解器方法为平衡计算精度与效率提供了一种有前景的策略。该方法在近体区域和远体区域分别采用不同的网格和求解技术。例如，Lakshminarayan等（Lakshminarayan等，2018年）研究了机身附近粘性网格层数对直升机旋翼流场模拟的影响，他们在近体区域使用非结构网格FVM，在远体区域使用自适应笛卡尔网格有限差分法。Kamkar等（Kamkar等，2011年）探索了涡流识别方法，在近体区域使用低阶格式，在远体区域使用高阶格式。Sitaraman等（Sitaraman等，2011年）使用双网格/双求解器框架研究了直升机的流固耦合问题。然而，双网格/双求解器方法仍处于起步阶段，在风电机组气动特性和流场计算研究中的应用相对有限。

（2）CSD计算方法

作为风电机组的关键组成部分，叶片随着更大规模风电机组的出现变得越来越柔性，需要精确的结构动力学模型。建立这些模型的主要方法包括模态分析法（MAM）和有限元法（FEM）。MAM计算效率高，但通常忽略叶片变形中的几何非线性，且专注于低阶模态。相比之下，FEM计算精度更高，其建模方法可分为基于梁的模型和三维有限元模型。

风电机组叶片以其细长结构为特征，轴向长度显著大于其他维度，可视为典型的梁结构。在几何非线性梁模型中，绝对节点坐标法（ANCF）（Shabana和Yakoub，2001年；Yakoub和Shabana，2001年）和几何精确梁理论（GEBT）（Simo，1985年；Simo和Vu-Quoc，1986年）是突出的代表。ANCF引入广义坐标，导致节点自由度数量增加——两节点空间梁单元包含24个节点自由度——这影响了计算效率。相比之下，GEBT在计算效率和精度方面均具有优势。它对有限旋转进行参数化处理，以建立大旋转下的应变-位移关系，有效地描述梁截面的旋转。然而，GEBT对插值和积分方法提出了严格要求。与基于梁的模型相比，三维有限元方法采用壳单元或实体单元对风电机组叶片进行详细模拟，提供了更精确的计算结果。然而，更高的精度也意味着更高的计算成本。

由于GEBT具有高精度以及相比三维有限元法显著增强的效率，它得到了广泛应用（Zhang和Huang，2011年）。Sabale（Sabale和Gopal，2019年）研究了风电机组叶片的动态响应；Tang等（Tang等，2021年）探讨了浮动式海上风电机组的水-气-弹-控耦合问题。Li等（Li等，2020年）建立结构模型研究了非稳态来流条件下（包括剪切、塔影和偏航）的流固耦合效应。GEBT也是工程应用中的基础模型。Wang等（Wang等，2015年）基于GEBT开发的BeamDyn代码已被集成到OpenFAST等风电机组模拟软件中（Jonkman和Kelsey，2021年）。SHARPy（Carre等，2019年）是一款面向高柔性飞行机翼流固耦合问题的开源软件，同样采用了GEBT。GEBT根据独立变量和依赖变量的选择存在多种形式，Hodges等（Hodges，2006年）提出的混合形式几何精确梁是风电机组中非线性计算的主要方法。

（3）分区式流固耦合方法

分区式耦合方法可进一步分为强耦合和弱耦合方法。强耦合方法在每个时间步内执行多次迭代计算，从而实现更高的时间精度并减少误差累积。然而，该方法显著增加了计算成本。相比之下，弱耦合方法每个时间步仅需一次迭代，计算成本和时间精度较低，但仍能满足计算要求（Ahn和Kallinderis，2006年）。在CFD的双时间推进方案中，强耦合方法需要在每个虚拟时间步进行耦合，而弱耦合方法只需在每个物理时间步进行耦合。为了提高弱耦合方法的时间精度，研究人员在双时间推进框架下开发了改进版本（Farhat等，2003年；Kim等，2005年）。这些修改后的方法在每个物理时间步进行两次耦合，代表了强耦合与弱耦合之间的中间形式。该改进方案已证明能达到二阶时间精度。此外，它消除了在每个伪时间迭代中进行数据交换的需要（这是强耦合方法的要求），从而在不牺牲精度的情况下显著提高计算效率。图1展示了弱耦合和强耦合策略的比较流程图。

在CFD/CSD耦合模拟中，数据插值方法是关键的计算组成部分。主要方法包括表面跟踪法（Guruswamy，2002年）和表面拟合法（Smith等，2000年）。表面跟踪法通过将流体网格节点正交投影到结构网格单元上来操作，并通过等参变换或面积坐标公式建立插值关系，包括主要用于三角形单元的恒体积变换法（Cebral，1996年）、快速映射算法（Goura等，2001年）和主要用于四边形单元的反向等参映射技术（Murti和Valliappan，1986年）。相比之下，表面拟合法通过对结构网格点进行插值来重建结构表面，以确定相应的流体网格位置，如在薄板样条插值方案中所实现的。

（4）研究空白与目标

当前风电机组流固耦合模拟工具大致可分为三类：低保真工程软件、开源多保真度耦合框架和自研高保真度分区式流固耦合方法。DNV Bladed、DTU HAWC2和NREL OpenFAST是典型的低保真工程模拟工具。在这些工具中，流体域依赖叶素动量（BEM）理论、涡格法或动量线模型进行气动计算，而结构域采用基于Timoshenko梁或几何精确梁理论（GEBT）的多体动力学模型，流体和结构域之间的数据交换通过弱交错耦合机制实现。此类工具表现出卓越的计算效率，被广泛用于整个风电机组的快速载荷评估；然而，其主要缺点在于无法解析叶片边界层内的细流细节和尾流场。

为了弥补低保真工具气动精度的不足，世界各地的研究人员开发了分区耦合模拟系统，如EllipSys3D-HAWC2和OF2（OpenFAST + OpenFOAM）。

在 EllipSys3D-HAWC2 中，有限体积不可压 Navier-Stokes 求解器 EllipSys3D（由 DTU 开发）以分区耦合方式与气动弹性求解器 HAWC2 进行耦合。HAWC2 输出的节点位移驱动 CFD 网格的变形，随后 EllipSys3D 计算更新后的流场。另一方面，OF2 专为浮式海上风机（FOWTs）的高保真气动-水动-伺服-弹性耦合仿真而设计。它采用 OpenFAST 求解转子气动载荷、结构柔性响应和控制系统，同时使用 OpenFOAM 模拟水动力和平台运动，从而有效处理风暴、偏航漂移和系泊线断裂等大位移、强非线性场景。此外，开源多物理场耦合库（如 preCICE）可以灵活地将 OpenFOAM 流体求解器与 CalculiX 结构求解器进行耦合；其官方适配器支持多种数据映射方案和耦合策略，已被多位研究人员应用于模拟风机转子的 FSI 问题。

对现有仿真工具的综合比较揭示了不同方法之间清晰的分层结构：低保真工具在计算效率方面具有主要优势，非常适合工程级的快速载荷计算；高保真 CFD-CSD 耦合框架旨在高保真地重现详细的流场和结构变形响应。然而，传统的高保真耦合方案通常存在固有的局限性，难以同时实现叶面附近的边界层分辨率、远场尾流分辨率以及大叶片变形条件下的计算稳定性。

本研究的目的是开发一种平衡计算效率与精度的 FSI 模型，强调精确的气动计算、改进的尾流分辨率以及叶片结构的几何非线性特征。对于 CFD 部分，采用基于重叠网格方法的双网格/双求解器方案。这一策略背后的理由是充分利用不同类型网格和数值格式的互补优势，同时减轻各自的缺点。

具体而言，对于近体区域，采用非结构混合网格 FVM 求解 RANS 方程。这种方法适合以适度的计算代价捕捉复杂的叶片几何形状和边界层流动，因为非结构网格可以高效地适应曲面而不需要过多的网格单元数。对于远场区域，在自适应笛卡尔网格上采用高阶数值方法（五阶 WENO）以提高远场尾流和涡结构的分辨率。笛卡尔网格的规则性允许高效实现低数值耗散的高阶格式，而自适应网格加密（基于 Q 判据）将计算分辨率仅集中在感兴趣的区域——如叶尖涡和尾流剪切层——而非对整个域进行均匀加密。与全局高阶方法或均匀细网格 LES/DNS 方法相比，这种局部化的方法显著减少了总网格单元数和计算成本。

在 CSD 部分，基于 GEBT 为柔性叶片开发了非线性结构动力学模型。此外，还采用了有限旋转方法来处理 CFD 网格变形。重叠网格框架进一步促进了网格变形，允许近体网格在背景笛卡尔网格内独立移动，消除了通常单片网格系统所需的全局网格重新生成或复杂变形算法的需要。这一特性不仅提高了大变形条件下的计算稳定性，还减少了 FSI 迭代中与重复网格更新相关的开销。

章节片段
CFD 方法
(1) 近体：非结构混合网格 FVM
非结构混合网格 FVM 基于 RANS 方程的积分形式（Blazek, 2001）开发：

?/?t ∫Ω W dΩ + ∮?Ω (Fc ? Fv) dS = ∫Ω Q dΩ

在该方程中，W 表示守恒变量，t 表示积分时间，Fc 表示对流通量，Fv 表示粘性通量，Q 表示源项。在非结构混合网格 FVM 中，采用 Harten-Lax-van Leer-Contact（HLLC）格式（Toro 等，1994）计算对流

NREL Phase VI 风机的空气动力学特性
由美国国家可再生能源实验室（NREL）开发的 NREL Phase VI 风机（Hand 等，2001）提供了全面可靠的实验测量数据，使其成为验证风机气动方法和仿真工具的国际基准。表 1 总结了转子的基本几何参数，更多规格请参阅参考文献（Hand 等，2001）。来流速度为

计算模型与网格
(1) 计算模型
本研究关注 NREL 5 MW 风机（Jonkman 等，2009），这是一种概念性的三叶片模型，包含全面的气动外形和结构特性参数。该风机模型已被广泛应用于风机应用中的方法验证和工程研究。风机每片叶片长 63 米，最大叶尖偏转可达叶片长度的 10%。基本形状参数详见表 2。

(2) 结论
本研究开发了一种平衡计算精度与效率的柔性风机叶片流固耦合数值模型。该模型采用松散耦合策略，将双网格/双求解器 CFD 模块与基于 GEBT 的 CSD 模块集成，以实现流体与结构动力学域的协同仿真。在 CFD 部分，建立了近体/远体分区框架。近体区域采用非结构混合网格有限

作者贡献声明
Pan He：初稿撰写、可视化、软件、方法论、调查研究。Jian Xia：概念构思、方法论、软件、资源、项目管理、指导、审阅与编辑。

Pan He | Jian Xia