李燕玲 胡进 周锟

(武汉科技大学省部共建耐火材料与冶金国家重点实验室，武汉430081)

据统计[1]，全球可再生能源(风能，太阳能，以及生物质能等)于2019年占据全部一次能源的5%，并且实现了12.1%年度增长。其中，风能占全部可再生能源的51%，并且以30%的年度速度增长。而中国是可再生能源增长速度最快的国家。风力机是将自然界的风能转换成机械能或电能的装置。其重要工作部件是叶片，决定了风能转换的效率。因此，提高叶片的气动性能是研发风力机的基础与关键，也是难点之一。

风力机也隶属于叶轮机械。叶轮机叶片气动性能研究主要包括理论研究、数值模拟和风洞实验。理论研究发展较早[2-3]，如常用的Betz理论、叶素理论、动量理论以及叶素−动量理论等发展相对成熟。然而理论研究的计算模型、求解条件都比较复杂，并且实际工作环境中存在许多不可预测因素，通常需要对模型简化以得到相关解[4-5]。风洞试验虽然能够准确地控制实验条件，受气候条件影响小，但需要花费大量的时间和资金。随着研究的深入和计算机技术的发展，越来越多的研究人员倾向用数值模拟方法对旋转叶片流场特性进行研究分析。而当前的计算流体力学(computational fluid dynamics,CFD)方法，可分为以下几类。(1)多重旋转坐标系方法[6-7]：在叶片周围使用旋转坐标系，设置旋转域和叶片转速相同，模拟过程中旋转域和叶片之间不发生网格的相对运动，将非定常问题转化为定常问题，可大幅节省模拟计算的时间，但计算精度不是非常高。(2)滑移网格法[8]：需要划分包含风机叶轮的旋转域和剩余的固定域，在旋转的过程中两域在交界面处相对滑动，但此方法需要对交界面进行特殊设置。(3)重叠网格法[9-10]：将复杂的流动区域分成几何边界比较简单的子区域，各子区域中的计算网格独立生成，彼此存在着重叠关系，但在计算域运动的过程中，需要不断检测网格的重合区域，流场信息通过插值在重叠区边界进行匹配和耦合。(4)动网格法[11]：动网格法被运用在瞬态模拟中，风机叶片结构在旋转过程中不再是一个刚体，而会由于受力导致变形，因此每一时间步中均需对计算域中的网格进行变形与重构，确保网格质量满足模拟计算要求。

基于浸入边界法[12]，本文发展了一种简单的建模，网格离散，与数值模拟统一方法框架，来模拟风力机的气动问题。在此框架下，整个风力机(包括叶片，轮毂，机舱，以及塔架)与周围的空气流场所占据的空间都用统一的均匀网格来离散划分。而风力机与空气的界面，采用散点来离散。在这些离散点上，通过对空气流场施加合适的作用力来模拟流固耦合作用。对于叶片构型，本文提出一种简单的同伦变形，来光滑融合不同的二维翼型从而生成更加复杂多样的三维叶片。同时，采用仿射变换，来处理叶片的扭转与渐缩问题。对于所有的风力机部件，本文都给出很简单的离散点生成方法。相较于基于贴体网格的各类模拟方法而言，本文采用的方法在模型离散时更加强健，可以非常容易地处理模型各部件之间的狭缝，区域重叠等常见问题，且数值稳定性强。此方法框架的另一个显著特点，是将建模，网格离散，数值模拟融为一体，这使得有可能在此框架下发展设计选型(关于翼型，尺寸，扭转角等参数)自动优化算法--而这是本课题组的长期研究方向之一。

为了检验方法与程序的正确性，本文模拟了二维典型翼型在不同攻角下的升阻力问题，并且与高精度的谱元方法模拟结果进行了对比。针对浸入边界法中阻力过估明显的问题，提出了一种简单有效的理查森线性外推误差修正方法。同时，模拟研究了拱曲度以及厚度对二维翼型升阻力的影响。随后，模拟研究了单风力机(包含塔架)在不同尖速比下的功率系数，并对塔架与叶片间的相互气动作用进行了初步分析。最后，模拟研究了双风力机在风场中不同前后间隔距离下的气动干涉问题。

1 研究方法

1.1 流场模拟方法——浸入边界法

浸入边界法是一种很常见的处理流固耦合问题的方法。在此方法中，包含固体与流体的全部空间，都采用欧拉坐标下的纳维−斯托克斯方程来描述

式中，u，p，ρf和v分别是流体的速度、压力、密度和运动黏度，∇为哈密顿算子，fIBM表示固体对流体的作用力。

另一方面，对于流体中固体的运动，则是采用拉格朗日坐标系下的牛顿运动方程来描述。比如固体的平动，其方程为

式中，ρp与up分别表示固体的密度与平动速度，而FIBM则表示作用于固体上的作用力。当固体进行旋转时，还要建构类似的转动方程。

理想状态下，固体对流体的作用力−FIBM与流体对固体的作用力FIBM应满足牛顿第三定律，即大小相等，方向相反。因为固体对流体的作用力−FIBM与流体对固体的作用力FIBM是分别定义在两种不同坐标下(即欧拉坐标与拉格朗日坐标)，一般情况下，两者空间位置并不重合，因而需要采用不同坐标间的相互插值来进行交互。最常用的则是采用离散δ函数。

浸入边界法中最后一步，是构建流固间作用力模型。对于刚体，最流行的方式是直接作用力方法[12]

式中，˜U是在固体边界位置处流体的插值速度，∆t是时间步长，即流固间的相互作用力等于流体与固体在界面处速度的差异除以时间步长。换言之，直接作用力是通过施加一个正比于流固速度差的回复力，而使得流固间的不可滑移边界条件得到满足。

本研究中采用的程序，是基于经典谱方法槽道流动求解器，并在其中引入额外的作用力项来模拟流固耦合作用。关于浸入边界法的技术细节以及程序实现，可参照文献[12-14]。

1.2 风力机叶片建模及网格离散

风力机叶片是收集风能的核心部件，其作用原理与机翼非常类似。空气流过叶片时，在叶片上同时产生升力与阻力。风力机叶片多采用高升阻比翼型，亦即更大比例地利用升力作用来驱动叶片转动，从而获得更高的风能捕获效率。出于对效率和结构强度等多方面的考虑，叶片形状非常多样。一般来说，为了提高结构强度，在叶片根部采用更厚重的翼型；而在尖部采用更轻薄的翼型。一方面，在改变翼型的相对厚度的同时，朝着尖端方向叶片的弦长也不断减小，从而使得整个叶片呈现尖锥状。另一方面，沿着翼展方向叶片也呈现不同角度的扭转。这是因为叶片的旋转线速度沿翼展方向不断增大，从而使得空气来流与叶片的相对速度沿翼展方向也不断改变。为了保证相对来流的攻角在整个翼展方向都处于理想范围，需要沿翼展方向适当扭转叶片。

此处，介绍一种简单地生成光滑渐变叶片的方法。假设在叶片最大弦长处(靠近根部)采用翼型B，而在最小弦长处(叶片尖端)采用翼型C。翼型B和C可以从通用翼型库中选择(国家可再生能源实验室官方网站上提供了几十种二维翼型供选择)。每种二维翼型，都可以从其前缘到后缘沿着弧长进行参数化，即以前缘为零点到前缘的表面弧长为坐标。这些坐标，都采用总弧长作为归一化参数，从而使得坐标从0到1之间变化(本文将此坐标记为l)。在翼型B和C之间，采用同伦变形来生成中间的翼型，使得翼型从B光滑变化到C。此同伦变化数学上可表示为

式中，D表示位于B与C之间s处的翼型，而s是沿翼展方向距离B的归一化距离(即s=0表示B的位置，而s=1则是C的位置)。上面的同伦变形，给出了在任一位置s，弧长坐标l处的翼型坐标，从而唯一确定了B与C之间所有的翼型。这些翼型从B沿翼展方向光滑变化到C。这种构建复杂翼型的方法非常简单，且可直接推广到更复杂的情形，即在叶片多处位置指定不同的二维翼型，而在任意两段翼型之间进行同伦变形来生成光滑过渡翼型。

此处介绍的同伦变形翼型构建方法，还可以非常便捷地生成叶片的表面离散网格。作者关于浸入边界法的前期研究表明[14]，当固体的表面离散网格点呈现均匀分布，与整个空间的均匀欧拉网格相仿时，将具有更好的数值精度。所以本文研究模拟的叶片(包括后面讨论的风力机其他组件)，都采用尽可能均匀的表面离散点。生成叶片表面离散点时，首先给定离散点之间的需要间隔dx，此间隔应该与欧拉网格的尺寸相仿。然后沿着翼展方向，以此间隔dx，垂直于翼展方向截取叶片形成一系列的截断面。对于每个截断面，先计算出其总周长，再以总周长除以间隔dx来确定所需的离散点数目，最后沿着周向等间距生成此数目的离散点。

对于真实叶片，还有两个参数需要考虑：锥度与扭转角。这两个参数都可以采用仿射变换来统一处理。仿射是包含平移、旋转、镜像和放缩等一系列几何相似变形的组合。在数学上可以通过将原始几何坐标左乘一个变换矩阵来实现。对于叶片而言，所需要的变换是在沿着翼展方向的每个截断面上，对每个翼型周围的离散点进行适当的放缩和旋转，从而获得所需的锥度和扭转角。这种仿射变换表达式为(不失一般性，假设二维离散点处于x−y平面上)

式中，S为由离散点x，y坐标以及数字1填充而成的3×N矩阵(N是离散点数目)，S′则是对应的变形之后的离散点坐标矩阵。方程右边三个矩阵分别表示离散点需要的平移量xt和yt，旋转角度θ，以及放缩系数W。平移量直接由叶片最终的空间位置确定，旋转角度就是所需的扭转角，而放缩系数则由锥度确定。一般来说，叶片的弦长沿翼展方向大多呈线性变化(固定锥度)，此时放缩系数沿整个翼型方向都是固定值；但旋转角θ沿翼展方向通常是变化的。值得特别说明的是，此处采用了扩展的变换矩阵来统一处理平衡、旋转与放缩问题。在最后结果S′矩阵中的第三行数据是冗余的，只有第一与第二行数据分别包含所需的变换后离散点x与y坐标。

1.3 风力机轮毂，机舱，塔架建模及网格离散

风力机轮毂，机舱以及塔架等的气动力直接作用远小于叶片。在本文中，轮毂、机舱和塔架模型分别采用简单的半椭球体、长方体和细长圆台。对于半椭球及圆台这两种中心对称的几何体，其离散方法与前面的叶片沿翼展方向的离散完全相同，而对于长方体表面则可直接采用正交均匀网格离散。

2 风力机几何模型及离散点示例

本文将风力机置于经典槽道流动中来模拟风力机的气动力问题。槽道横向、流向以及垂直方向分别记为x，y和z。槽道的底部设置为不可滑移的零速度边界条件来模拟地面。顶部采用给定的速度边界来模拟稳定的空气来流。在其他方向都采用周期性边界(因为采用的傅里叶谱方法)。对于叶片，选择NREL S825和S805A两种翼型(分别对应图1中B和C)，按照前文介绍的同伦变形方法来生成(其三维投影见图2)。图1中O对应于30%弦长位置，是翼型扭转中心，同时也是叶片的装配中心。两种翼型的方向对应于各自的扭转角(叶片旋转平面为x−z)。翼型S825具有更大的相对厚度，使其位于叶片近根部具有最大弦长的位置(图2中B位置)。而相对厚度较小的翼型S805A则位于叶片尖端(图2中C位置)。翼型S805A的弦长设置为最大弦长的一半，使得整个叶片形成锥状。沿着翼展方向叶片不断扭转，其扭转角(翼弦与叶片转动平面的夹角)采用文献[15]中给出的非线性分布。在叶片根部扭转角最大为30°；而在尖端扭转角较小为−1.9°。除了B至C段的渐变翼型，叶片模型的圆形的装配底座(图2中A位置)以及从A到B的过渡段的离散点同样采用同伦变形方法来生成。整个叶片，从转轴到翼尖为5(本文所有长度都采用无量纲化数值，以叶片的最大弦长为标准长度)。以旋转轴为零点，装配底座沿翼展方向位于0.5处，最大弦长位置为1.25。图2给出了整个叶片的真实比例结构图和离散点间隔dx=0.1条件下的表面散点图。在此条件下单个叶片大约有860个表面离散点。

图1 叶片最大弦长及最小弦长处的二维翼型

图2 叶片模型尺寸以及表面离散点

整个风机(轮毂，机舱以及塔架)的装配模型以及离散点如图3所示(全部风力机表面离散点约为6100个)。不同于贴体网格类型的CFD方法，浸入边界法对几何模型的表面离散要求非常低，并不要求相临离散单元之间满足任何拓扑相洽关系。换言之，不同部件离散单元之间重合和穿透等，都不会给数值方法本身造成任何困难。对于风力机这样复杂的结构，可以通过本文介绍的方法，在不依赖任何复杂的几何建模和网格生成软件的前提下，进行成功模拟。

图3 风力机整机装配模型以及表面离散点

3 模拟结果

3.1 二维机翼绕流

为了验证数值模拟结果的精度，本文模拟了不同网格精度、不同来流攻角α条件下二维翼型NREL S805A的绕流问题。分析攻角分别为0°，5°，10°，15°以及20°，雷诺数Re(雷诺数定义为来流速度乘以弦长并除以运动黏性)为10和100的升阻力结果。值得指出的是，此处模拟的雷诺数非常低，远低于风力机工程实际应用中的雷诺数。因为本研究所有结果都是采用直接数值模拟，而众所周知，直接数值模拟方法的计算量近似正比于雷诺数的三次方[16]，使得无法真正模拟实际工况。此算例主要是用来便捷地检验算法的网格依赖性与精度。在后文模拟的三维风力机，雷诺数更加接近工程实际，并将雷诺数对结果的影响进行了探讨。

作为对比标准，本文还采用了基于贴体网格的高阶精度谱元方法Nek5000[17]，模拟了相同工况下翼型的升阻力。图4给出了Re=10条件下的阻力系数(Cd=2Fd/(ρfu2∞c)，其中Fd是阻力，u∞是来流的速度，c是弦长)与升力系数(Cl=2Fl/(ρfu2∞c)，其中Fl是升力)随攻角的变化。在每一幅子图中，都包含5组数据，分别对应于4种不同的网格精度与谱元方法的结果。这4种网格精度分别为dx=0.05，0.025，0.012 5，0.006 25。对于升力系数，浸入边界法与谱元方法升力结果吻合很好。这说明浸入边界法对于升力(也就是翼型上下表面的压差)模拟非常准确。但是对于阻力，则与谱元方法结果相差明显，这种差距明显依赖于网格大小。

图4 不同攻角α，不同网格尺寸dx条件下二维翼型S805A阻力与升力与谱元方法Nek5000结果对比

为了进一步研究阻力与网格尺寸之间的关系，图5给出了不同网格大小下浸入边界法阻力与谱元方法结果的相对误差表示谱元方法阻力)。图中还包括Re=10条件下的结果。从图中可以清楚看出，不论是在何种攻角或者雷诺数条件下，相对误差都几乎正比于网格尺寸dx。即浸入边界法对于阻力的模拟，相对网格尺寸具有较严格的一阶精度。这与作者之前的平板黏性阻力理论研究结果[14]完全相符，并进一步表明之前的理论结果对于翼型这种复杂的构型也适用。一方面，因为风力机属于大升阻比类机械，主要是依赖升力驱动，而升力的精确模拟表明浸入边界法可以较精确地模拟风力机的功率。另一方面，可以利用不同网格尺寸下的阻力数值，通过线性外推法来预估理想条件下(即网格尺寸趋近于零)的阻力。具体来说，可以选定两种不同尺寸网格重复模拟同一问题。以图5中Re=10结果为例，可以使用dx=0.2，0.1这两组最粗网格下的模拟结果来构建线性的误差曲线。而此曲线与纵轴的相交点，即为使用IBM方法在极限条件dx=0时的误差。采用此线性外推方法，得到的结果误差不到1%(图6中数据线性外推，与坐标原点非常接近)。对于Re=100的结果，采用dx=0.05，0.025这两组数据，线性外推的结果误差约为3%；如若采用dx=0.025，0.012 5这两组网格数据线性外推，误差则下降到1%以下。总的来说，雷诺数增大时需要更精细的网格来保证精度；而采用线性外推法，可以使用相对稀疏的两种不同尺寸网格，得到精度远高于采用单一精细网格进行模拟所得的结果。因为对于二维问题，网格加倍时计算量上升4倍；而对于三维问题，计算量上升8倍。就图5中的二维问题，如果采用极其细密网格来得到与线性外推法同样精度的结果(误差1%以下)，所需要的计算量则上升2到3个数量级。由此可以看出，这种简单的线性外推法对于黏性阻力的模拟误差纠正，是极其有效的。而这种方法能成功的理论根源，则在作者前期工作中进行了深入阐述[14]。此处的模拟结果，则进一步验证了前期理论的普适正确性。值得指出的是，这里的翼型模拟结果仍然限制于边界层产生分离之前的流动情况。在更大雷诺数下，边界层产生局部分离时，前述的线性外推法是否依然精确高效，将是进一步要研究的工作。

图5 不同网格尺寸dx以及攻角条件下，二维翼型阻力相对误差

图6 不同叶尖速比下力矩系数与功率系数

二维翼型的拱曲度(即拱起来的程度)与厚度(即最大内切圆的直径)是影响其升阻力的两个主要形状参数。此处，模拟了三个不同拱曲度的翼型，即S825，S805A，以及S809(三者拱曲度从大到小)。模拟结果表明，翼型的拱曲度越大，升力系数越大，这与文献[20-21]中的研究结果一致。比较三种翼型的厚度，S809厚度最大，其次是翼型S825，最小的是翼型S805A，且翼型S809的厚度明显大于S825和S805A，S825翼型厚度略大于S805A。翼型S805A属于薄翼型，而S825和S809属于厚翼型。三种不同厚度翼型在不同雷诺数及来流攻角下的模拟结果表明，当雷诺数越小时，厚度较大的翼型升阻比越大。雷诺数变大时，在小攻角下，翼型的厚度对升阻比的影响不敏感，但随着攻角的增大，厚度较大的翼型首先出现失速现象。

3.2 单风力机流场模拟

基于前文介绍的方法与模型，模拟了三叶片风力机在不同叶尖速比λ(翼尖转动线速度与来流风速之比)下的流动。本研究中，基于来流与叶片长度的雷诺数为2.5×104，这比工程应用中的实际雷诺数要低两个数量级左右。这是由于本研究中采用直接数值模拟，无法实现对应于工程应用中的高雷诺数。尽管如此，也有理由相信此处关于风力机力矩与功率的数值模拟结果与实际工况比较吻合。这主要是因为风力机升力为主型机械，而在本文模拟的雷诺数下，叶片绕流的整体结构与工程实际非常类似，从而使得叶片上模拟的压力分布具有较高的可信度(尽管阻力误差较大)，最终保证叶片上的升力模拟所具有的精度。事实上，关于类似翼型的风洞试验结果表明[22]，在湍流绕流条件下，不同攻角下翼型升力系数在整个实验范围内(雷诺数从6.0×104至4.6×105)几乎完全相同，而只有阻力系数才对雷诺数呈现出较强的依赖性。

图6给出了不同叶尖速比下叶片上的力矩系数(Cm=2M/(ρV2SR)，其中M，S，R分别是风力机所受力矩、风轮扫掠面积和风轮半径，ρ和V分别是流体的速度和来流风速)与功率系数(CP=Cmλ，其中Cm和λ分别为力矩系数和叶尖速比)。虽然叶片阻力随着叶尖速比增大近似线性增大(此处没有显示模拟结果)，但对力矩系数以及功率系数的影响则完全不同。在叶尖速比为1∼3时，风力机的功率系数加速增加，叶尖速比为3∼5时，功率系数随叶尖速比的增加变缓且在叶尖速比为5时达到最大值0.43，然后随着叶尖速比的升高，风力机发生失速现象导致功率系数下降，即此条件下风力机的最佳叶尖速比为5，变化趋势与文献[18-19]的结果相同。

为了研究塔架对风力机的影响，本文还模拟了没有塔架的风力机。对比两者结果发现，就本文模型而言，塔架所受阻力约为总风力机的30%，而叶片所受升力为7%。而塔架的作用，是使得叶片的上力矩减小了0.6%。总的来说，因为塔架本身较大的风载荷，必需考虑这部分载荷的强度要求；而塔架对风力机功率的影响则较小，可以忽略不计。

3.3 双风力机流场模拟

在风力场中，多风力机的气动干涉问题，是决定风力机布局的关键因素之一。一方面，希望风力场占地尽可能小，要求风力机布局尽可能紧密；但另一方面，为了保证单风力机的风能效率，又要求风力机之间相互影响尽可能小。本文采用的浸入边界法，可以很便捷地模拟多风力机的气动干涉问题。只需将前文介绍的风力机表面离散点，复制放置于风力场中不同位置，便可模拟多风力机的流动问题。在此过程中，无需对均匀的背景欧拉网格进行任何处理。

图7给出了在叶尘速比为1，雷诺数为2.5×104条件下，两风力机前后间隔为3D时(D为叶片直径)垂直截面瞬时速度云图。从图上可以清楚看出地面，塔架等对空气来流的阻碍作用。为了分析前风力机尾流的影响，图8给出了双风力机在前后间隔分别为3D与5D时，沿叶片旋转轴平均流动速度图。第一台风力机后，轴线处的速度从零(近轮毂处)开始迅速增加，然后趋于稳定；在接近第二个风力机时快速减小到零，而在离开第二个风力机后又快速增加。在风力机间隔为3D条件下，第一个风力机后尾流的最大速度恢复到空气来流的63%。而在间隔为5D的条件下，最大速度恢复到来流的79%。因为风力机的最大理想输出功率正比于空气来流速度的三次方，那么基于轴线上来流风速进行简单估计，在风力机间隔为3D时，第二个风力机功率仅为第一个风力机的25%；在间隔为5D时，第二个风力机功率仅为第一个风力机的一半。所以，为了保证风力场的整体风能效率，要尽可能避免风力机位于其他风力机的近尾流区。

图7 双风力机速度云图(垂直剖面)

图8 前后双风力机沿轴向气流平均速度

4 结论

浸入边界法是一种非常常见的多相颗粒流直接数值模拟方法。本文将浸入边界法应用到风力机的气动力学数值模拟中，建立了从建模，网格离散，到数值模拟分析的统一框架模型。首创性地提出采用同伦变形和仿射变换的方法便捷地处理复杂叶片的建模问题。在模拟二维翼型的升阻力过程中，通过与高精度谱元方法结果对比，发现浸入边界法相对空间离散具有较严格的一阶精度。这一方面数值检验并扩展了作者的前期理论研究成果[14](即浸入边界法在平板边界层模拟中具有严格的一阶精度)，证明在具有复杂压力梯度的边界层中，浸入边界法仍然具有一阶精度。另一方面，基于不同网格下的升阻力模拟结果，提出了一种简单的理查森外推法来修正结果，从而得到远高于一阶精度的模拟结果。

对不同拱曲度以及厚度的二维翼型模拟结果表明：翼型的拱曲度越大，升力系数越大；当雷诺数越小时，厚度较大的翼型升阻比越大。雷诺数变大时，在小攻角下，翼型的厚度对升阻比的影响不敏感，但随着攻角的增大，厚度较大的翼型首先出现失速现象。

对包含塔架的三叶片风力机整机模拟表明，塔架本身能承受比较大比例的风阻力载荷(此处为30%)，但塔架对风力机叶片的转运力矩的气动影响很小，可以忽略不计。而通过模拟不同叶尖速比下的风力机运行效率，发现在叶尖速比为5左右时，其有最大功率系数。

通过模拟双风力机在不同前后间隔下的空气流场，发现在风力机场中，前后风力机间有很强的气动干涉。间隔越短，干涉越强。在间隔为5倍叶片直径条件下，空气下游风力机的功率只有上游风力机功率的一半。因而，风力机就尽可能避免处于上游风力机的近尾流区。

本研究的最主要意义是提出一种用来模拟风力机各种气动作用的统一框架模型，可以简单地处理叶片几何构型(拱曲度，厚度，扭转角，锥度等)、部件(塔架，轮毂，机舱等)、以及多风力机相互气动干涉等众多参数。在这种简单的统一框架下，有可能实现关于风力机气动性能的自动优化选型--这是正在进行的深入研究工作。