基于dffd网格变形技术的椭圆对称翼型多目标优化方法
技术领域
1.本发明属于气动外形优化领域，具体涉及了基于dffd网格变形技术的椭圆对称翼型多目标优化方法。
背景技术：2.在气动外形优化过程中的一个重要问题是如何针对新的气动外形快速有效地生成满足cfd计算需要的高质量计算网格。基于网格生成软件如icemcfd等录制脚本，再调用软件生成网格的方式是最直观的方式，但是耗费大量的网格生成时间，网格生成效率低下，直接后果是造成优化代价大。因此，有必要开发一套稳健快速的网格运动变形方法去克服优化过程中的这一难点。常用的网格运动方法有：基于delaunay背景网格的直接插值法、径向基函数插值法，无限插值法以及rbf/体样条插值和tfi的组合模式。这些方法对于气动优化中的小变形设计具有很好的鲁棒性，但当物体出现大幅度弯曲、扭转变形时，极容易出现物面网格边界线与物体交叉，造成网格负体积的出现。为了提高气动优化设计中的网格大幅度运动变形能力，本发明提出了一种基于直接控制自由造型变形(dffd)方法的网格运动变形技术。
3.自由造型变形(ffd)技术起源于计算机图形领域，它允许在二维或者三维空间内对物体进行变形操作，主要用于生成复杂的外形或使其变形。自由造型变形通过建立一个空间控制框架定义一个变形区域将物体包围，并使用经典bezier函数对物体表面进行参数化，直接移动控制框架的顶点对物体表面进行间接变形。由于物体变形是通过控制顶点来实现的，因此理论上它可以适用于对任意复杂的物体进行变形操作。在这种背景下，自由造型变形技术被广泛应用于工业优化设计领域。duvigneau使用自由造型变形设计了一个控制体来对气动计算网格进行变形，这种方法使用控制体的顶点作为设计变量来对机翼进行气动优化设计；menzel使用相似的方法对一个涡轮叶片进行设计；morris借鉴自由造型变形方法提出了域元法，将变形区域嵌入ffd控制框架中，通过移动控制框架定点，并对物面网格和空间网格进行rbf插值。上述气动优化设计方法中，通过移动控制点驱使气动分析网格变形，并且在优化过程中任意改变网格点，这种方法减少了设计变量、效率高、可以适应任意类型的计算网格。但是，任意改变控制点可能会造成气动外形不光滑，进而使气动性能降低；此外，任意移动的控制点不具有明晰的设计意义，以至于在优化之后很难反设计出真实外形。为了避免生成不光滑的气动外形，li将控制点组合起来，在这些组合点与几何尺寸之间建立运动联系，并将这种方法应用于高速火车的气动优化设计。gagnon等发展了一种双层自由造型变形法对机翼进行优化设计，并使用轴曲线对ffd控制体的控制点进行控制，使机翼的后掠、平面形状、扭转以及尖梢比及其计算网格发生变化。上述基于ffd方法的气动外形优化设计都考虑了几何尺寸和控制点运动之间的关系，为参数化网格运动提供了足够的支持，但是他们都没有形成气动设计参数和网格运动之间的联系。为了建立气动设计参数和网格运动之间的联系，li提出了一种新型的ffd基的参数化网格运动映射方法，这种方法减少了设计变量的数量，增强了设计效果。
4.受以上研究成果的启发，同时为了改进上述方法中的不利因素，本专利提出基于直接控制自由造型变形(dffd)的网格运动方法。
技术实现要素：5.目前技术是基于网格生成软件如icemcfd等录制脚本，再调用软件生成网格的方式是实现网格运动最直观的方法，但是耗费大量的网格生成时间，网格生成效率低下，直接后果是优化代价大、时间长。现在常用的网格运动方法能很好地适应小变形气动优化设计，但对于在大幅度弯曲、扭转变形极容易出现物面网格边界线与物体交叉，造成网格负体积的出现。基于直接控制自由造型变形(dffd)的方法可以实现几何形状的大幅度变形，将其应用于网格运动变形具有无可比拟的优势。
6.标准ffd算法具有很多优良特性，但是也存在着一些局限性：ffd方法只是一种变形工具，并不能直接控制物体本身，也就是设计变量与物体外形之间没有直接的映射关系，设计变量的物理意义不够明晰，不能准确地生成设计人员想要的几何形状；要想进行复杂外形的变形，需要大量的控制点，造成设计变量太多，设计效率低下；bernstein多项式的局部变形能力不强，对局部变形能力捕捉不足等。
7.本发明提出的基于dffd方法的网格运动算法能够直接建立设计变量和几何外形以及所带来网格运动之间的联系，以极少的设计变量就能实现外形参数化和空间网格运动变形，通过调整dffd算法内b样条函数的次数可以实现大幅度几何外形变形和空间cfd计算网格的运动变形。主要用于气动优化中使用最小计算代价生成高质量的cfd计算网格。
8.本发明的技术方案：
9.基于dffd网格变形技术的椭圆对称翼型多目标优化方法，步骤如下：
10.1)选定网格变形区域，并布置ffd控制点，如图2所示，其中实心点为ffd控制点，四个相邻的控制点与连接它们的粗实线组成一个控制体，所有控制体组合起来即为网格变形区域。网格变形区域可以是网格的局部区域，也可以包含整个网格；
11.2)在外形上选定若干直接操作点，设置b样条函数的阶次，初始化直接操作点和网格点的局部坐标，对于不规则布置的ffd控制框架，需要将它们转换成一个优化问题来求解近似的局部坐标s，t和u值；
12.3)使用直接操作点局部坐标计算b样条系数矩阵；
13.4)用系数矩阵和直接操作点位移反计算ffd控制点的位移；
14.5)利用移动后的控制点计算变形区域空间网格变形后的位置。
15.dffd方法与基本型ffd的流程图比较见图1。
16.当设计人员给目标物体上的直接控制点施加一定的位移δx时，会在控制框架上产生相应的位置变化。设有q个直接控制点，其初始位置为s
f
，f＝1,2,
…
,q，与之相对应的局部坐标为(s
f
,t
f
,u
f
),需要将其移动到目标位置t
f
，f＝1,2,
…
,q，那么有以下约束成立：
[0017][0018]
其中，p
ijk
为控制框架中控制顶点的初始坐标，r
ijk
(s
f
,t
f
,u
f
)＝b
il
(s
f
)b
jm
(t
f
)b
kn
(u
f
)，δ
ijk
是控制框架中每个控制顶点的位移量。b
il
(s
f
)为l次bernstein多项式基函数，b
jm
(t
f
)为m次bernstein多项式基函数，b
kn
(u
f
)为n次bernstein多项式基函数。l、m和n分别为三个坐标方向上的控制点数量减1。
[0019]
此问题的目标是求解δ
ijk
，满足公式(1)的δ
ijk
可能有很多种，需要在最小二乘条件下使各个控制顶点的位移平方和达到极小值，因此定义拉格朗日函数：
[0020][0021]
其中λ
f
＝(λ
f1
，λ
f2
，λ
f3
)
t
，f＝1,2,
…
,q是拉格朗日乘子。因此，有
[0022][0023]
将上式写成矩阵形式：
[0024][0025]
其中d＝[t1‑
s1,
…
,t
q
‑
s
q
]
t
，δ＝[δ
000
,δ
001
,
…
,δ
lmn
]
t
，t
f
＝[s
f
,t
f
,u
f
]，λ
f
＝[λ
f1
,λ
f2
,λ
f3
]。
[0026]
将公式(4)中的下式代入上式，可以得到d＝r
t
rλ，因此λ＝(r
t
r)
‑1d。最终可以得到控制顶点的位移为：
[0027]
δ＝r(r
t
r)
‑1d＝r(r
t
r)
‑1(t
‑
s)
ꢀꢀꢀ
(5)
[0028]
求得控制框架各控制点的位移之后，将其代入下式，即可以计算出目标物体上及空间网格所有点的新坐标。
[0029]
x
′
＝x+δx＝b(p+δ)
ꢀꢀꢀ
(6)
[0030]
其中，x为目标物体上的初始坐标，δx为目标物体位移量；向量矩阵b为b
q
×
(l+1)(m+1)(n+1)
，是b样条基函数的矩阵向量，p为控制顶点初始坐标。
[0031]
本发明的关键点是将直接控制自由造型变形(dffd)方法应用在气动优化设计中的网格运动变形上。基于dffd方法的新型网格运动技术的核心思想是：选择物体表面的若干点作为直接控制点，移动至所要求的位置，反算出计算网格控制框架顶点的位置变化，计算出控制框架内其它空间网格点的新坐标。本发明实现了外形参数化和网格运动一体化实施，并通过调整dffd中的b样条次数实现网格变形能力和变形区域的控制。
[0032]
本发明的有益效果：此技术解决了翼型气动优化过程中传统网格运动方法不能适应外形剧烈变形的问题，大幅度提高了网格的运动变形能力。基于直接控制的自由造型变形方法的最大的优点是用户只需要在变形目标周围定义一个控制框架就能实现复杂外形的变形，并且可以在物体表面特定区域实现精确变形，方便用户对重点区域的关注；此外，
网格的变形过程与网格的空间分布特性并不相关，这种特性对于具有多种结构类型的网格实现运动变形非常有利，因此可以对具有任意空间结构的网格实现网格运动；然后，对于网格发生重叠的区域，可以通过加入一个直接控制点实现网格质量的改善；最后，基于直接控制点的网格变形技术可以方便加入几何约束，具有更强的适应性和实用性。通过网格变形可以看出，基于dffd方法的网格运动技术具有强大的变形能力，对于复杂构型或者网格运动，具有很强的稳健性和更高的变形效率。
附图说明
[0033]
图1直接控制ffd与基本型ffd的流程图比较；
[0034]
图2网格变形区域及ffd控制框架；
[0035]
图3两点移动示意图；
[0036]
图4两点运动外形及网格变形，其中，(a)为三次b样条，(b)为五次b样条，(c)为七次b样条；
[0037]
图5三次b样条时七点运动外形及网格变形；
[0038]
图6五次b样条时七点运动外形及网格变形；
[0039]
图7七次b样条时七点运动外形及网格变形；
[0040]
图8优化设计空间及直接控制点；
[0041]
图9椭圆对称翼型的网格变形区域及ffd控制框架；
[0042]
图10两种极端变形情况下的参数化翼型和空间网格运动，其中，(a)为情况一，(b)为情况二；
[0043]
图11椭圆翼型优化后的pareto解；
[0044]
图12优化翼型与基准翼型外形对比；
[0045]
图13设计点一翼型压力分布比较；
[0046]
图14设计点二翼型压力分布比较；
[0047]
图15优化翼型和基准翼型的cl,cd和cm对比(ma＝0.6,re＝4.8e6)，其中，(a)为cl，(b)为cd，(c)为cm；
[0048]
图16优化翼型和基准翼型的cl,cd和cm对比(ma＝0.5,re＝4.0e6)，其中，(a)为cl，(b)为cd，(c)为cm；
[0049]
图17基准翼型和优化翼型极曲线对比一；
[0050]
图18基准翼型和优化翼型极曲线对比二；
[0051]
图19椭圆翼型四目标优化后的pareto解；
[0052]
图20paretoa翼型与基准翼型外形对比；
[0053]
图21椭圆翼型五目标优化后的pareto解；
[0054]
图22paretoa翼型与基准翼型外形对比。
具体实施方式
[0055]
以下根据实施例和附图对本发明的技术方案进行进一步说明。
[0056]
使用一个二维翼型的网格运动来详细说明dffd方法的网格运动变形过程和能力。其中，采用b样条函数来获取系数矩阵。图2显示了一个绕单位弦长的椭圆翼型生成的四边
形网格，将需要变形的网格区域和翼型置于一个简单的网格控制框架内，控制框架的控制点沿坐标轴均匀分布，见图2所示正方形点。这里重点介绍网格变形运动的实现过程以及b样条次数对网格运动变形能力的影响。
[0057]
在优化设计过程中，一般会直接控制两个及以上的点进行变形设计，这样才能为外形设计提供足够的设计变量。为了研究基于dffd方法的网格运动变形技术的效果，在图3中定义了四个直接控制点，其中椭圆翼型前缘和后缘直接控制点保持固定，翼型中间两个控制点各自向上移动0.1倍单外弦长。图4分别给出了使用3次、5次和7次b样条基函数时的网格变形情况。可以看出，3次b样条产生的翼型表面有凹凸不平的形状；样条次数增大为5时，翼型表面的凹凸变小；样条次数继续增大为7，翼型上下表面都只有向上的凸起，变得更光滑。与之相对应的网格变形区域和控制框架变化区域变大，发生位移的控制顶点数量增多，单一控制顶点的位移幅度减小。利用网格生成软件icemcfd的网格质量检查功能对运动网格的质量进行了检查，网格质量对比数据列于表1中。可以看出网格运动之后的网格质量都有一定的降低，但是最小角度都大于50度，最小quality指标也高于0.8，而且随着b样条次数的增高，网格质量指标也随之改善。
[0058]
表1网格质量检查
[0059][0060]
为了更好地产生网格大变形运动，在翼型的上下翼面布置七个可以移动的直接控制点，以及前、后缘各一个位置约束点，产生如图5所示的翼型外形变化。可以看出，在弦向0.6位置处，翼型外形发生了0.4倍弦长的剧烈变化，随之空间网格也发生了剧烈变化，在使用三次b样条基函数时，空间网格出现了大面积的重叠区域；随着样条增加为5次时(图6所示)，网格重叠区域急剧变小；当样条次数为7时，重叠区域完全消失，如图7所示。
[0061]
实施例：
[0062]
椭圆翼型的气动设计需要考虑许多的设计因素和设计指标，传统的设计方法很难在多个方面取得气动性能的提升，有必要开展椭圆翼型气动多目标优化设计。选择一款16％厚度的椭圆翼型作为基准翼型进行多目标气动优化设计。
[0063]
在椭圆翼型的上、下表面各布置11个可变直接控制点，前缘和后缘各布置一个固定的直接控制点，一共24个直接控制点，直接控制点具体位置见图8中的实心圆点。鉴于椭圆翼型需要保持前后对称，选择上、下翼面的前六个可移动直接控制点作为设计变量，总设计变量数为12，每个设计变量的变化区间为其z方向坐标的
±
30％，具体的设计空间范围由图8中的虚线表示的区间确定。
[0064]
图9中给出了dffd的控制框架，其中实心较大的顶点为控制顶点，数量为45x31。由于本发明优化设计中采用的直接控制点较多，对翼型附近的控制顶点进行了加密处理，防止出现龙格现象；为了在直接控制点移动后，翼型的前后对称变形，所有控制顶点坐标关于(0.5,0)位置左右对称。当图8中所示的直接控制点发生移动后，控制框架顶点产生相应位
移，使控制框架内的非直接控制点产生移动。此过程同时实现了翼型的参数化和空间网格运动，简化了优化步骤。
[0065]
采用6次b样条dffd方法对翼型进行参数化和网格运动，对于优化过程中可能出现的两种极端变形情况的翼型参数化和空间网格运动进行了验证。最终翼型参数化和网格运动结果见图10，从图中可以看出翼型外形没有发生龙格现象，而且翼型附近和空间的网格分布良好，这说明基于dffd方法的网格运动变形技术适用于本优化算例。
[0066]
综合考虑椭圆对称翼型的主要工作状态，归纳出如表2列出的2个设计点和2个设计目标，其中最小化目标一是旨在降低悬停时翼型的阻力系数；最小化目标二是为了减小固定翼巡航状态下的阻力系数。建立优化数学模型为：
[0067]
min f1(x)＝cd1[0068]
min f2(x)＝cd2[0069]
s.t.x∈(x
l
,x
u
)
[0070]
g
m
(x)≥0，m＝1,2,
…
,n
m
[0071]
其中下标1和2对应表2中的设计点序号，表2中的“*”表示基准椭圆翼型的气动力值，对翼型进行最大厚度约束。
[0072]
采用lmpso算法作为优化算法，lmpso的参数保持与基准函数测试时相同，种群规模为60，循环次数为50，调用目标函数次数为3000次。多目标优化之后获得的pareto前沿如图11所示，从图中可以看出所有的pareto解在两个设计目标上相较于原始翼型都有很大改进，这说明对椭圆对称翼型的多目标优化设计是有效的。图中pareto解分布表明目标二上的优化改善量大于优化目标一，因此选择左数第三个pareto解作为最优解进行比较分析，下文中此翼型皆称为优化翼型。
[0073]
表2 16％厚度椭圆翼型的两个优化设计目标
[0074][0075]
图12给出了基准椭圆翼型和选择的优化翼型的外形对比。可以看出优化翼型前后缘都变薄，上翼面和下翼面形状不再对称，上翼面中间凸起，最大厚度变大；下翼面最大厚度略有减小，还出现了两个对称的加载区；相对于基准翼型的中心线为一条水平直线，优化翼型中心线向上凸起。表3给出了两个设计点的优化前后气动力对比，与基准翼型相比，优化翼型的最大厚度略有增加，气动力特性都有改善。设计点一的优化目标阻力系数减小了11.60％(
‑
20.6counts)，减阻效果非常明显；同等升力下迎角减小了0.695度，这也是阻力减小的主要原因；力矩系数的绝对值减小了84.86％。优化翼型设计点二的阻力系数降低了14.21％(
‑
21.4counts)；同等升力下迎角降低了0.767度；力矩系数的绝对值降低了66.53％。可以看出，优化翼型在两个设计点时的阻力减小的同时，俯仰力矩系数绝对值的减小量都很大，这对飞行控制非常有利。
[0076]
图13和图14为设计点一和设计点二时优化前后的翼型表面压力系数分布对比，图
中显示优化翼型上翼面前缘和后缘附近吸力减小，中部吸力增大；下翼面两个气动力加载区域的效果明显，后缘鱼尾状交叉区域面积变小，翼型后半段的升力增加，这是优化翼型两个设计点力矩系数减小的主要原因。图15和图16给出了两个设计点的升力系数、阻力系数和力矩系数随迎角的变化曲线，可以看出设计点一时优化翼型的升力特性明显优于基准翼型，最大升力系数增大，但是失速迎角略有减小；优化翼型的力矩在小迎角时变为低头力矩，在更大的迎角范围内具有更小的绝对值；优化翼型的阻力系数在小迎角时优于基准翼型，但是大迎角时劣于基准翼型。设计点二时优化翼型的升力系数曲线线性段性能明显优于基准翼型，但是最大升力系数和失速迎角都劣于基准翼型；阻力系数和力矩系数的表现和设计点一类似。在大迎角时优化翼型升、阻力系数变差主要是因为翼型前缘半径变小引起的，这也说明椭圆翼型在大迎角时升、阻力特性可能与小迎角时的力学特性冲突，应将其作为优化目标进行优化设计。
[0077]
表3优化翼型与基准翼型设计点一气动性能对比(cl＝0.6)
[0078][0079]
图17和图18分别给出了悬停状态和直翼前飞状态下优化前后翼型的升力、阻力极曲线。从图中可以看出，在升力系数为0.6时和升力系数为0.4时优化翼型阻力系数都有明显降低。
[0080]
1)椭圆翼型四目标优化设计结果
[0081]
表4给出了四目标优化的设计点和优化目标，以此为基础建立数学模型：
[0082]
min f1(x)＝cd1[0083]
min f2(x)＝cd2[0084]
min f3(x)＝cd3[0085]
min f4(x)＝cd4[0086]
s.t.x∈(x
l
,x
u
)
[0087]
g
m
(x)≥0，m＝1,2,
…
,n
m
[0088]
表4 16％厚度椭圆翼型的四个优化设计目标
[0089][0090]
表5优化翼型与基准翼型四个设计点的气动性能对比
[0091][0092][0093]
2)椭圆翼型五目标优化设计结果
[0094]
表6给出了五目标优化的七个设计点和优化目标，以此为基础建立数学模型：
[0095]
min f1(x)＝cd1[0096]
min f2(x)＝cd2[0097]
min f3(x)＝cd3[0098]
min f4(x)＝cd4[0099]
min f5(x)＝[(cd5‑
cd6)2+(cd6‑
cd7)2]/(δma)2[0100]
s.t.x∈(x
l
,x
u
)，g
m
(x)≥0，m＝1,2,
…
,n
m
[0101]
表6 16％厚度椭圆翼型的五个优化设计目标
[0102][0103]
表7优化翼型与基准翼型七个设计点的气动性能对比
[0104]
[0105]