2. 中国航空研究院, 北京 100012
近些年来，随着大型运输机、大展弦比长航时无人机的快速发展，弹性变形对飞机气动特性的影响不容忽视。静气动弹性变形十分重要是因为它们主宰着定常飞行状态下的载荷、升力分布、阻力、操纵面效率、飞机配平特性以及静态操稳特性。起飞和着陆构型的气动弹性变形对全机的起降性能有着直接影响，进而影响飞机的安全性、经济性和运营成本。波音公司的Garner等人认为[1]，对一架典型的双发喷气式运输机，在同一迎角下升力系数增加0.1相当于进场迎角减小1°，即可缩短起落架而使飞机重量减小635 kg；在同一进场速度下，最大升力系数增加1.5%相当于载重量可增加3 000 kg；起飞状态的升阻比提高1%相当于载重可增加1 270 kg或航程增加约280 km。弹性变形可使增升装置的缝道参数发生明显变化[2]，而这方面的相关研究还较少，现阶段静气弹的研究大多是针对飞机巡航状态，但是起飞和着陆阶段同样决定着飞机设计的成功与否。对于增升装置等复杂构型的气动弹性问题在具体型号研制过程中已经引起研制单位的关注。因此寻求一种可以应用在详细设计阶段或是适航审定阶段的相对精确的静气动弹性分析方法十分必要。
准确预测静气动弹性结构响应特性是飞机设计中的一个技术难点。它关系到准确确定飞行器的操纵性、安定性的空气动力学数据、飞行器的空气动力载荷数据、保证飞行范围内不出现发散和确定机翼的型架外形等一系列复杂的设计问题[3, 4]。对于起飞、着陆等增升构型由于其流动特性非常复杂，而且几何细节流动对增升装置流动具有重要影响[5, 6]。所以传统的面元法、Euler方程等方法已不能满足精细设计阶段的计算精度需求，而准确的流体域分析是开展精确气动弹性分析的基础。随着计算机硬件的发展，多CPU工作站和大规模并行机群的研制和应用，以及并行计算的成熟发展，使得飞行器气弹的求解转向并行计算方向，使得计算时间大大缩短，效率大大提高，也使得基于Navier-Stokes方程求解器的气动弹性分析成为可能。
本文针对大型客机增升装置构型，发展了一种可以应用于飞行器详细设计阶段、适航审定阶段的相对精确的静气动弹性分析方法。通过耦合通过风洞试验验证的N-S方程求解器和结构静力学方程，
应用RBF数据交换技术和基于RBF插值的动网格方法，引入加速收敛因子，建立静气弹分析方法。除此之外，为了解决复杂构型的数据交换问题，建立了一种混合局部/全局的数据交换方法。最后，以第二届阻力预测会议的DLR-F6标模和带前缘缝翼、后缘襟翼的某大型客机着陆构型为例进行静气动弹性分析，验证文中静气动弹性分析方法的可行性和对模型的鲁棒性。
1 理论基础以下针对于静气动弹性分析涉及的几方面问题分别进行研究，建立分析增升装置复杂构型弹性变形的工作基础。
1.1 CFD理论基础精确的计算流体力学数值模拟技术是进行准确静气动弹性分析的基础。增升装置模型复杂，通常包含很多几何细节，例如包含襟翼、襟翼滑轨、缝翼、缝翼滑轨等。为了提高增升装置构型的静气动弹性数值模拟精度，对其流场进行准确的数值模拟十分必要[7]。本文拟采用能RANS方程进行流场数值模拟，应用有限体积法进行空间离散，隐式时间推进求解，采用k-ω SST 湍流模型[8]。为了提高计算效率，采用多重网格加速收敛技术和并行计算的流场求解器。
积分形式的N-S方程[9]形式：
图 1为某大型客机着陆构型CFD计算与风洞试验结果对比曲线。
 |
| 图 1 升力系数对比 |
图中CFD计算得到的最大升力系数和升力线斜率均与试验结果吻合较好,可得以此CFD求解器为基础建立的静气弹方法是可信赖的。
1.2 结构有限元建模与分析建立正确且合理的有限元模型,是取得可信计算分析结果的前提。在建模工作中,首先应遵循力学等效原则,同时兼顾计算精度、计算效率和经济性。对于复杂构型来讲,直接采用梁元素模拟结构产生一个问题,这样计算得到的刚度分布较为粗糙。因此,采用更加精细的结构有限元建模方法十分必要。对于复杂结构模型,往往需要采用多种有限元素组合的方式描述。当计算模型中选用多种元素组合的形式建模时,所有元素在节点处必须具有相同类型和相同数量的自由度,另外还应注意相邻边界上元素形状函数的协调,否则会影响计算的收敛性。
1.3 流固耦合数据交换方法气动弹性分析过程是不断进行气动和结构耦合的过程,为了保证气动弹性分析的稳定性,一般采用气动、结构各自较成熟的分析方法,针对于不同建模原则的气动、结构模型,在耦合界面处需要进行数据交换[9]。研究表明,在各种散乱数据插值方法中,径向基函数(RBF)的结果最能使人满意[10, 11]。最早提出了将径向基函数插值应用到飞行器跨音速飞行时的静气动弹性变形中[11],提出用于满足平动和转动的线性插值多项式。RBF插值方法[10, 11]可以用于离散点之间的数据转换,具有实现简单、结果准确和耗时少等优点。
三维RBF方法插值公式为:
式中:基函数选取自由,根据对插值基函数的研究[12],文中选取基函数CP C0:φ(‖x‖)=(1-‖x‖)2,其中‖x‖=√x2+y2+z2。
以位移插值为例,只要确定(2)式中的N+4个未知数a0,a1,a2,a3,F1,F2,…,FN,可以确定数据转换矩阵。
数据交换方式主要有2类：①局部插值,这种插值方式有时会造成某些局部失真;②全局插值,此方法的特点为距离近的点分的多一些,距离小的点影响少一些,但是针对于复杂构型来讲,对建模要求较高。
本文通过对2种不同数据交换方式的研究,结合2种插值方式的优点,建立一种混合局部／全局的插值方式。具体来讲,对于全体模型采用分区域的数据交换方式,即将复杂构型分成不同部件,而在同一部件上采用自身全局插值。这样既可以降低对有限元模型的要求,也可以保证在局部作用载荷传递不失真。
以下通过采用某增升装置构型对全局数据交换方式和混合局部／全局数据插值方式进行研究。气动表面网格与结构有限元模型分别如图 2、图 3所示,采用不同的数据交换方式进行气动载荷传递,进而将结构弹性变形传递给气动表面网格,插值结果分别如图 4所示。
 |
| 图 2 气动表面网格 |
 |
| 图 3 结构有限元模型 |
 |
| 图 4 2种数据交换方式对比 |
由图 4 2种不同数据交换方式得到的气动表面网格位移云图可知,针对相同的结构有限元模型,采用不同的数据交换方式将气动载荷传递给结构模型,得到不同的变形结果。从云图上可以得出,尤其对于相邻较近的2个部件会出现明显的不同位移。在图 4左部放大部分可以得出,这部分弹性变形不光滑。而由图 4右部可以得出,采用文中改进的混合局部／全局数据交换方式得到光滑的气动表面网格。
1.4 动网格技术为了节省新网格生成的时间和解决大变形问题,需要发展一种高效的动网格方法,这种方法必须能很好的适应大展弦比弹性机翼产生的大变形情况。通过对多种动网格方法的研究,应用RBF插值技术开展一种兼顾变形能力、变形效率的动网格方法。RBF插值首先被应用于流固耦合分析中的界面参数传递[13],因此,将这种传递扩展到全体流体域内的就可以实现气动网格的更新。应用径向基函数实现网格更新时,仅需要各网格节点的坐标,不需要气动网格连接的信息,这样数据结构简单,操作方便,并且适用于分布式存储和并行计算等。以下将对本文采用的RBF动网格方法进行详细介绍,并通过算例验证此动网格方法的变形能力。
本文应用在网格更新方面的RBF插值基函数中范数的定义与进行流固耦合界面进行数据交换时有所不同,选用基函数CP C0、范数具体表达式如下[13]:
根据紧支函数的特点,函数值随着离中心距离的增加而减小,该距离达到某一特定值r(称为紧支半径support radius)后函数值恒为零,即当‖x‖>1时,函数值为0。紧支函数的特点与动网格计算的要求一致,在计算过程中的求解矩阵呈稀疏、带状分布,有利于大型问题的求解。本文选取的基函数,方程组的系数矩阵对称正定,应用共轭梯度法实现快速求解[14]。
值得注意的是,网格变形和CFD／CSD数据交换需要满足的要求是不同的,CFD／CSD数据变换矩阵需要保证准确的力和力矩平衡准则而包含线性项,而网格变换矩阵为了避免边界网格的移动和非物理力的传递而不需要线性项。所以要分别建立CFD／CSD数据变换矩阵和动网格变换矩阵[13, 15]:
式中：L为动网格变换矩阵(X代表网格点坐标,下标mesh代表网格点,下标a代表表面网格点)。
为了验证动网格方法的变形能力和对增升装置等复杂构型的适应性,选取增升构型网格进行动网格程序测试,图 5为网格变形前后的网格拓扑对比,图中,粗实线代表变形前网格拓扑,细实线代表变形后网格拓扑。由图可得,网格变形前后网格拓扑结构不变,并且本文动网格程序可以在各部件间缝道的变形情况下依然得到质量较好的计算网格。
 |
| 图 5 某增升装置变形前后网格拓扑 |
本文应用RBF混合全局／局部数据交换技术进行气动载荷和结构弹性变形之间的数据交换,并应用基于RBF插值方法的动网格方法,开展了一种适合于起飞、着陆等复杂构型的静气动弹性变形分析方法。
2.1 静气动弹性分析流程本文建立的静气动弹性分析方法思路为:计算定常流场直至收敛;将收敛的气动载荷传递给结构有限元模型;利用结构静力学求解程序求解结构变形;将结构变形传递给气动表面网格节点;更新气动网格,得到新的气动分析模型,并计算流场,如此循环直至满足预设的收敛条件,得到最终结构变形、变形后新的气动载荷分布。
分析流程如下图所示:
 |
| 图 6 静气弹分析流程图 |
图 6中,CFD solver为经过可靠性验证的流体域N-S方程求解器;CSD solver为结构静力学求解器;Program 1为气动载荷传递程序,根据虚功原理利用RBF插值方法完成;Program 2为结构弹性位移插值程序,利用RBF插值方法完成;Mesh Motion代表动网格程序,采用1.4节中介绍的动网格方法。
2.2 加速收敛技术进行静气动弹性分析是不断进行气动、结构耦合直至收敛的过程,为了加速收敛速度,文中引入加速收敛因子c,(6)式中Δ代表每次迭代弹性变形量,下标i代表迭代次数,c为加速收敛因子。
经过研究本文采用不同的迭代之间采用不同的加速收敛因子c,i=1时c=0.75,i>1时c=1.0。
3 算例分析为了验证本文建立的静气动弹性分析方法的可行性和准确性,以下选用DLR-F6翼身组合体为例进行静气动弹性分析,并与第二次阻力预测会议的结果进行对比。进一步选用某经过刚性风洞试验的着陆构型验证本文方法的可行性和鲁棒性,得到了一些很有意义的结论。
3.1 算例一:程序验证选取DLR-F6翼身组合体模型,对其进行静气动弹性分析。机翼采用参考文献[16]中提到的结构参数进行结构有限元建模。材料的弹性模量E=2.05×1011,泊松比ν=0.3。此算例中气动计算时考虑了机身对气动力的影响,在弹性结构有限元分析时认为机身为刚性的,只考虑机翼的弹性变形。气动网格、结构有限元模型分别如图 7、图 8所示。
 |
| 图 7 气动表面网格及空间拓扑 |
 |
| 图 8 结构有限元模型 |
本算例的计算状态为:Ma=0.75,Re=3.0×106,Cl=0.5。选取收敛条件为:节点位移收敛条件ε1=Δi+1-Δi≤10-5,升力系数收敛条件ε2=Cli+1-Cli≤10-4,i代表迭代次数,i=1,2,3…。
图 9、图 10分别是后缘挠度与展向扭转角和第二次阻力预测会议的参考文献[16]的对比图。由图可知,弹性变形引起的负扭转角随着展向的增加而增大,在翼梢处达到最大值。通过对比可得后缘挠度和展向扭转角分布捕捉都比较准确,5个站位的后缘挠度和展向扭转角都与参考值吻合得比较好,95％展向位置的挠度与参考值相差6.6％,95％的扭转角与参考值相差3.7％,表明本文的计算方法是可信的。
 |
| 图 9 后缘挠度对比图 |
 |
| 图 10 展向剖面扭转角对比 |
为了研究弹性变形对大展弦比机翼增升装置的气动效率影响,本算例对某三段着陆构型进行了静气动弹性数值模拟,分析弹性变形对着陆构型的影响。并根据模拟结果分析了弹性变形对增升效率的影响。由于着陆型较复杂,机翼部分由前缘缝翼、主翼、后缘襟翼以及缝翼滑轨、襟翼滑轨等组成,因此,本算例的弹性变形分析过程中采用文中发展的混合全局／局部的数据交换方式。为了提高静气弹收敛效率,引入加速收敛技术。
气动网格与结构有限元模型分别见1.3节中的图 2和图 3。本算例的计算状态为:Ma=0.2,Re=1.982×107,α=8.0°。收敛条件为:节点位移收敛条件ε1=Δi+1-Δ≤10-5,升力系数收敛条件ε2=Cli+1-Cli≤3×10-4。
图 11为Z方向弹性变形量云图。
 |
| 图 11 弹性变形量云图 |
可以得出，机翼的弹性变形主要表现为弯曲变形和负扭转变形，这是典型的后掠机翼的静气弹现象。从图中可以看出,机翼内翼段变形量很小，从kink处沿展向变形量逐渐增大，并在翼尖达到最大值0.62 m。由于弹性变形的影响，机翼沿流向各个剖面产生负的扭转角。
图 12为机翼展向剖面扭转角曲线图。
 |
| 图 12 正则化展向剖面弹性扭转角及变形量 |
可以得到，随着展向的增加，机翼弹性变形产生的剖面负扭转角越大。而且，随着展向的增加，各剖面的弹性扭转角呈现非线性变化。
图 13为弹性变形前后外侧后缘襟翼缝道搭接量对比曲线。
 |
| 图 13 弹性变形前后外侧襟翼缝道搭接量对比 |
可以得到，后缘外侧襟翼的缝道搭接量与初始构型比较产生明显的变化。对于着陆构型，缝道参数的变化将直接影响着着陆构型的气动效率，进而影响其运行的安全性和经济性。
本算例验证，弹性变形使此着陆构型缝道参数变化，而这些变化直接影响飞行工况下着陆特性。在本文的计算状态下，考虑弹性变形后的升力系数比刚性构型较小了1.0%。因此，在飞机详细设计阶段开展考虑弹性变形影响的增升装置设计十分必要。
4 结 论通过本文的研究得到以下结论：
1) 对于由较少部件组成的简单构型进行静气动弹性分析时，可以采用全局数据交换方式；而对于复杂构型，全局数据交换方式传递的数据会出现明显的误差。经过分析，本文建立的全局/局部混合方式的数据交换方式可以很好地实现复杂构型静气动分析过程中的数据传递。
2) 基于RBF插值技术发展了动网格技术，并引入并行技术，提高动网格的效率和变形能力。
3) 静气弹分析过程中根据计算收敛特点引入了加速收敛因子。
4) 以某一大型客机着陆构型为例，研究了弹性变形对增升装置增升效率的研究。在本文的计算状态下，弹性变形使着陆构型的升力系数减小了约1.0%。
| [1] | Garner P L, Meredith P T, Stoner R. Areas for Future CFD Development as Illustrated by Transport Aircraft Applications[R]. AIAA-1991-1527 |
| [2] | Vanderburg J W, Eliasson P, Delille T, et al. Geometric Installation and Deformation Effects in High-Lift Flows[J]. AIAA Journal, 2009, 47(1): 60-70 |
| Click to display the text | |
| [3] | Liu Y, Bai J Q, Hua J. Static Aeroelasticity Analysis of High Aspect Ratio Wing and the Jig-Shape Design Based on Multi-Block Structural Grid[J]. Advanced Engineering and Materials, 2013, 302: 377-383 |
| Click to display the text | |
| [4] | 万志强, 邓立东, 杨超, 等. 基于非线性试验气动力的飞机静气动弹性响应分析[J]. 航空学报, 2005, 26(4): 439-455 Wan Z Q, Deng L D, Yang C, et al. Aircraft Static Aeroelstic Response Analysis Based on Nonlinear Experimental Aerodynamic Data[J]. Acta Aeronautica et Astronautica Sinica, 2005, 26(4), 439-455 (in Chinese) |
| Cited By in Cnki (14) | Click to display the text | |
| [5] | Livne E, Terrence A. Weisshaar. Aeroelasticity of Nonconvertional Aircraft Configurations-Past and Future[J]. Journal of Aircraft, 2003, 40(6): 1047-1065 |
| Click to display the text | |
| [6] | 邱亚松, 白俊强, 李亚林, 等. 复杂几何细节对增升装置气动性能影响研究[J]. 航空学报, 2012, 33(3): 421-429 Qiu Y S, Bai J Q, Li Y L, et al. Study on Influence of Complex Geometry Details on the Aerodynamic Performance of High-Lift System[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(3): 421-429 (in Chinese) |
| Cited By in Cnki (5) | Click to display the text | |
| [7] | 白俊强, 刘南, 邱亚松, 等. 大型民用运输机短舱涡流片增升效率以及参数影响研究[J]. 西北工业大学学报, 2013, 31(4): 522-529 Bai J Q, Liu N, Qiu Y S, et al. Investigation on Influence of Nacelle Chine of Large Civil Transport Aircraft on High-Lift Efficiency and on Influence of Relevant Parameters[J]. Journal of Northwestern Polytechnical University, 2013, 31(4): 522-529 (in Chinese) |
| Cited By in Cnki (1) | Click to display the text | |
| [8] | John D Anderson, 吴颂平, 刘赵淼. 计算流体力学基础及其应用[M]. 北京:机械工业出版社, 2008 John Anderson, Wu S P, Liu Z M. Computational Fluid Dynamic Basic and Application[J]. Beijing, Mechanical Industry Press, 2008 (in Chinese) |
| [9] | Frank R. Scattered Data Interpolation: Tests of Some Methods[J]. Math Comp, 1982, 38(175): 191-200 |
| Click to display the text | |
| [10] | Ahrem R, Beckert A, Wendland H. A Meshless Spatial Coupling Scheme for Large-Scale Fluid-Structure Interpolation Problems[J]. Computer Modeling in Enginnering and Sciences, 2006, 12: 121-136 |
| Click to display the text | |
| [11] | Beckert A, Wendland H. Multivariate Interpolation for Fluid-Structure-Interaction Problems Using Radial Basis Functions[J]. Aerospace Science and Technology, 2001, 5: 125-134 |
| Click to display the text | |
| [12] | 刘艳, 白俊强, 华俊, 等. 非线性流固耦合分析方法研究[C]//中国力学大会2013, 2013 Liu Y, Bai J Q, Hua J, et al. An Approach to the Nonlinear Fluid-Structure Coupling Technical[C]//CCTAM2013, 2013 (in Chinese) |
| [13] | Allen C B, Rendall T C S. Unified Approach to CFD-CSD Interpolation and Mesh Motion Using Radial Basis Functions[R]. AIAA-2007-3804 |
| [14] | Buhmann M. Radial Basis Functions[M]. Cambridge, Cambridge University Press, 2005 |
| [15] | 白俊强, 刘南, 邱亚松, 等. 基于RBF动网格方法和改进粒子群优化算法的多段翼型优化[J]. 航空学报, 2013, 34(12): 2701-2715 Bai J Q, Liu N, Qiu Y S, et al. Optimization of Multi-Foil Basing on RBF Mesh Deformation Method and Modified Particle Swarm Optimization Algorithm[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(12): 2701-2715 (in Chinese) |
| Cited By in Cnki (4) | Click to display the text | |
| [16] | Rakowitz M. 2nd AIAA Drag Prediction Workshop DLR-F6 Geometry & Issues[EB/OL].(2003-09-01).http://www.awc.larc.nasu.gov/tsab/cfdlarc/aiau-dpw/workshop2/workshop2.html. |
| Click to display the text |
2. Chinese Aeronautical Establishment, Beijing 100012, China