风力机叶片流固耦合计算的降阶模型研究

孙芳锦，祝东涵，张大明

(1.广西岩土力学与工程重点实验室，广西 桂林 541004；2.桂林理工大学 土木与建筑工程学院，广西 桂林 541004；3.辽宁工程技术大学 土木工程学院，辽宁 阜新 123000；4.桂林理工大学 信息科学与工程学院，广西 桂林 541006)

风能具有蕴藏量丰富、可再生、无污染等特性，成为新能源重要的发展方向。风力发电作为风能利用的最主要形式，近年来得到快速发展。风力机叶片是风力机的发电来源，其性能直接影响发电机的效率。目前风力机叶片柔性大，有较强的几何非线性，在运行工况下，如阵风作用等都会引起风轮叶片发生变形和振动，而振动产生的变形又会改变空气流场分布，产生了风荷载与风力机叶片的流固耦合作用。流固耦合作用会引起叶片的气弹变形，可能导致叶片振动过大，或振动失稳，甚至导致疲劳破坏[1]。因此建立准确方法研究风荷载与风机叶片的流固耦合作用，并分析叶片气弹振动的变化规律，对于为风力机叶片的合理设计、保证风力机的安全高效运行具有重要意义。

随着计算机软硬件技术的迅速提高，数值模拟方法已发展成为研究风荷载与风力机叶片流固耦合作用的重要工具。由于风力机叶片外形复杂、非线性强，其流固耦合计算考虑的自由度多，通常会耗费大量机时。降阶模型是研究流固耦合问题的重要手段之一[2-3]，降阶模型属于经验谱方法，假定流场包含一系列特征值和特征向量，其最重要的特性可以通过降阶模型获得。建立降阶模型可以有效减小模型尺寸，节省计算机时和计算资源，因此采用降阶模型进行风力机叶片的流固耦合计算，对于准确进行风力机叶片设计有着重要的意义。流固耦合计算方法可以分为强耦合法和弱耦合法，它们指的是空气流体控制方程和结构控制方程间的耦合关系，弱耦合法是在每一时间步内先对流体控制方程和结构控制方程分别独立求解，然后将作用在流体域模型上的气动力荷载通过流体和结构的交界面传递给结构域模型，用来预测结构的位移，然后再将结构的位移作为新的荷载传递给流体域，此过程如此反复，直到结果收敛于指定值。强耦合方法是指在每一时间步内对流体控制方程和结构控制方程同时联立求解，即在每一时间步内流体域和结构域互相传递气动力和结构位移的过程是同时进行的。对于风力机这种大型柔性结构的流固耦合计算，弱耦合法更加实用，本文即针对弱耦合法，建立风力机叶片流固耦合计算的降阶模型。

建立降阶模型的第一步是进行流场分解，分解方法包括正交分解法(POD)，Vorterra级数方法和谐波平衡法等。其中正交分解法(POD)是比较常用的方法，它是从时域或频域计算中提取结构的动力信息[4-6]，为了从试验或数值模拟中提取流场特征向量，Sirovich[7]提出了“快照”的概念构建降阶模型。目前国内外对降阶模型在风力机叶片流固耦合计算中的研究还非常有限，降阶代理模型曾被用于预测失速条件下风力机的非稳态气动力荷载[8]，本征正交分解的降阶模型曾用于复杂流体动力系统的低维建模[9]，文献[10]提出了风力机叶片非定常气动力的降阶模型方法，建立了叶段振动状态下的非定常气动力模型,用以模拟叶片变形与气流耦合作用下的附加非定常气动力,给出多叶段模型的气动弹性响应历程，文献[11]基于POD方法，结合边界有限元法和降阶模型对水平轴风力机的气动特性进行了研究。

针对风力机叶片流固耦合计算复杂，耗时量大的局限性，本文采用降阶模型对风力机叶片进行流固耦合计算。针对风力机叶片设计参数多，参数空间维数大的特点，采用本征正交分解(POD)方法建立流固耦合计算中空气流体方程的降阶模型，利用POD方法将流体控制方程的速度和压力等未知量进行基模态分解并用一小部分空间模态表示，对流体方程进行时间离散，再在上述缩减空间上采用Galerkin投影方法获得流体控制方程的最小残差投影，构造流体方程的最小残差降阶模型。再将流体降阶模型与风力机叶片进行流固耦合计算，达到减少计算时间，提高计算效率的目的。

1 基于本征正交分解的最小残差降阶模型

本文研究的是弱耦合问题，这里将采用本征正交分解(POD)建立空气流体方程的降阶模型。其主要思路是利用POD方法将流体控制方程的速度和压力等未知量进行基模态分解，用一小部分空间模态表示，再对各基向量进行投影重构，构造最小残差的降阶模型，达到减少计算时间，提高计算效率的目的。

流体控制方程由为不可压黏性Navie-Stokes方程(N-S方程)，即连续方程和动量方程，

(1)

(2)

流体的边界条件为

(3)

(4)

(5)

(6)

则对于t∈[0,T]，未知量对应的基向量空间可以写作，

(7)

在POD的基向量空间上将流体方程的未知量可以写作，

(8)

为推导流体降阶模型，需要对流体方程进行时间离散，本文采用一阶Euler时间离散[13]，则有

在上述缩减空间上采用Galerkin投影方法获得流体控制方程的最小残差投影，进而获得降阶模型，那么对于∀t∈[0,T],n=1,2,…,N，在时间步n时，流体方程的降阶模型可以写作，

其中ζn+1={an+1,bn+1}T

那么上述函数的变分形式可以写作，

且对于时间系数，该函数应满足，

则最小残差应满足如下关系，

将上述降阶模型与风力机叶片有限元模型进行耦合计算，便可以计算出风力机叶片周围风场分布、风力机叶片的气动性能等参数。

图1给出了本文降阶模型的流程和主要步骤。

图1 降阶模型计算流程图

2 算例分析

本文采用本文降阶模型对NREL V大型风力机叶片[14]进行了流固耦合计算，NREL V风力机是一个数值风力机模型，参考了多种商用风力机原型机，美国国家可再 生能源实验室曾对其流固耦合作用进行过试验研究。风力机为上风向3 叶片风力机，采用变速变桨运行控制方式。风力机叶片的参数如表1所示。

表1 风力机叶片参数[14]

首先采用POD快照方法，组成流场状态矩阵，利用奇异值分解,可得未知量对应的基向量空间。这里基于FLUENT平台，按照上述NREL大型风力机叶片的初始条件，计算流体速度场和压力场的POD基，得到t=0 s到20 s的速度场和相应压力场，共200组即为快照。再应用本征正交分解(POD)，得到基向量及对应的特征值，得到的每一个特征值可以用来刻画所对应的模态所捕捉到的能量，也就是说特征值的大小可以反应该基向量在整个系统中所占的权重。表2分别给出了速度场和压力场的模态特征值和能量累计比例。

从表2数据可以看出，速度场的前15阶和压力场的前14阶特征向量能够捕获系统约99%的能量，因此选择速度场的前15阶和压力场的前14阶特征向量作为降阶模型的POD正交基，作为原模型的替代。

表2 速度场和压力场的特征值和能量比

为证明本文提出的最小残差降阶模型的正确性，这里分别采用本文方法和经典Galerkin法计算了时间离散系数，这里分别给出了速度场和压力场前4阶POD模态的时间系数对比，如图2和图3所示。

从图2和图3可以看出，采用本文降阶模型和经典Galerkin方法计算得到的速度场和压力场POD模态时间系数变化趋势一致，且比较接近，说明本文提出的最小残差降阶模型在计算中是正确的。

(a)一阶模态

(a)一阶模态

采用本文提出的降阶模型计算了风力机叶片流固耦合计算中的受力和变形(风速为11.4 m/s)，分别如图4和图5所示，并将计算结果与已有文献[15]进行了对比，以验证本文降阶模型在风力机叶片流固耦合计算中的准确性。

从图4可以看出，采用本文降阶模型得到的叶片轴向力和切向力与已有文献[11]变化趋势一致，且结果比较接近，可以看出，轴向力在叶片的中部靠后的位置达到最大值，切向力在叶片的中后部达到最大值。观察图5，可以看出，降阶模型计算得到的叶片在挥舞方向和摆振方向的变形与已有文献结果变化趋势一致，看出叶片的主要变形区域发生在叶片的中上部分，且越接近中上部，变形越大。因此，由流固耦合作用引起的风力机叶片大变形必须在大型风力机叶片设计中引起足够重视，并在必要时限制相应的风力机叶片变形。

(a)轴向力

(a)挥舞方向变形

同时本文还采用提出的降阶模型计算了叶尖的变形情况(风速为11.4 m/s)，迎角为30°，图6给出了叶尖随时间变化的变形情况，图7给出了叶尖随方位角变化的变形情况(0°方位角是指叶片垂直向上)，并与CFD计算结果进行了对比。

观察图6和图7可以看出，本文降阶模型计算的叶尖随时间和方位角变化的变形与CFD计算结果变化趋势一致，说明了本文降阶模型的正确性，以及在风力机叶片流固耦合计算中的有效性。从图6看出，随着时间推进，风力机叶尖挥舞运动方向的变形要大于摆振运动方向的变形，叶尖挥舞运动方向的变形随时间增加而增大，在约30 s达到最大时，变化较小，变形基本保持不变。摆振运动方向的变形随时间增加而减小，在约15 s达到最小值后变形基本保持不变，因此在设计时应对叶尖挥舞运动方向的变形加以重视。从图7可以看出，随方位角变化，风力机叶尖挥舞运动方向的变形和摆振运动方向的变形变化趋势比较接近，基本都是在90°和270°方位角时达到最大值。

图6 风力机叶尖变形随时间变化对比

图7 风力机叶尖变形随方位角变化对比

为分析本文方法的计算效率和计算收敛情况，表3对比了不同时间步采用不同模态数(前10个模态)时计算容差收敛的情况。

从表3中可以看出，本文的降阶模型方法收敛速度较快，使用较少的时间步就能达到较小容差标准；且采用较少模态就可以达到较高计算精度，可以看到，当采用6个模态就可将计算容差平均减少到10-3以下，说明本文降阶模型的计算效率高，且计算精度随时间步的增加而增加，反映了本文的降阶模型具有较高精度。

表3 不同时间步时采用不同模态数时计算容差收敛情况

为了进一步说明本文降阶模型的计算效率，本文对比了采用降阶模型和全阶模型时的相关计算参数对比，如表4所示，这里主要对比了在不同计算收敛容差时，计算叶片和叶尖变形时的平均耦合迭代次数和平均耗费计时。

表4 采用不同方法计算叶片变形的平均耦合迭代次数和计算计时

分析表4中的数据可以看出，达到同样的计算容差时，采用本文的降阶模型的平均迭代次数和平均耗时都要比全阶模型少，其中平均耦合迭代次数比全阶模型平均减少约42%，平均耗费机时比全阶模型平均减少了约40%，主要原因在于构造最小残差的降阶模型只是用一小部分空间模态表示各参数，在缩减空间上获得流体控制方程的最小残差投影，因此达到了降低计算成本的目的，说明了本文降阶模型在达到计算准确性的同时还提高了计算效率，是正确可靠的。

3 结 论

本文提出了基于POD方法，建立了风力机叶片的最小残差投影降阶模型，将该降阶模型应用于风力机叶片的流固耦合分析中，取得了较好效果，得到的主要结论有：

(1)采用本文降阶模型和经典Galerkin方法计算得到的速度场和压力场POD模态时间系数变化趋势一致，说明本文提出的最小残差降阶模型在计算中是正确的。

(2)本文降阶模型对风力机叶片的流固耦合计算的受力和变形与已有文献符合良好，证明了本文降阶模型的准确性。同时结果表明考虑流固耦合作用时，风力机叶片的变形呈非线性分布，流固耦合作用对叶片变形有重大影响。

(3)随着时间推进，风力机叶尖挥舞运动方向的变形要大于摆振运动方向的变形，叶尖挥舞运动方向的变形随时间增加而增大，在约30 s达到最大。摆振运动方向的变形随时间增加而减小。在设计时应对叶尖挥舞运动方向的变形加以重视。随方位角变化，风力机叶尖挥舞运动方向的变形和摆振运动方向的变形基本都在90°和270°方位角时达到最大值。

(4)本文的降阶模型方法收敛速度较快，使用较少的时间步就能达到较小容差，且采用较少模态就可以达到较高计算精度。且与全阶模型相比，采用本文的降阶模型进行流固耦合计算时，平均耦合迭代次数比全阶模型平均减少约42%，平均耗费机时比全阶模型平均减少了约40%，大大提高了计算效率。