轴流风机被广泛应用在工业、消防以及各种公共场所的强制对流换热和通风换气中．为了实现风机全压升和效率等设计指标，须要根据实际工况对风机叶片的叶型进行气动优化设计．由于叶片内部三维流场的复杂性，因此其物理和数学上的简化模型二维基元叶栅成为人们研究分析叶轮机械内部流动情况的基础．在来流工况和安装角等条件相同的情况下，基元叶栅进出口速度三角形由基元叶栅型线决定[1]．当总压超过一定的范围时，叶栅流道内的逆压梯度会加剧，叶片吸力面上发生流动分离，从而导致流动损失迅速增大．这表明风机的总压和效率这两个气动性能参数在一定范围内是相互矛盾的，采用单目标优化设计方法难以获得同时具有高效率和高压升的叶片．传统的叶型优化主要利用效率、总压损失、熵参数或熵增等参数[2-3]来衡量叶轮机械流动损失大小，但这些参数只能反映损失总体累积的大小，不能反映流道内局部损失的具体来源和增加过程．熵产理论基于热力学第二定律可作为描述不可逆过程的量度．熵产是由不可逆过程产生的，它反映了系统的不可逆性及流动中能量损失的大小，进而反映了过程中能量转换与利用的完善程度，因此在不可逆现象比较集中的叶轮机械领域得到了广泛深入的研究．文献[4-5]通过引入损失功(irreversibility)即㶲损的概念，并基于压气机内流动为绝热、定常、稳态假设，对压气机内局部区域的损失进行了分析．文献[6]推导了含耗散函数项的总压输运方程，利用总压损失、熵增、耗散函数等参数分析某扩压叶栅流动与损失特性，探讨其损失产生机理．文献[7]推导出能量耗散率的分解式，并结合总压损失分析了耗散各组分在前缘损失、叶表损失和通道损失中的主导因素，但耗散函数和能量耗散率并不代表流动损失的全部，耗散函数只能表征黏性及湍流耗散产生的损失，而能量耗散率只能表征由湍流引起的耗散．文献[8]利用控制体方法，以低速压气机为研究对象，进行熵产的量化研究，开展间隙近乎于零范围的内动叶端区流动损失机理研究，并分析了引起流动损失的流动结构及其损失大小和变化规律．文献[9]应用雷诺应力模型对壁面绝热的旋风分离器的流场进行数值分析，对模拟结果采用熵产分析法和㶲分析法计算分离器的㶲损，其结果证实了热力学第二定律研究旋风分离器能量损失的可行性．文献[10]同样将熵产方法应用于水轮机，突出了该方法在水轮机能量损失分析中定量分析与定位分析的优势．文献[11-12]将熵产引入气动优化设计中，进一步表明熵产理论具有一定的工程应用前景．本研究从热力学第二定律及流体力学基本方程出发，在㶲分析方法的框架提出了一种基于熵产的基元叶栅多目标气动性能优化．通过建立叶栅参数化、网格生成、流场数值计算和数值优化方法一体化的平面叶栅优化平台，并以特定工况下最大全压升和最小熵产为优化目标，降低流场熵产损失，进而实现提高叶栅性能的优化目标．在总压提高的条件下，减少了流场熵产率．1 湍流熵产数值计算选定熵产参数作为衡量流动损失的关键物理量进行研究，根据高乌-史多台拉(Guoy-Stodola)公式EX,des=T0(Sgen)．(1)可见开口系统的各种内、外㶲损失的总和EX,des与环境绝对温度T0与系统熵产率Sgen有关[13]，因此计算熵产可以得到系统由不可逆过程引起的能量损失．根据能量方程可得熵的平衡方程为[13]                ∂(ρs)/∂t+div(ρUs)=-T-1div(q)+Φ/T-T-2q·∇T，(2)式中：ρ为密度；U=ui+vj+wk为速度矢量；T为温度；q为导热流密度矢量；s为比熵；Φ=∇·τ·U-∇·τ·U为耗散函数，τ为黏性应力张量．又根据热力学第二定律，开口系统的熵增等于熵产与熵流之和，因此开口系统单位体积熵产率为                s˙gen=∂(ρs)/∂t+div(ρsU)+T-1div(q)-T-2q·∇T，(3)式中：-div(ρsU)为质熵流；-T-1div(q)+T-2q·∇T为热熵流．联立方程(2)和(3)可得s˙gen=-T-2q·∇T+Φ/T．(4)式(4)表明：引起流场熵产的因素是有限温差引起的传热和黏性引起的耗散，其大小取决于流场的温度梯度场和速度梯度场．当湍流流动时，采用RANS方法并引入涡黏性假设和傅里叶定律，有s˙gen¯=[(γ+γt)/T¯2](∂T¯/∂xj)2+[(μ+μt)/T¯]τji¯(∂ui¯/∂xj)，(5)式中湍流黏性系数μt和湍流导热系数数γt是空间坐标函数，取决于流动状态而不是物性参数．对式(5)进行体积分即可知流场总熵产率为S˙gen¯=∫Ωs˙gen¯dΩ．(6)在对平面叶栅进行数值模拟之后，可以求出任意节点的熵产率以及任意区域的熵产率总和，并将熵产率作为后处理模式计算其分布来评价叶栅的流动损失．2 流场数值计算及优化方法2.1　数值计算方法校核二维叶栅几何模型的参数取自一种单转子低压轴流风机的截面参数．此轴流风机设计转速为 2 800 r/min，直径为0.3 m，轮毂比为0.35，叶片数均为6，设计全压升pt为180 Pa，设计流量为3 000 m3/h．其中叶片叶型是NACA6510，二维模型截面在风机相对叶高10%处，叶片安装角为55°，叶片弦长为80 mm．由于叶栅是沿周向均匀分布的，每个叶栅通道内的流动情况也是相同的，因此计算模型只要取部分叶栅通道即可．计算区域的几何模型与网格如图1所示．将整个叶栅简化为两流道进行模拟计算，叶栅进口和出口延伸的长度均为3倍弦长，计算区域分为三部分：进口域、转子域和出口域．采用非结构化网格，在叶片近壁面处进行网格加密．10.13245/j.hust.210910.F001图1计算域网格示意图计算湍流模型求解雷诺时均N-S方程，对于不可压缩流场，选择压力基求解器，设置速度进口，压力出口，叶片壁面采用无滑移边界条件，上下边界采用周期性边界条件．速度与压力耦合采用SIMPLEC算法求解，单元中心梯度则通过格林-高斯公式得到．为了充分捕捉近壁区信息，在近壁区必须细化网格并采用适应低雷诺数的SST k-ω湍流模型．当计算残差值小于1.0×10-6时认为达到收敛．经验证，在保证计算精度的前提下，选用约1.32×105网格可以保证数值计算结果与网格数的无关性，此时近壁面第一层网格y+值小于5，其边界层网格共布置8层，第一层边界层网格高度取1×10-4倍弦长．2.2　多目标优化算法采用改进的NSGA-II多目标遗传算法完成叶片气动优化设计．在采用实数编码的交叉操作中，传统的NSGA-II算法采用模拟二进制交叉(simulated binary crossover，SBX)算子，对于父代染色体x1和x2，子代染色体c1和c2为：c1=[(x1+x2)+β(x1-x2)]/2;c2=[(x1+x2)-β(x1-x2)]/2,式中β为随机变量，其表达式为β=(2u)1/(η+1)    (u≤0.5);[2(1-u)]-1/(η+1)     (u0.5),其中u为区间(0，1)上均匀分布的随机数．因为该算子搜索范围有限，使得算法易产生局部最优、进化过程不稳定等问题．在SBX算子中引入正态分布，得到正态分布交叉(NDX，normal distribution crossover)算子，此时随机变量β=1.481×N(0,1)，其中N(0,1)为正态分布随机变量．传统NSGA-II算法采用拥挤算子改善了种群的多样性，但同时由于其密度信息的估计仅仅限于同一级非支配个体集中，因此不能很好地反映个体周围的密度信息．采用动态拥挤算子DDCi=DCi/logVi-1，式中Vi=1Nobj∑m=1Nobjfi+1m-fi-1m-DCi2为个体i在各维目标相邻个体聚集距离的方差，能反映各维目标聚集距离差异程度，即不同维目标上拥挤距离的差异程度较大的个体在种群维护中有更多的机会被保留，其中，Nobj为优化目标数，fi+1m和fi-1m分别为第i+1和i-1个体的第m个目标值，DCi为原NSGA2第i个体的拥挤距离，即DCi=∑m=1Nobjfi+1m-fi-1mfmaxm-fminm，式中fmaxm和fminm为第m个目标值的最大和最小值．由于引入了正态分布交叉算子和动态拥挤算子，因此可以提高其对搜索空间的探索与开发能力，增强算法的鲁棒性，对多目标问题而言更容易得到完整的Pareto解集．通过每一层非支配个体的状态进行动态调节，对每一层个体分布密集的区域进行修剪，具有一定的自适应性，分布性好的个体有更大的生存机会．3 分析与讨论A．叶栅几何参数对性能的影响为了研究几何参数对叶栅性能的影响，分别计算了不同安装角、翼型厚度和叶栅稠度下叶栅的气动性能．采用如下4位数NACA翼型的中弧线形式y=f(x/d2)(2d-x/c)     (0≤xdc);f[(c-x)/(1-d)2](1+x/c-2d)    (pcx≤c),式中：c为叶型弦长；f为最大弯度；d为最大弯度位置．叶栅的翼型d=0.3，f=0.06．厚度分布采用如下形式：y=5t(0.296  90x-0.126  00x-0.351   69x2+0.284   30x3-0.101   50x4),式中t为最大相对厚度．图2(a)和(b)为不同翼型厚度的叶栅气动性能随安装角β的变化，此时叶栅稠度τ=1.0．由图可知：随着安装角的增大，叶栅的全压先增大后减小，而叶栅流场的熵产率先减小后增大．叶栅的全压随着翼型的厚度的增大而减小，在大的安装角下，厚度对翼型全压的影响逐渐变小；在小的安装角下，翼型越薄，对应的流场熵产率更低，随着安装角的增大，趋势却相反，翼型越厚，对应的流场熵产率反而更低．10.13245/j.hust.210910.F002图2不同厚度下叶栅性能随安装角变化图3(a)给出了两种厚度(t=10%，20%)在不同安装角(β=50°，70°)下的相对流线及压力云图．当安装角β=50°时，对于薄翼型(t=10%)，叶型吸力面的低压区较小，而且前缘处的高压区也较小，尾缘处未发生分离；随着厚度的增大，吸力面的低压区逐渐扩大，导致尾缘处发生流动分离．随着叶型的安装角增大，进气攻角增大，因此尾缘发生分离，而且随着安装角的增大，分离点前移．当β=70°时，对于薄翼型(t=10%)，叶型前缘压力面侧发生流动分离，因而导致更大的流动损失．图3(b)分别给出了对应的熵产云图以明确流动机制和损失机理．由于熵产的量级范围为1×10-10~1×104，因此对其进行对数处理以清楚显示损失来源．由熵产云图可见：当安装角β=50°时，叶型越厚，由于湍尾流场的损失越大，导致整个流场的熵产也越大；而当安装角β=70°时，叶型越厚，反而吸力面的流动损失更小，导致整个流场的熵产也越小．10.13245/j.hust.210910.F003图3不同厚度叶栅在各安装角下的压力云图和熵产云图图4(a)和(b)为不同翼型厚度的叶栅气动性能随叶栅稠度的变化，此时选取的叶型安装角β=55°．由图可知：随着叶栅稠密度的增大，叶栅的全压先增大后减小，而叶栅流场的熵产率也先减少后增大．10.13245/j.hust.210910.F004图4不同厚度下叶栅性能随叶栅稠度变化图5(a)给出了两种厚度(t=10%，20%)在不同叶栅稠度下(τ=2.0，1.0)下的相对流线及压力云图．由图可知：当栅距过大时(τ=1.0)，翼型尾缘发生流动分离，翼型厚度越大，分离越严重．而当栅距过小时(τ=2.0)，由于流道变窄，因此吸力面的低压区与压力面的高压区产生干涉，若此时翼型厚度过大，也会导致尾缘处的流动分离．选择合适的叶栅稠度可以较好地控制叶栅的流动分离．图5(b)分别给出了对应的熵产云图，可以看出：栅距过大时，尾流区的熵产明显变大，而栅距过小时，叶栅通道内的压力梯度变大，速度梯度也变大，导致流动损失增大．10.13245/j.hust.210910.F005图5不同厚度叶栅在各稠度下的压力云图和熵产云图 B．翼型截面的优化 前面研究了给定来流下叶型厚度、叶栅稠度以及安装角对叶栅气动性能的影响，这里在安装角β=55°、叶栅稠度τ=1.6的情况下对叶栅的叶型进行优化．为了在优化过程中灵活调整叶型型线，并且尽量减少型线参数变动的自由度，弧线与厚度分布均采用如图6所示的4控制点的3次B样条曲线．从图6(a)可以看出：B样条曲线具有良好的切矢性，曲线在P0点与线段P1P0相切，在P3点与线段P2P3相切，利用该特性，只须调整P1，P2的位置就可以使得χ1为前缘方向角，χ2为后缘方向角．从图6(b)可以看出：当线段P1P2与线段P0P3平行时，两线段的距离为B样条最大厚度的4/3．因此叶型的中弧线的参数化需要3个参数表示，即前缘方向角χ1、后缘方向角χ2、最大相对弯度f；同样厚度分布曲线也需要三个参数表示，即进口构造角ξ1，出口构造角ξ2，最大相对厚度t．10.13245/j.hust.210910.F006图6叶型参数化构造叶型优化目标是在设计进口来流下，以全压pt最大和熵产率S˙gen¯最小为优化目标，以上述6个参数为优化变量．其约束条件分别为15°≤χ1,χ2≤55°，15°≤ξ1,ξ2≤75°，0.05≤f≤0.20 ，0.06≤t≤0.24．优化算法参数设置如下：种群规模为40，交叉概率为0.7，变异概率为0.12，变量采用实数编码．优化过程中的种群分布如图7所示．由图7可知：优化初始种群分布杂乱，会出现全压小、熵产率大的叶型，随着优化的进行，通过交叉、变异，淘汰了部分劣势个体，种群的分布朝好的方向进化，但仍呈现不规律性．第10代后有收敛的趋势，在以后不断优化出全压更大、熵产率更小的个体，说明优化的有效性．到第20代，优化基本收敛．选择第30代的优化种群为最终的优化结果．10.13245/j.hust.210910.F007图7优化过程中的种群分布图8为优化后叶型目标函数的Pareto分布，图中每个数据点代表一种可行的优化翼型，它是全压和熵产率不断调和的结果，可以根据设计需要选择合适的最优叶型．由于优化的目标是最大全压和最小熵产率，因此由Pareto分布图可以看出：叶型的全压越大，对应流场的熵产率越大．取Pareto解集中任意三个解，所对应的优化翼型按分别命名为Opt_1，Opt_2和Opt_3，其中Opt_1叶型对应的全压最大，Opt_2叶型对应的全压最小，Opt_3对应一种中间性能叶型．10.13245/j.hust.210910.F008图8三种优化叶型流场压力云图与熵产云图图9(a)和(b)分别给出了三种优化翼型的中弧线以及厚度分布．可以看出：三种优化翼型的厚度分布近乎一致，趋于设计变量给定的最小厚度，因此在设计工况下，降低翼型厚度可以很好改善叶栅性能．三种最优叶栅随着总压的增大，对应的叶型弯度也随之增大．10.13245/j.hust.210910.F009图9三种优化叶型中弧线以及厚度分布图8(a)和(b)给出三种优化叶型流场的压力云图和熵产云图．可以看出：随着全压的增大，三种优化叶型吸力面的低压区明显增大，尾流区的熵产也随之增大．通过研究主要得出以下结论．a.定量计算了叶栅内熵产大小及来源，叶栅内流动损失主要是由耗散引起的熵产．熵产能够被用来方便、定量地研究局部区域内由于不可逆过程引起的损失，从而能够更准确地衡量局部损失的大小，为风机的设计、优化和改型提供了一个有效工具．b.随着安装角的增大，叶栅的全压先增大后减小，而叶栅流场的熵产率先减小后增大．在小的安装角下，翼型越薄，对应的流场熵产率更低，随着安装角的增大，趋势却反过来，翼型越厚，对应的流场熵产率反而更低．c.当栅距过大时，尾缘发生流动分离，翼型厚度越大，分离越严重；而当栅距过小时，由于流道变窄，导致吸力面的低压区与压力面的高压区产生干涉，若翼型厚度过大时，也会导致尾缘处的流动分离．选择合适的叶栅稠度，可以较好地减小叶栅的流动损失．d.给出的二维叶栅优化算例相对简单，对于更为复杂的三维叶轮机叶片优化，只须要考虑参数化方法及CFD求解方法，便可将熵产方法推广到更为复杂的气动优化设计中．