自激振荡脉冲射流减阻特性的多目标优化设计

汪朝晖，冯亚楠，高全杰

（武汉科技大学机械自动化学院，湖北 武汉 430081）

1 引言

高速运动的流体会在壁面边界层产生强烈的摩擦阻力，是高速流体能量损耗的主要原因。研究表明：高速流射流形态对流道壁面剪切力有较大影响，并且调整自激振荡腔室结构参数能够改变射流在边界层中的拟序结构，可实现流体减阻效果［1-2］。

利用壁面振动引起的湍流减阻效应，控制射流形态，可有效减少流体在壁面处的阻力。文献［3］利用雷诺数在（9700～19200）的射流，同时对固体壁面进行激振研究，实现了（10～15）%的减阻效果。然后利用仪器对该壁面展翼方向振动时的湍流边界层变化进行详细探讨。结果表明：射流流体近壁区部分受到壁面振动的影响，其平均速度梯度降低，壁面平均剪切应力下降，从而发现壁面的振动可以有效改善壁面阻力。文献［4］利用DNS 数值研究发现，流向涡旋具备循环再生的现象。旋涡结构从主射流中获得能量，这一能量通过旋涡形成了近壁条纹波，同时受到近壁条纹波的影响，又会诱发产生新的旋涡结构。因此，他们认为形成湍流壁面摩擦阻力的关键因素是射流中的流向涡旋，而抑制射流中流向涡旋的产生就可以到达降阻的效果。文献［5］在研究柔性涂层减阻时，将壁面模型设置为多个弹性体，并建立了壁面随流体波振的环形系统，从而发现振动频率、射流流速对减阻效果的影响，其减阻率随着振动频率和射流流速的提高而变大。文献［6］人研究了简谐振动的壁面下的流体运动，利用直接数值模拟，将阻力变化周期分为三个区域进行讨论，提出振动引起的负展向涡的产生以及条带滑移加宽和流向涡减弱是壁面振动产生减阻的主要原因。高速流体管道内沿程阻力主要来源于近壁湍流边界层，边界层中的湍流拟序结构会导致管道内的流阻大幅增加，而这种拟序结构主要包括：条带结构和流向涡结构，条带结构的向上抬升和向下运动均会导致流阻的急剧增加，而流向涡结构的周期性破碎是导致条带结构抬升、下降的关键原因，因此，抑制湍流拟序结构的周期猝发现象是提高管道内射流减阻效果的关键因素。

自激振荡腔室在一定条件下能够产生脉冲射流，形成周期性的湍流拟序结构［7］。有效调整射流形态，抑制湍流结构的周期猝发现象，可以实现减阻。文献［8］研究了自激振荡腔室的压力和速度的时均分布特征，并依据这些特征划分了腔室流域，提出了自激振荡腔室的脉冲射流形态与自激振荡腔室结构参数密切相关。文献［9］通过建立LES模型对自激振荡腔室的流动特性进行研究，并设计实验验证了LES模型的准确性，对比了不同振幅下的射流脉冲性能，发现了自激振荡腔的腔室长径比对自激振荡射流形态有着重要影响。因此，利用自激振荡脉冲装置特殊的流体及壁面自振特性，将自激振荡脉冲理论和壁面振动减阻理论结合起来，通过优化自激振荡腔室的关键结构参数，来抑制湍流流体近壁区的拟序结构的周期猝发现象，能够有效减少管道流阻。

为了获得减阻效果最佳的自激振荡腔室流道结构参数，需要设计有效的多目标优化方法。文献［10］采用正交试验法建立初始样本点，由多目标正交试验矩阵分析得出各因素权重，通过BP神经网络与遗传算法进行协同优化，得出了最优结构参数。文献［11］采用拉丁式立方法设计实验生成初始样本点，近似响应面法拟合出近似模型，运用NSGA-II多目标优化，求解出Pareto最优解集，实现了对最小压降和最大换热的多目标优化；文献［12］采用田口设计，Kriging响应面法构建模型，NSGA-II对板翅式换热器组进行多目标优化，讨论了四因素的二阶交互作用对换热的影响。优化效果关键在于两点：近似模型和优化算法。由于自激振荡腔室的设计变量与射流减阻特性之间没有直接精确的数学关系式，CFD模拟与响应面拟合精度取决于计算成本，同时优化算法存在一定计算误差，加之优化结构的误差随优化过程不断累计。因此，探索更高优化精度的优化方法，研究针对实际问题的有效多目标优化方法极其关键。

鉴于此，为了行之有效的优化自激振荡腔室结构，达到自激振荡脉冲射流的最佳减阻效果。这里以自激振荡腔室脉冲射流为研究对象，采用大涡模拟的方法分析了不同结构参数下的自激振荡腔室在气体管道运输中的减阻效果，并通过面心复合设计（FCCD）方法与响应面法结合来建立近似模型，采用NSGA-II算法对自激振荡腔室无量纲结构参数进行多目标优化，得到最优自激振荡腔室结构，为自激振荡脉冲减阻增输装置设计提供理论基础。

2 自激振荡腔室助推涡减阻效应

流体在自激振荡腔室出口内高速运动的过程中，流体阻力的形成主要可归纳为两部分，一部分是流体与固体壁面之间的摩擦阻力，另一部分是流体之间具有相对运动产生的剪切应力。这些阻力会导致流体的自身动能产生耗散，而耗散的动能转化为流体内的热能。流体的粘性和惯性和固体不光滑表面的阻碍作用是导致这一现象发生的重要原因。因此，在分析流体流动阻力时，必须对流体的剪切应力以及壁面特性进行着重分析。流体在沿壁面流动时，靠近壁面的粘性底层处流体的速度为零。流体流速沿着垂直方向从零开始逐步增加，然后逐步减缓直至基本不变。故流体在管路中流动可以分为边界层区域和主流区域。其中，流动主要阻力都集中在边界层区域，该区域速度变化较大，而主流区域离壁面较远，速度基本上不变的区域，其流动阻力很小可以忽略。自激振荡腔室出流管道中单位面积上的剪切应力为［13］：

式中：τ—单位面积上的剪切应力；μ—内摩擦系数（黏滞系数）；du∕dy—速度梯度。其中，速度梯度表示速度沿垂直于速度方向y的变化率，同时它也是直角变形速度，称为剪切变形速度。因此，式（1）可以理解为：剪切应力与速度梯度成正比。剪切速度应变率可以用来表征剪切应力，同时可以根据模拟计算得到壁面剪切应力推导计算壁面的摩擦阻力［13］：

式中：τ—单位面积上的剪切应力；d—管道直径；l—管道长度。

在管道直径和管道长度一定的条件下，壁面摩擦阻力和壁面剪切应力呈正比关系，因此，该摩擦阻力也和速度梯度呈正比关系，这里利用速度梯度变化来表征流体之间的内摩擦力变化。

在经过自激振荡腔室后，这种基本稳定的速度分布发生改变。自激振荡腔室出流管道水平射流剪切层涡环，如图1所示。靠近壁面的反向助推涡是由剪切层分离后形成的。这些反向助推涡夹带下的高速主射流同固体壁面没有直接相接触，起到了类似空气轴承的作用，将高速流体同固体壁面进行分层。

图1 自激振荡室出流管道水平射流剪切层涡环Fig.1 Self-Excited Oscillation Chamber Outflow Pipe Horizontal Jet Shear Layer Vortex Ring

3 数值计算

3.1 控制方程

与其他湍流模型相比，大涡模拟能够更加精准的捕捉瞬态流动，因此可以采用大涡模拟有效求解自激振荡腔室周期性瞬时的剪切应力和速度梯度变化。大涡模拟采用滤波函数将湍流瞬时运动中尺度比滤波函数小的涡环滤掉，从而得到大涡流动的运动方程。在非定N-S方程的基础上进行Favre滤波后得到可压缩理想气体的连续性方程、动量方程、能量方程［7］。

3.2 有限元模型

3.2.1 物理模型及关键参数

自激振荡脉动特性主要由自激振荡腔室结构和关键参数决定［7］，如图2所示。关键结构参数如下：喷嘴进出口直管道管径d；上喷嘴入口流道管径d1；下喷嘴出口流道管径d2；圆锥形收缩管夹角α1；下喷嘴碰撞壁夹角α2；圆锥形扩散管夹角α3；自激振荡脉冲腔室长度L；自激振荡脉冲腔室直径DT。依据先前的实验和仿真结果，自激振荡腔室结构参数优选范围和关键参数［13］，如表1所示。

图2 自激振荡腔室结构Fig.2 Self-Oscillating Chamber Structure

表1 自激振荡腔室无量纲结构参数Tab.1 Self-Excited Oscillation Chamber Dimensionless Structural Parameters

3.2.2 模型验证及网格无关性

计算方法选用大涡模拟，计算模型为可压缩理想气体，时间项为二阶隐式格式，对流项和扩散项均为二阶迎风格式。模型边界条件，如图3所示。根据先前实验［16］，计算模型边界条件采用压力入口和压力出口，上游入口压力设为101325Pa，下游出口压力设为84300Pa。

图3 网格模型及边界条件Fig.3 Grid Model and Boundary Conditions

由于3D-LES 计算成本较高，为了有效提升运算效率，采用2D 轴对称大涡模拟是可行的［14-15］。通过与文献［16］的实验数据进行比较，测试了上游入口压力为90500Pa条件下的平均流量，测试结果，如表2所示。采用2D-LES模型的误差为1.41%，验证结果表明数值计算结果和实验结果之间具有良好的一致性。因此，该数值模型对于本研究是可靠的。

表2 2D-LES模型的平均流量误差Tab.2 Average Flow Error of the 2D-LES Model

分别使用116282网格、194300网格及415721 网格，对自激振荡腔室中心结构的2D-LES 计算模型进行网格无关性测试。轴向流速的网格无关性检测，如图4所示。综合权衡计算效率和求解精度，最终确定选用194300网格模型。

图4 网格无关性验证Fig.4 Grid Independence Verification

4 多目标优化

通过面心复合设计、响应面法和NSGA-II算法进行多目标协同优化，获取自激振荡腔室无量纲结构参数的最优解，具体的多目标优化流程图，如图5所示。

图5 自激振荡腔室多目标优化流程图Fig.5 Multi-Objective Optimization Flow Chart of Self-Excited Oscillation Chamber

4.1 约束处理

自激振荡腔室减阻特性的关键要求是降低流阻，其中流阻主要包含两部分，一部分是流体与固体壁面之间的摩擦阻力，另一部分是流体之间相对运动的内摩擦。因此，出口壁面平均摩擦阻力F1和出口流道平均内摩擦力F2作为多目标优化的两个目标，计算公式如下：

式中：τ—单位面积上的剪切应力；A—自激振荡腔室出口流道网格总面积；μ—内摩擦系数（黏滞系数）；du∕dy—速度梯度；V—自激振荡腔室出口流道网格总体积。

4.2 试验设计

中心复合设计（CCD）和面心复合设计（FCCD）是响应面分析的常用试验设计方法之一。试验设计的实验点按照类别可以分为中心点、立方点及轴向点三种，不同类别实验点位置示意图，如图6所示。

图6 CCD和FCCD设计点分布对比Fig.6 Comparison of CCD and FCCD Design Points

可以看出在CCD 的设计点中，存在超出原定水平的数据。由于样本点数据超出原定水平范围违背了实际工况要求，因此，面心复合设计（FCCD）在本次优化中被采用，既保留了CCD试验设计在拟合各个因素与输出指标间的响应关系的优势，又避免了样本点数据超出原定水平。设计各因素水平及编码值，如表3所示。

表3 试验设计因素水平Tab.3 Trial Design Factor Level

4.3 响应面近似模型

响应面表面形式族的选择对响应面的分析结果具有重大影响，因此，所选择的响应面表面形式应满足特定光滑度要求。对于不同的实际应用，最优的响应面表面形式的选择也不相同。由于多项式的计算简单，并且能够生成封闭形式的表达式。在计算流体力学中，多项式成为表示响应面的常用形式。鉴于多目标优化的设计变量数为3，选择二次多项式响应面的形式，其表达式为［12］：

式中：β0，βi，βij—回归系数。响应面y表示为设计变量在设计空间内的函数。

在提出响应面近似模型用于多目标优化之前，应检测响应面模型是否能够满足精度要求。如果不满足精度要求，则需要增加样本点，从而调整响应面的可变性和偏差之间的平衡，使得偏差减小，响应面可变性增强。R2准则和修正复相关系数R2adj通常被用于响应面的验证，其表达式如下［12］：

式中：SSR—指回归平方和，表示因回归方程所引起的y的不均匀程度；SSY—指总体平方和，表示观测值y的不均匀程度；SSE—误差的平方和，表示因随机误差引起的y的不均匀程度；k—优化变量个数；n—样本点个数。

4.4 多目标优化算法

在NSGA-II中，优化的最终结果是一组Pareto解集，也被称为Pareto前沿。SBX作为交叉和突变的操作机制被用于生成下一代个体，包含交叉操作和变异操作，其表达式分别为［12］：

交叉操作：

变异操作：

式中：xi(1，t)和xi(2，t)—父代个体；xi(1，t+1)和xi(2，t+2)—子代个体；βqi—随机变量；xiUB和xiLB—变量上下限；δq—变异算子；u—［0，1］之间的随机数；ηm—指定的分布指数；δ—干扰因子。多目标优化算法的设置，如表4所示。

表4 NSGA-II优化算法参数设置Tab.4 NSGA-II Optimization Algorithm Parameter Settings

5 优化结果与分析

5.1 响应面模型构建

FCCD设置因素数为3，由于CFD模拟过程中不存在实验环境等其他不可控因素影响，因此中心试验重复次数为5，因变量个数为2，分别为出口壁面平均摩擦阻力F1和出口流道平均内摩擦力F2。为了消除两个目标之间量纲的影响，对CFD 数值结果归一化处理，试验设计样本点集及数值模拟所得结果，如表5所示。

表5 试验设计和CFD归一化结果Tab.5 Experiment Design and CFD Normalized Results

这里采用Quadratic 模型对目标变量F1，F2和设计变量d2∕d1，DT∕d1，L∕DT之间的函数关系进行拟合，将设计变量d2∕d1，DT∕d1，L∕DT分别记为x1、x2、x3，得到响应面回归方程，如式（10）、式（11）所示：

5.2 模型显著性分析及误差检验

通过对响应面回归模型进行检测，方差分析结果，如表6所示。Prob＞F表示无显著影响的概率，并且F值越大P值越小说明模型的相关性越显著。当P小于0.01时，表示模型非常显著。当P值大于0.01且小于0.05时，模型为显著。当P值大于0.05表示模型不显著。

表6 响应面回归模型的方差分析结果Tab.6 Analysis of Variance by Response Surface Regression Model

对响应面回归模型进行误差统计分析，相关系数R2值和Radj2值越接近1，说明相关性越好。一般情况下，实际工程应用要求相关系数R2值和Radj2值大于0.9。响应面回归模型误差统计分析，如表7所示。统计分析结果表明，这里所得的响应面近似模型符合上述检验原则，具有较好的适应性。同时通过图7可以看出，预测值与实际值基本分布在一条直线上，说明回归模型具有较高的精度。

表7 响应面回归模型误差统计分析表7 Statistical Analysis of Errors in Response Surface Regression Models

图7 预测值与实际值分布Fig.7 Distribution of Predicted and Actual Values

5.3 响应面关系及Pareto解集分析

响应面法能够通过样本点拟合出直观的三维响应面，并且可以通过等高线图的特征反映出设计变量之间的交互作用。通过控制某个因素不变的条件下，得到另外两个因素的交互作用对出口壁面平均摩擦阻力F1和出口流道平均内摩擦力F2的影响。其中，等高线为椭圆形，表示交互作用显著。d2∕d1，DT∕d1和L∕DT对出口壁面平均摩擦阻力F1的响应关系，如图8 所示。可以看出，当L∕DT处于0 水平时，出口壁面平均摩擦阻力F1随d2∕d1的增大而减小，随着DT∕d1的增大而减小，并且无明显交互作用。当DT∕d1处于0 水平时，出口壁面平均摩擦阻力F1随着d2∕d1的增大而减小，并且随L∕DT的增大出口壁面平均摩擦阻力F1先增大后减小，并且无明显交互作用。当d2∕d1处于0 水平时，出口壁面平均摩擦阻力F1随着DT∕d1的增大而减小，并且随L∕DT的增大出口壁面平均摩擦阻力F1先增大后减小，并且无明显交互作用。无论DT∕d1，L∕DT处在什么水平，F1随着d2∕d1的增大都呈现出了明显的降低趋势，并且与DT∕d1，L∕DT相比d2∕d1起主导作用。从等高线图可以看出并非椭圆形说明d2∕d1与DT∕d1，L∕DT对于出口壁面平均摩擦阻力F1的响应无明显交互作用。d2∕d1，DT∕d1和L∕DT对出口流道平均内摩擦阻力F2的响应关系，如图9所示。

图8 d2∕d1，DT∕d1，L∕DT对出口壁面平均摩擦阻力的响应关系Fig.8 Response Relationship of d2∕d1，DT∕d1，L∕DTto the Average Frictional Resistance of the Outlet Wall

图9 d2∕d1，DT∕d1，L∕DT对出口流道平均内摩擦阻力的响应关系Fig.9 Response Relationship of d2∕d1，DT∕d1，L∕DTto the Average Internal Frictional Resistance of the Outlet Flow Channel

由图9可以看出，当L∕DT处于0水平时，出口流道平均内摩擦阻力F2随d2∕d1的增大而先增加后减小，随着DT∕d1的增大而增大，并且无明显交互作用。当DT∕d1处于0水平时，出口流道平均内摩擦阻力F2随着d2∕d1的增大而先增大后减小，随L∕DT的增大而先增大后减小，并且从等高线投影结果可以看出有明显交互作用，出口流道平均内摩擦阻力F2存在极大值点。当d2∕d1处于0水平时，出口流道平均内摩擦阻力F2随着DT∕d1的增大而增大，并且随L∕DT的增大出口流道平均内摩擦阻力F2先增大后减小，并且无明显交互作用。

NSGA-II 算法结合响应面模型进行多目标优化计算，得到Pareto最优前沿，如图10所示。由于在Pareto前沿中，没有任何一个解能够优于其他解，因此任何一个解都是可以被选择Pareto最优解。

图10 目标F1和F2的Pareto最优前沿Fig.10 Pareto Front for Objectives F1 and F2

5.4 中心结构与优化结构数值模拟对比

通过多目标优化，充分考虑加工精度和设计标准，引入1：1的决策权重，从Pareto前沿中选出1个符合设计条件的候选设计点进行CFD仿真，并与中心设计点原始结构参数进行对比，结果如表8所示。

表8 中心结构与优化结构对比Tab.8 Comparison of Central Structure and Optimized Structure

从优化结果可以得出，优化后的模型与CFD 模拟值之间存在一定的误差，分别为2.18% 和4.58%，满足精度要求。与中心原始结构相比，出口壁面平均摩擦阻力和出口流道平均内摩擦力分别减少了43.01%和54.52%。

6 结论

这里对自激振荡腔室的无量纲结构参数进行了多目标优化，以获得最优减阻效果。采用了2D大涡模拟的数值分析方法，结合面心复合设计（FCCD）、响应面法和NSGA-II算法进行了多目标协同优化。分析了不同无量纲结构参数对于自激振荡腔室减阻特性的影响。得到如下结论：

（1）通过CFD与面心复合设计（FCCD）、响应面法和NSGA-II算法的多目标协同优化方法，能够有效获取自激振荡腔室无量纲结构参数的最优解，实现自激振荡腔室的最佳减阻。

（2）分析响应面关系，得出了d2∕d1，DT∕d1，L∕DT与出口壁面平均摩擦阻力F1和出口流道平均内摩擦力F2之间的关系。出口壁面平均摩擦阻力F1随着d2∕d1的增大都呈现出了显著降低，并且与DT∕d1，L∕DT相比d2∕d1起主导作用，d2∕d1，DT∕d1，L∕DT无明显交互作用。出口流道平均内摩擦阻力F2s随着DT∕d1的增大都呈现出了显著升高，并且与d2∕d1，L∕DT相比DT∕d1起主导作用，d2∕d1与L∕DT之间存在明显交互作用，并存在F2极大值点。

（3）充分考虑加工精度和设计标准，引入1：1的决策权重，从Pareto前沿中选出1个符合设计条件的候选设计点作为优化后的自激振荡腔室无量纲结构。与中心原始结构进行对比，结果表明，优化后的自激振荡腔室减阻性能得到了显著提高。与中心原始结构相比，出口壁面平均摩擦阻力和出口流道平均内摩擦力分别减少了43.01%和54.52%。