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

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 1 、r 2 、r 3 ，令第一控制点的坐标 令第二控制点的坐标 令第三控制点的坐标 将C 1 、C 2 、C 3 加入控制点坐标序列，令第一控制点位移 令第二控制点位移 令第三控制点位移 将ΔC 1 、ΔC 2 、ΔC 3 加入控制点位移量序列；9.3、从外边界点坐标序列 随机选择3个坐标，令这3个坐标对应的网格点r 4 、r 5 、r 6 ，令第四控制点坐标 令第五控制点坐标 令第六控制点坐标 将C 4 、C 5 、C 6 加入控制点序列，令第四控制点位移ΔC 4 ＝(0,0)，令第五控制点位移ΔC 5 ＝(0,0)，令第六控制点位移ΔC 6 ＝(0,0)，将ΔC 4 、ΔC 5 、ΔC 6 加入控制点位移量序列；由此得到初始控制点坐标序列[C 1 ,C 2 ,C 3 ,C 4 ,C 5 ,C 6 ]以及对应的控制点位移量序列[ΔC 1 ,ΔC 2 ,ΔC 3 ,ΔC 4 ,ΔC 5 ,ΔC 6 ]，令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 1 ＝[0,0,0,0,0,0] T ,λ0 2 ＝[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 n 作为初始解，采用预条件投影的共轭梯度方法求解公式(8)的线性方程组，得到维度为N c ×1的权重系数矩阵λ n ；9.10、由权重系数矩阵λ n 得到第n维度如式(9)的RBF函数s n (X)：φ(||X-C j1 ||)是以||X-C j1 ||为自变量的基函数，||X-C j1 ||表示坐标X与坐标C j1 之间的欧式距离， 是权重系数矩阵λ n 的第j1个元素；9.11、令n＝n+1，若n≤2，转到第9.8步；否则转到第9.12步；9.12、令内边界点变量i＝1，令当前最大位移误差σ now ＝-1；9.13、将内边界点坐标X i 分别代入RBF函数s 1 (X)和s 2 (X)，得到临时位移量(s 1 (X i ),s 2 (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 1 (X)和s 2 (X)，得到临时位移量(s 1 (Y j ),s 2 (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 1 ＝[λ 1 ,0] T ,λ0 2 ＝[λ 2 ,0] T 分别作为两个坐标方向扩充后线性方程组的初解，转9.7；9.23、令内部节点变量k＝1；9.24、将第k个内部节点的坐标Z k 分别代入RBF函数s 1 (X)和s 2 (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 1 (X)和s 2 (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倍。