CN110610050B - 基于改进径向基函数变形算法的翼型气动减阻方法 - Google Patents

Abstract

本发明公开了一种基于改进径向基函数变形算法的翼型气动减阻方法，目的是提高翼型减阻设计的自动化程度并进一步降低翼型的风阻系数并提高减阻优化效率。技术方案是将翼型几何边界上的所有节点作为设计变量，建立关于阻力变量d的翼型气动阻力函数，采用序列最小二乘规划算法求解气动阻力函数f(d,S)关于设计变量组的最小值问题，得到网格内边界上N_w个内边界点最佳的位移量序列，采用增量式求解的RBF网格变形算法对翼型网格进行气动减阻变形，得到最佳的风阻系数或在最大迭代步数内翼型减阻达到最好结果的翼型网格和翼型几何。本发明可适用于多种类型的网格，可显著提高翼型减阻设计的自动化程度，且减阻优化效率更高，结果更为精确。

Description

基于改进径向基函数变形算法的翼型气动减阻方法

技术领域

本发明属于航空飞行器气动减阻技术领域，具体涉及一种基于改进径向基函数(Radial Basis Function，RBF)变形算法的翼型气动减阻方法。

背景技术

在飞行器的气动减阻设计过程中，机翼作为承受升力和阻力的主要部件，机翼剖面的翼型是气动设计中非常重要的环节，基本翼型的设计关系着机翼和机身的压力分布，从而对飞行器的流动形态有着重要影响。传统的飞行器及其翼型的设计方法多以风洞试验或计算流体力学(Computational Fluid Dynamics，CFD)分析结果为基础，采用依赖经验的人工试凑和反复修正为主。随着科学技术的发展，对飞行器性能指标的要求越来越高，因此亟需高自动化、高效率的减阻方法。好的气动减阻方法可以提高翼型设计质量并降低风阻系数，同时缩短设计周期；另一方面，随着对翼型性能要求的增加，设计过程中所面临的约束和矛盾越来越多，利用自动化的减阻技术有助于在这些复杂约束条件下得到良好的设计结果，并进一步降低翼型的风阻系数。

近年来，基于CFD数值模拟的气动外形自动化减阻技术快速发展，并且在航空领域得到了广泛应用。大体而言，气动外形自动化减阻技术主要分为两类：一类是随机型减阻优化方法，主要基于遗传算法和模拟退火算法等[郄永军等，基于多岛遗传-模拟退火算法的气动外形优化，飞行力学，2017年第5期]；另一类是梯度型减阻优化方法，基本思想是利用目标函数关于设计变量的梯度作为在整个可行域内确定寻优方向的标准，然后采用最小二乘法等经典算法进行寻优，其中梯度计算是该方法的关键。其典型代表是基于控制理论的伴随优化方法[Michael B.Giles等，An Introduction to the Adjoint Approach toDesign，Flow,Turbulence and Combustion，Vol.65，2000；Michael B.Giles等，伴随设计方法简介，流动、湍流与燃烧，2000年第65卷]。伴随优化方法的特点是将设计问题转化为寻找满足约束的最优控制问题，将流场方程作为约束条件引入目标函数中，控制函数则是物体的边界形状。梯度的求解需要分别计算流场方程和伴随方程，而伴随方程是由目标函数变分得到，与设计变量的多少几乎没有关系，所以求解梯度的计算量只是一般流场计算量的两倍，与其他计算方法相比，大大减小了计算量。

决定气动设计结果的主要因素包括CFD数值模拟、优化方法以及网格变形(meshdeformation)算法等，这些因素之间相互影响相互制约。网格变形算法是在保持网格拓扑结构和节点总数不变的情况下，只将网格节点移动到新的位置，以符合计算区域的变化。在翼型气动优化设计中，当翼型几何改变后，网格变形算法将内部网格节点重新移动调整以获得新的高质量网格，这对CFD数值模拟结果的准确性至关重要。该类算法能够保持网格关联信息不变，避免每个时间步的流场变量插值；此外还能处理任意单元类型的网格，因此也特别适合单元类型多变的非结构网格。

基于RBF插值的网格变形算法，由于精度高、适用性强，因而在动网格问题中得到了广泛应用。给定基函数φ()，RBF变形算法中坐标为X的任意网格点在第n坐标分量上的位移可以通过位移函数s(X)表示：

其中，n表示二维(三维问题是三维，下同)坐标方向中的任意一个分量，

表示第p个边界点的二维坐标，φ()是以欧式距离

为自变量的径向基函数，N_b是网格边界点的数量，λ_p是关于第p个边界点的权重系数。在这N_b个边界点中，动边界点的位移可在气动减阻应用中由优化算法计算得到，而静边界点的位移始终为0。因此，可以将这N_b个边界点的位移ΔX_b当作已知条件，反过来求取位移函数中的权重系数。ΔX_b是由这N_b个边界点的位移构成的维度为N_b×1的矩阵，如公式(2)所示：

通过求解公式(3)，可以计算得到维度为N_b×1的权重系数矩阵λ：

其中，Φ_b,b ^-1是Φ_b,b的逆矩阵，Φ_b,b是维度为N_b×N_b的距离函数矩阵，包含了这N_b个边界点之间的RBF插值：

Φ_b,b(p,q)表示矩阵Φ_b,b第p行第q列的元素，

是第p个边界点坐标，

是第q个边界点坐标，1≤p,q≤N_b。获得权重系数后，网格第l个内部点的空间坐标

在第n坐标分量上的位移

可通过公式(5)计算得到：

径向基函数φ()一般可以分为两类：全局型函数和紧凑型函数。紧凑型函数的控制点只作用在支撑半径内的网格点，即只有在半径范围内的网格点位移需要计算该控制点的影响。因此，紧凑型函数可以产生一个稀疏的插值矩阵，计算效率较高。全局型函数则是产生一个全域插值矩阵，耗时较大，且不同径向基函数的网格变形效果也相差甚远，实际中需要选择变形效果好、鲁棒性强的基函数。当网格规模较大时，RBF变形算法会产生规模巨大的线性方程组，导致整体计算效率的严重下降。为了减轻这一压力，研究者提出了一种基于最小化误差函数的数据减缩算法[T.C.S.Rendall等，Efficient Mesh Motion UsingRadial Basis Functions with Data Reduction Algorithms，Journal ofComputational Physics，Vol.228，2009；T.C.S.Rendall等，基于数据缩减的径向基函数高效网格运动算法，计算物理杂志，2009年第228卷]。该算法只选取部分边界点作为控制点，代替选取所有边界点来计算RBF插值函数，此方法已被证明能够大幅度降低计算时间和存储方面的成本。虽然RBF变形算法已得到较广泛的使用，但其还有进一步优化的空间，而且由于RBF方法涉及参数众多，目前还缺乏系统性的参数配置策略指导实际应用。此外，目前也未见有基于RBF变形算法的翼型减阻优化方法的公开报道。

发明内容

为克服现有技术的上述缺点，本发明要解决的技术问题是：提供一种基于改进径向基函数变形算法的翼型气动减阻方法，通过改进RBF网格变形算法(即变为增量式求解的RBF网格变形算法)，采用CFD数值模拟和伴随优化方法解决飞行器翼型设计中的自动化气动减阻问题，进一步降低设计翼型的风阻系数并提高减阻优化效率。

本发明技术方案如下：

第一步，读入初始翼型对应的二维网格文件，所述二维网格的内边界对应初始翼型S(即初始翼型的几何形状，如图1中由实心圆点围成的二维翼型)，远场外边界距离内边界20倍翼型弦长(翼型前后缘的距离，如图1中的虚线所示)以上；令初始翼型S的二维网格边界点的总数量为N_b，其中内边界点总数量为N_w，内边界点的坐标序列为

其中X_i是一个内边界点的二维矢量坐标(二维网格上的点的坐标均为二维矢量)，

分别为内边界点坐标X_i两个方向的坐标分量，1≤i≤N_w；令外边界点总数量为N_o，外边界点的坐标序列为

分别为外边界点坐标Y_j两个方向的坐标分量，1≤j≤N_o；令网格内部点的总数量为N_i，内部点的坐标序列为

分别为网格内部点坐标Z_k两个方向的坐标分量，1≤k≤N_i。

第二步，建立关于阻力变量d的翼型气动阻力函数

(f(d,S)为翼型S所受的总阻力，等于S上阻力变量d的积分，总阻力越小表明翼型风阻系数越小)；设置减阻收敛阈值ε(ε为浮点数且0＜ε＜1)和减阻过程的最大迭代步T_max；当迭代过程中翼型气动阻力函数值的变化比值(完全收敛时无变化，比值为1)大于减阻收敛阈值ε时，说明迭代收敛；最大迭代步T_max表示迭代超过T_max步后直接结束，用于控制计算总量和防止无法收敛的情况，为正整数，一般设置为100以内。

第三步，令当前迭代步T＝0，阻力函数旧值V0＝-1。

第四步，采用CFD求解软件(如SU2开源软件4.3.0以上版本中的SU2_CFD软件，https://su2code.github.io)计算流动控制方程R(W,S)＝0，得到流场变量W；根据翼型阻力公式g(W)[周伟等，翼型阻力计算方法的数值模拟研究，科学技术与工程，2011年第11期]计算阻力变量d＝g(W)，并将阻力变量d代入翼型气动阻力函数

得到阻力函数新值V，V＝f(d,S)。

第五步，若T≥1且V/V0＞ε，表明翼型的气动阻力函数值收敛到最小值，即翼型已达到最佳的风阻系数设计；或者T≥T_max，表明已达到最大迭代步数内翼型减阻的最好结果，这两种情况下转第十一步；否则转到第六步。

第六步，基于流场变量W，采用伴随方程求解软件(如SU2开源软件4.3.0以上版本中的SU2_CFD及SU2_DOT软件，https://su2code.github.io)计算连续伴随方程

(引入伴随变量ψ，使得气动阻力函数的导数只与翼型的几何变化相关)[陈帅等，基于连续伴随方程和FFD技术的气动外形优化设计，江苏航空，2018年第1期]，得到伴随变量ψ。

第七步，将S上的N_w个内边界点坐标作为设计变量组。

第八步，采用SLSQP(Sequential Least Squares Programming，序列最小二乘规划)算法(如SU2开源软件4.3.0以上版本中的SLSQP优化器，https://su2code.github.io)求解气动阻力函数f(d,S)关于设计变量组的最小值问题

(当翼型变形到使气动阻力函数的导数为最小值的位置时，由于该导数为负，变形后翼型的气动阻力值会减小，从而实现减阻优化)，得到网格内边界上N_w个内边界点最佳的位移量序列

其中ΔX_i为第i个内边界点最佳的位移量，

分别表示第i个内边界点在两个方向的位移分量，1≤i≤N_w。

第九步，根据翼型网格边界点的位移，采用增量式求解的RBF网格变形算法对翼型网格进行气动减阻变形：

9.1、由于远场外边界点始终静止不动，将N_o个外边界点的位移量(点的位移量是指该点相对于原来位置的偏移量)均设为(0,0)；令控制点坐标序列为空，令控制点位移量序列为空，令控制点数量N_c＝0；

9.2、从内边界点坐标序列

随机选择3个坐标，令这3个坐标对应的网格点为r₁、r₂、r₃，令第一控制点的坐标

令第二控制点的坐标

令第三控制点的坐标

将C₁、C₂、C₃加入控制点坐标序列，令第一控制点位移

令第二控制点位移

令第三控制点位移

将ΔC₁、ΔC₂、ΔC₃加入控制点位移量序列；

9.3、从外边界点坐标序列

随机选择3个坐标，令这3个坐标对应的网格点r₄、r₅、r₆，令第四控制点坐标

令第五控制点坐标

令第六控制点坐标

将C₄、C₅、C₆加入控制点序列，令第四控制点位移ΔC₄＝(0,0)，令第五控制点位移ΔC₅＝(0,0)，令第六控制点位移ΔC₆＝(0,0)，将ΔC₄、ΔC₅、ΔC₆加入控制点位移量序列；由此得到初始控制点坐标序列[C₁,C₂,C₃,C₄,C₅,C₆]以及对应的控制点位移量序列[ΔC₁,ΔC₂,ΔC₃,ΔC₄,ΔC₅,ΔC₆]，令N_c＝6，其中

分别表示第v个控制点在两个方向的位移分量，1≤v≤N_c；

9.4、采用公式(6)作为RBF网格变形算法的基函数：

其中ξ＝||x_p1-x_q1||/R(||x_p1-x_q1||表示坐标x_p1与坐标x_q1之间的欧式距离)，x_p1,x_q1泛指两个网格点坐标，支撑半径R表示网格点的影响范围(若两个网格点之间的欧氏距离大于R，则这两个点的移动互不影响)，支撑半径R设为翼型弦长的3-5倍；

9.5、根据控制点坐标序列中的边界点，由公式(7)计算得到维度为N_c×N_c的距离矩阵Φ_b,b中第p2行第q2列的元素Φ_b,b(p2,q2)，并令Φ_b,b(q2,p2)＝Φ_b,b(p2,q2)，其中1≤p2≤N_c,p2≤q2≤N_c；

Φ_b,b(p2,q2)＝φ(||C_p2-C_q2||) (7)

公式(7)中，||C_p2-C_q2||表示坐标C_p2与坐标C_q2之间的欧式距离；

9.6、令两个坐标维度方向上关于RBF函数权重系数的线性方程组的初始解分别为λ0¹＝[0,0,0,0,0,0]^T,λ0²＝[0,0,0,0,0,0]^T；

9.7、令坐标维度变量n＝1；

9.8、由控制点位移量序列中第n维度的位移分量构成维度为N_c×1的矩阵

其中1≤i1≤N_c，形成公式(8)的线性方程组：

9.9、以λ0ⁿ作为初始解，采用预条件投影的共轭梯度方法[Z.Dostál，Conjugategradient method with preconditioning by projector，International Journal ofComputer Mathematics,Vol.23，1988；Z.Dostál，预条件投影的共轭梯度方法，国际计算机数学杂志，1988年第23卷]求解公式(8)的线性方程组，得到维度为N_c×1的权重系数矩阵λⁿ；

9.10、由权重系数矩阵λⁿ得到第n维度如式(9)的RBF函数sⁿ(X)：

φ(||X-C_j1||)是以||X-C_j1||为自变量的基函数，||X-C_j1||表示坐标X与坐标C_j1之间的欧式距离，

是权重系数矩阵λⁿ的第j1个元素；

9.11、令n＝n+1，若n≤2，转到第9.8步；否则转到第9.12步；

9.12、令内边界点变量i＝1，令当前最大位移误差σ_now＝-1；

9.13、将内边界点坐标X_i分别代入RBF函数s¹(X)和s²(X)(s¹(X)和s²(X)是两个坐标方向上不同的函数，第9.10步得到的RBF函数sⁿ(X)，其中n＝1,2)，得到临时位移量(s¹(X_i),s²(X_i))；与X_i的实际位移量

进行比较，计算得到临时位移误差

若σ_now＜σ_tmp，令σ_now＝σ_tmp，令临时控制点坐标C_tmp＝X_i，令临时控制点位移ΔC_tmp＝ΔX_i，转9.14；若σ_now≥σ_tmp，直接转9.14；

9.14、令i＝i+1，若i≤N_w，转到第9.13步；否则转到第9.15步；

9.15、令外边界点变量j＝1；

9.16、将外边界点Y_j分别代入RBF函数s¹(X)和s²(X)，得到临时位移量(s¹(Y_j),s²(Y_j))；与点Y_j的实际位移量(0,0)进行比较，计算得到临时位移误差

若σ_now＜σ_tmp，令σ_now＝σ_tmp，令临时控制点C_tmp＝Y_j，令临时控制点位移ΔC_tmp＝(0,0)，转9.17；若σ_now≥σ_tmp，直接转9.17；

9.17、令j＝j+1，若j≤N_o，转到第9.16步；否则转到第9.18步；

9.18、设置位移误差阈值σ＝L/T_max，其中L是初始二维网格中边长的最小值，T_max是减阻过程的最大迭代步；

9.19、若当前最大位移误差σ_now＜位移误差阈值σ，表明由当前控制点序列计算得到的RBF函数满足精度要求，转到第9.23步计算内部点位移量；否则转到第9.20步更新控制点坐标序列；

9.20、将C_tmp加入控制点坐标序列，将ΔC_tmp加入控制点位移量序列，令N_c＝N_c+1；

9.21、在距离矩阵Φ_b,b的最后增加一行一列：根据控制点序列采用公式(7)计算得到距离矩阵Φ_b,b中第p3行第N_c列元素Φ_b,b(p3,N_c)的值，并令Φ_b,b(N_c,p3)＝Φ_b,b(p3,N_c)，其中1≤p3≤N_c；

9.22、令λ0¹＝[λ¹,0]^T,λ0²＝[λ²,0]^T分别作为两个坐标方向扩充后线性方程组的初解；转到第9.7步重新计算权重系数；

9.23、令内部节点变量k＝1；

9.24、将第k个内部节点的坐标Z_k分别代入RBF函数s¹(X)和s²(X)，得到该点位移量

并更新该点坐标Z_k，令

9.25、令k＝k+1，若k≤N_i，转到第9.24步；否则转到第9.26步；

9.26、令内边界点变量i＝1；

9.27、将第i个内边界点的坐标X_i分别代入RBF函数s¹(X)和s²(X)，得到该点位移量

更新第i个内边界点坐标X_i，令

9.28、令i＝i+1，若i≤N_w，转到第9.27步继续对翼型几何S上的内边界点进行更新；否则说明对翼型网格进行气动减阻变形处理完毕，得到气动减阻后的新翼型网格，转第十步；

第十步，令T＝T+1，V0＝V，转到第四步，对气动减阻后的新翼型网格继续进行翼型气动减阻。

第十一步，输出达到最佳的风阻系数或在最大迭代步数内翼型减阻达到最好结果的翼型网格和翼型几何，结束。

采用本发明可以达到以下技术效果：

1、本发明结合CFD数值模拟和连续伴随方法，采用增量式求解的RBF网格变形算法对翼型网格进行气动减阻变形，可显著提高翼型减阻设计的自动化程度并进一步降低设计翼型的气动阻力和风阻系数。

2、本发明将翼型几何边界上的所有节点作为设计变量，且在减阻变形过程中始终保持高质量的网格，使得减阻优化的计算结果更为精确。

3、本发明采用增量式求解的RBF网格变形算法变形翼型网格，该算法利用预条件投影的共轭梯度方法增量式求解RBF函数权重系数方程(公式8)，提高了减阻方法的算法效率。

4、本发明具有通用性，可适用于结构、非结构等多种类型的网格。

附图说明

图1是本发明所处理的二维翼型几何示意图，二维翼型是机翼z方向的垂直剖面，由图中实心圆点所围成，虚线是二维翼型的弦长。

图2是本发明第一步读入的初始翼型二维网格可视化之后的示意图，图2(a)是初始翼型S二维网格可视化之后的整体网格示意图；图2(b)是图2(a)所示网格中心区域(即中心全黑部分)的局部放大图。

图3是本发明总体流程图。

图4是本发明第九步的流程图。

具体实施方式

如图3所示，本发明的基于改进径向基函数变形算法的翼型气动减阻方法流程包括如下步骤：

第一步，读入初始翼型对应的二维网格文件，所述二维网格的内边界对应初始翼型S(即初始翼型的几何形状，如图1中由实心圆点围成的二维翼型)，远场外边界距离内边界20倍翼型弦长(翼型前后缘的距离，如图1中的虚线所示)以上，读入初始翼型S二维网格如图2所示，图2(a)是初始翼型S二维网格可视化之后的整体网格示意图，X，Y坐标轴不属于网格的一部分，这里用坐标仅用于说明网格的大小和位置，图2(a)中最外面的圆形是初始翼型网格的外边界；图2(b)是图2(a)所示网格中心区域(即中心全黑部分)的局部放大图，图2(b)中梭形是初始翼型网格的内边界(即翼型几何)；令二维网格边界点的总数量为N_b，其中内边界点总数量为N_w，内边界点的坐标序列为

其中X_i是一个内边界点的二维矢量坐标(二维网格上的点的坐标均为二维矢量)，

分别为内边界点坐标X_i两个方向的坐标分量，1≤i≤N_w；令外边界点总数量为N_o，外边界点的坐标序列为

分别为外边界点坐标Y_j两个方向的坐标分量，1≤j≤N_o；令网格内部点的总数量为N_i，内部点的坐标序列为

分别为网格内部点坐标Z_k两个方向的坐标分量，1≤k≤N_i。

第二步，建立关于阻力变量d的翼型气动阻力函数

(f(d,S)为翼型S所受的总阻力，等于S上阻力变量d的积分，总阻力越小表明翼型风阻系数越小)；设置减阻收敛阈值ε(ε为浮点数且0＜ε＜1)和减阻过程的最大迭代步T_max；当迭代过程中翼型气动阻力函数值的变化比值(完全收敛时无变化，比值为1)大于减阻收敛阈值ε时，说明迭代收敛；最大迭代步T_max表示迭代超过T_max步后直接结束，用于控制计算总量和防止无法收敛的情况，为正整数，一般设置为100以内。

第三步，令当前迭代步T＝0，阻力函数旧值V0＝-1。

第四步，采用CFD求解软件计算流动控制方程R(W,S)＝0，得到流场变量W；根据翼型阻力公式g(W)计算阻力变量d＝g(W)，并将阻力变量d代入翼型气动阻力函数

得到阻力函数新值V，V＝f(d,S)。

第五步，若T≥1且V/V0＞ε，表明翼型的气动阻力函数值收敛到最小值，即翼型已达到最佳的风阻系数设计；或者T≥T_max，表明已达到最大迭代步数内翼型减阻的最好结果，这两种情况下转第十一步；否则转到第六步。

第六步，基于流场变量W，采用伴随方程求解软件计算连续伴随方程

(引入伴随变量ψ，使得气动阻力函数的导数只与翼型的几何变化相关)，得到伴随变量ψ。

第七步，将S上的N_w个内边界点坐标作为设计变量组。

第八步，采用SLSQP算法求解气动阻力函数f(d,S)关于设计变量组的最小值问题

得到网格内边界上N_w个内边界点最佳的位移量序列

其中ΔX_i为第i个内边界点最佳的位移量，

分别表示第i个内边界点在两个方向的位移分量，1≤i≤N_w。

第九步，根据翼型网格边界点的位移，采用增量式求解的RBF网格变形算法对翼型网格进行气动减阻变形，如图4所示：

9.1、由于远场外边界点始终静止不动，将N_o个外边界点的位移量(点的位移量是指该点相对于原来位置的偏移量)均设为(0,0)；令控制点坐标序列为空，令控制点位移量序列为空，令控制点数量N_c＝0；

9.2、从内边界点坐标序列

随机选择3个坐标，令这3个坐标对应的网格点为r₁、r₂、r₃，令第一控制点的坐标

令第二控制点的坐标

令第三控制点的坐标

将C₁、C₂、C₃加入控制点坐标序列，令第一控制点位移

令第二控制点位移

令第三控制点位移

将ΔC₁、ΔC₂、ΔC₃加入控制点位移量序列；

9.3、从外边界点坐标序列

随机选择3个坐标，令这3个坐标对应的网格点r₄、r₅、r₆，令第四控制点坐标

令第五控制点坐标

令第六控制点坐标

将C₄、C₅、C₆加入控制点序列，令第四控制点位移ΔC₄＝(0,0)，令第五控制点位移ΔC₅＝(0,0)，令第六控制点位移ΔC₆＝(0,0)，将ΔC₄、ΔC₅、ΔC₆加入控制点位移量序列；由此得到初始控制点坐标序列[C₁,C₂,C₃,C₄,C₅,C₆]以及对应的控制点位移量序列[ΔC₁,ΔC₂,ΔC₃,ΔC₄,ΔC₅,ΔC₆]，令N_c＝6，其中

分别表示第v个控制点在两个方向的位移分量，1≤v≤N_c；

9.4、采用公式(6)作为RBF网格变形算法的基函数：

其中ξ＝||x_p1-x_q1||/R(||x_p1-x_q1||表示坐标x_p1与坐标x_q1之间的欧式距离)，x_p1,x_q1泛指两个网格点坐标，支撑半径R表示网格点的影响范围(若两个网格点之间的欧氏距离大于R，则这两个点的移动互不影响)，支撑半径R设为翼型弦长的3-5倍；

9.5、根据控制点坐标序列中的边界点，由公式(7)计算得到维度为N_c×N_c的距离矩阵Φ_b,b中第p2行第q2列的元素Φ_b,b(p2,q2)，并令Φ_b,b(q2,p2)＝Φ_b,b(p2,q2)，其中1≤p2≤N_c,p2≤q2≤N_c；

Φ_b,b(p2,q2)＝φ(||C_p2-C_q2||) (7)

公式(7)中，||C_p2-C_q2||表示坐标C_p2与坐标C_q2之间的欧式距离；

9.6、令两个坐标维度方向上关于RBF函数权重系数的线性方程组的初始解分别为λ0¹＝[0,0,0,0,0,0]^T,λ0²＝[0,0,0,0,0,0]^T；

9.7、令坐标维度变量n＝1；

9.8、由控制点位移量序列中第n维度的位移分量构成维度为N_c×1的矩阵

其中1≤i1≤N_c，形成公式(8)的线性方程组：

9.9、以λ0ⁿ作为初始解，采用预条件投影的共轭梯度方法求解公式(8)的线性方程组，得到维度为N_c×1的权重系数矩阵λⁿ；

9.10、由权重系数矩阵λⁿ得到第n维度如式(9)的RBF函数sⁿ(X)：

φ(||X-C_j1||)是以||X-C_j1||为自变量的基函数，||X-C_j1||表示坐标X与坐标C_j1之间的欧式距离，

是权重系数矩阵λⁿ的第j1个元素；

9.11、令n＝n+1，若n≤2，转到第9.8步；否则转到第9.12步；

9.12、令内边界点变量i＝1，令当前最大位移误差σ_now＝-1；

9.13、将内边界点坐标X_i分别代入RBF函数s¹(X)和s²(X)，得到临时位移量(s¹(X_i),s²(X_i))；与X_i的实际位移量

进行比较，计算得到临时位移误差

若σ_no_w＜σ_tmp，令σ_now＝σ_tmp，令临时控制点坐标C_tmp＝X_i，令临时控制点位移ΔC_tmp＝ΔX_i，转9.14；若σ_now≥σ_tmp，直接转9.14；

9.14、令i＝i+1，若i≤N_w，转到第9.13步；否则转到第9.15步；

9.15、令外边界点变量j＝1；

9.16、将外边界点Y_j分别代入RBF函数s¹(X)和s²(X)，得到临时位移量(s¹(Y_j),s²(Y_j))；与点Y_j的实际位移量(0,0)进行比较，计算得到临时位移误差

若σ_no_w＜σ_tmp，令σ_no_w＝σ_tmp，令临时控制点C_tmp＝Y_j，令临时控制点位移ΔC_tmp＝(0,0)，转9.17；若σ_now≥σ_tmp，直接转9.17；

9.17、令j＝j+1，若j≤N_o，转到第9.16步；否则转到第9.18步；

9.18、设置位移误差阈值σ＝L/T_max，其中L是初始二维网格中边长的最小值，T_max是减阻过程的最大迭代步；

9.19、若当前最大位移误差σ_now＜位移误差阈值σ，表明由当前控制点序列计算得到的RBF函数满足精度要求，转到第9.23步计算内部点位移量；否则转到第9.20步更新控制点坐标序列；

9.20、将C_tmp加入控制点坐标序列，将ΔC_tmp加入控制点位移量序列，令N_c＝N_c+1；

9.21、在距离矩阵Φ_b,b的最后增加一行一列：根据控制点序列采用公式(7)计算得到距离矩阵Φ_b,b中第p3行第N_c列元素Φ_b,b(p3,N_c)的值，并令Φ_b,b(N_c,p3)＝Φ_b,b(p3,N_c)，其中1≤p3≤N_c；

9.22、令λ0¹＝[λ¹,0]^T,λ0²＝[λ²,0]^T分别作为两个坐标方向扩充后线性方程组的初解；转到第9.7步重新计算权重系数；

9.23、令内部节点变量k＝1；

9.24、将第k个内部节点的坐标Z_k分别代入RBF函数s¹(X)和s²(X)，得到该点位移量

并更新该点坐标Z_k，令

9.25、令k＝k+1，若k≤N_i，转到第9.24步；否则转到第9.26步；

9.26、令内边界点变量i＝1；

9.27、将第i个内边界点的坐标X_i分别代入RBF函数s¹(X)和s²(X)，得到该点位移量

更新第i个内边界点坐标X_i，令

9.28、令i＝i+1，若i≤N_w，转到第9.27步继续对翼型几何S上的内边界点进行更新；否则说明对翼型网格进行气动减阻变形处理完毕，得到气动减阻后的新翼型网格，转第十步；

第十步，令T＝T+1，V0＝V，转到第四步，对气动减阻后的新翼型网格继续进行翼型气动减阻。

第十一步，输出达到最佳的风阻系数或在最大迭代步数内翼型减阻达到最好结果的翼型网格和翼型几何(翼型网格是指图2(a)所示整体网格，翼型几何是指图2(b)梭形的网格内边界)，结束。

Claims (6)

1.一种基于改进径向基函数变形算法的翼型气动减阻方法，其特征在于包括以下步骤：

第一步，读入初始翼型对应的二维网格文件，二维网格的内边界对应初始翼型S，远场外边界距离内边界20倍翼型弦长以上；令二维网格边界点的总数量为N_b，其中内边界点总数量为N_w，内边界点的坐标序列为

X_i是一个内边界点的二维矢量坐标，

分别为X_i两个方向的坐标分量，1≤i≤N_w；令外边界点总数量为N_o，外边界点的坐标序列为

分别为外边界点坐标Y_j两个方向的坐标分量，1≤j≤N_o；令网格内部点的总数量为N_i，内部点的坐标序列为

分别为网格内部点坐标Z_k两个方向的坐标分量，1≤k≤N_i；

第二步，建立关于阻力变量d的翼型气动阻力函数

f(d,S)为翼型S所受的总阻力，等于S上阻力变量d的积分；设置减阻收敛阈值ε和减阻过程的最大迭代步T_max，ε为浮点数，T_max为正整数；

第三步，令当前迭代步T＝0，阻力函数旧值V0＝-1；

第四步，采用CFD求解软件计算流动控制方程R(W,S)＝0，得到流场变量W；根据翼型阻力公式g(W)计算阻力变量d＝g(W)，并将阻力变量d代入翼型气动阻力函数

得到阻力函数新值V，V＝f(d,S)；

第五步，若满足T≥1且V/V0>ε，或满足T≥T_max，转第十一步；否则转到第六步；

第六步，基于流场变量W，采用伴随方程求解软件计算连续伴随方程

得到伴随变量ψ；

第七步，将S上的N_w个内边界点坐标作为设计变量组；

第八步，采用SLSQP算法即序列最小二乘规划算法求解气动阻力函数f(d,S)关于设计变量组的最小值问题

得到网格内边界上N_w个内边界点最佳的位移量序列

其中ΔX_i为第i个内边界点最佳的位移量，

分别表示第i个内边界点在两个方向的位移分量，1≤i≤N_w；

第九步，根据翼型网格边界点的位移，采用增量式求解的径向基函数即RBF网格变形算法对翼型网格进行气动减阻变形：

9.1、将N_o个外边界点的位移量均设为(0,0)；令控制点坐标序列为空，令控制点位移量序列为空，令控制点数量N_c＝0；

9.2、从内边界点坐标序列

随机选择3个坐标，令这3个坐标对应的网格点为r₁、r₂、r₃，令第一控制点的坐标

令第二控制点的坐标

令第三控制点的坐标

将C₁、C₂、C₃加入控制点坐标序列，令第一控制点位移

令第二控制点位移

令第三控制点位移

将ΔC₁、ΔC₂、ΔC₃加入控制点位移量序列；

9.3、从外边界点坐标序列

随机选择3个坐标，令这3个坐标对应的网格点r₄、r₅、r₆，令第四控制点坐标

令第五控制点坐标

令第六控制点坐标

将C₄、C₅、C₆加入控制点序列，令第四控制点位移ΔC₄＝(0,0)，令第五控制点位移ΔC₅＝(0,0)，令第六控制点位移ΔC₆＝(0,0)，将ΔC₄、ΔC₅、ΔC₆加入控制点位移量序列；由此得到初始控制点坐标序列[C₁,C₂,C₃,C₄,C₅,C₆]以及对应的控制点位移量序列[ΔC₁,ΔC₂,ΔC₃,ΔC₄,ΔC₅,ΔC₆]，令N_c＝6，其中

分别表示第v个控制点在两个方向的位移分量，1≤v≤N_c；

9.4、采用公式(6)作为RBF网格变形算法的基函数：

其中ξ＝||x_p1-x_q1||/R，||x_p1-x_q1||表示坐标x_p1与坐标x_q1之间的欧式距离，x_p1,x_q1泛指两个网格点坐标，支撑半径R表示网格点的影响范围；

9.5、根据控制点坐标序列中的边界点，由公式(7)计算得到维度为N_c×N_c的距离矩阵Φ_b,b中第p2行第q2列的元素Φ_b,b(p2,q2)，并令Φ_b,b(q2,p2)＝Φ_b,b(p2,q2)，其中1≤p2≤N_c,p2≤q2≤N_c；

Φ_b,b(p2,q2)＝φ(||C_p2-C_q2||) (7)

其中||C_p2-C_q2||表示坐标C_p2与坐标C_q2之间的欧式距离；

9.6、令两个坐标维度方向上关于RBF函数权重系数的线性方程组的初始解分别为λ0¹＝[0,0,0,0,0,0]^T,λ0²＝[0,0,0,0,0,0]^T；

9.7、令坐标维度变量n＝1；

9.8、由控制点位移量序列中第n维度的位移分量构成维度为N_c×1的矩阵

其中1≤i1≤N_c，形成公式(8)的线性方程组：

9.9、以λ0ⁿ作为初始解，采用预条件投影的共轭梯度方法求解公式(8)的线性方程组，得到维度为N_c×1的权重系数矩阵λⁿ；

9.10、由权重系数矩阵λⁿ得到第n维度如式(9)的RBF函数sⁿ(X)：

φ(||X-C_j1||)是以||X-C_j1||为自变量的基函数，||X-C_j1||表示坐标X与坐标C_j1之间的欧式距离，

是权重系数矩阵λⁿ的第j1个元素；

9.11、令n＝n+1，若n≤2，转到第9.8步；否则转到第9.12步；

9.12、令内边界点变量i＝1，令当前最大位移误差σ_now＝-1；

9.13、将内边界点坐标X_i分别代入RBF函数s¹(X)和s²(X)，得到临时位移量(s¹(X_i),s²(X_i))；与X_i的实际位移量

进行比较，计算得到临时位移误差

若σ_now<σ_tmp，令σ_now＝σ_tmp，令临时控制点坐标C_tmp＝X_i，令临时控制点位移ΔC_tmp＝ΔX_i，转9.14；若σ_now≥σ_tmp，直接转9.14；

9.14、令i＝i+1，若i≤N_w，转9.13；否则转9.15；

9.15、令外边界点变量j＝1；

9.16、将外边界点Y_j分别代入RBF函数s¹(X)和s²(X)，得到临时位移量(s¹(Y_j),s²(Y_j))；与点Y_j的实际位移量(0,0)进行比较，计算得到临时位移误差

若σ_now<σ_tmp，令σ_now＝σ_tmp，令临时控制点C_tmp＝Y_j，令临时控制点位移ΔC_tmp＝(0,0)，转9.17；若σ_now≥σ_tmp，直接转9.17；

9.17、令j＝j+1，若j≤N_o，转9.16步；否则转9.18；

9.18、设置位移误差阈值σ＝L/T_max，其中L是初始二维网格中边长的最小值；

9.19、若当前最大位移误差σ_now<位移误差阈值σ，转9.23；否则转9.20；

9.20、将C_tmp加入控制点坐标序列，将ΔC_tmp加入控制点位移量序列，令N_c＝N_c+1；

9.21、在距离矩阵Φ_b,b的最后增加一行一列：根据控制点序列采用公式(7)计算得到距离矩阵Φ_b,b中第p3行第N_c列元素Φ_b,b(p3,N_c)的值，并令Φ_b,b(N_c,p3)＝Φ_b,b(p3,N_c)，其中1≤p3≤N_c；

9.22、令λ0¹＝[λ¹,0]^T,λ0²＝[λ²,0]^T分别作为两个坐标方向扩充后线性方程组的初解，转9.7；

9.23、令内部节点变量k＝1；

9.24、将第k个内部节点的坐标Z_k分别代入RBF函数s¹(X)和s²(X)，得到该点位移量

并更新该点坐标Z_k，令

9.25、令k＝k+1，若k≤N_i，转到第9.24步；否则转到第9.26步；

9.26、令内边界点变量i＝1；

9.27、将第i个内部节点的坐标X_i分别代入RBF函数s¹(X)和s²(X)，得到该点位移量

更新第i个内边界点坐标X_i，令

9.28、令i＝i+1，若i≤N_w，转到第9.27步；否则转第十步；

第十步，令T＝T+1，V0＝V，转第四步；

第十一步，输出达到最佳的风阻系数或在最大迭代步数内翼型减阻达到最好结果的翼型网格和翼型几何，结束。

2.如权利要求1所述的基于改进径向基函数变形算法的翼型气动减阻方法，其特征在于所述减阻收敛阈值ε满足0<ε<1，所述减阻过程的最大迭代步T_max设置为100以内。

3.如权利要求1所述的基于改进径向基函数变形算法的翼型气动减阻方法，其特征在于第四步所述CFD求解软件采用SU2 4.3.0以上版本中的SU2_CFD软件。

4.如权利要求1所述的基于改进径向基函数变形算法的翼型气动减阻方法，其特征在于第六步所述伴随方程求解软件采用SU2 4.3.0以上版本中的SU2_CFD及SU2_DOT软件。

5.如权利要求1所述的基于改进径向基函数变形算法的翼型气动减阻方法，其特征在于第八步所述SLSQP算法指SU2 4.3.0以上版本中的SLSQP优化器。

6.如权利要求1所述的基于改进径向基函数变形算法的翼型气动减阻方法，其特征在于第9.4步所述支撑半径R设为翼型弦长的3-5倍。

Images