改进的动弹网格方法在航空气弹计算中的应用

胡 凡 ,范锐军

(1.陕西省科技资源统筹中心，陕西 西安 710077)

(2.陕西西安现代控制技术研究所，陕西 西安 710065)

改进的动弹网格方法在航空气弹计算中的应用

胡 凡1,范锐军2

(1.陕西省科技资源统筹中心，陕西 西安 710077)

(2.陕西西安现代控制技术研究所，陕西 西安 710065)

气动弹性是现代航空气动力计算中一个突出的问题。主要研究基于Delaunay图映射方法的动弹网格的欧拉方程CFD计算及其在航空标模M6机翼上的静气动弹性应用。以Delaunay图映射方法为基础，针对三维非结构运动网格技术进行了研究、开发和改进，同时利用计算流体力学的方法, 开发了一套适用性较好的非结构网格欧拉方程流场求解器，进一步通过流固耦合的力学方法，对航空标模M6机翼的静气动弹性问题进行了研究和分析，给出了CFD并行计算的设计方法及算例。

Delaunay图映射动弹网格方法；欧拉方程；CFD/CSD耦合技术；静气动弹性

目前针对基于Euler方程或Navier-Stokes(N-S)方程等流动控制方程的CFD计算流固耦合问题，普遍采用有限体积的空间离散方法。在众多的网格方法中，非结构网格的应用非常广泛，它对于不同拓扑结构的复杂气动外形具有公认的普适性和灵活性。由于在航空气弹计算中，气动模型的外形随着气动弹性现象不断变化，网格边界上的网格点也需要根据气动模型的外形变化加以变形和调整，甚至适时进行重构，所以网格的运动技术成为了提高计算方法适应性的关键。

本文以航空标模M6机翼作为算例，通过自行开发的气动弹性程序进行静气动弹性的计算及分析。通过使用基于Delaunay图映射的动弹网格技术，成功地对多块混合非结构网格系统进行了运动变形，实现了三维非结构网格的合理运动。在此基础上，结合结构有限元分析，耦合结构静平衡方程，进行了数值模拟流固耦合计算。

1 动弹网格技术

随着计算模型几何外形复杂度的增加，网格的变形和重构、气弹耦合计算的效率以及计算程序的代数鲁棒性，是动弹网格技术中同时作用又相互制约的几个主要因素。针对大形变和多块混合网格(例如复杂外形的多块划分网格)的运动变形，现有的几种网格运动方法，例如TFI代数网格重构生成方法、弹簧近似方法Spring Analogy Method等[1]，都具有计算效率低、鲁棒性差等缺陷，而且在运动过程中还易出现网格的交叉和产生负体积[2]，严重限制了航空气动弹性求解方法的进步。

本文中作者自主开发并使用Delaunay图映射方法来进行动弹网格的构造和变形[3]，在一定程度上，保证了变形前后网格的拓扑结构和密度分布的一致性，同时也避免了网格交叉重叠现象的出现，运动变形之后的网格具有良好的网格品质，具备一定的高效性和良好的鲁棒性。

1.1动弹网格系数的确定

图1解释了动弹网格中一个重要的系数——面积比系数的确定方法。为了简化说明，本文采用二维情况，三维情况可以类推得到。

如图1中所示，网格节点G包含在三角形ABC中，面积比系数定义如下：

ei=Si/Si=1,2,3

式中：Si和S为三角形面积。

1.2动弹网格的计算方法

动弹网格的计算流程，主要包括以下几步[4]：

步骤1，将位于网格边界上的网格节点，定义为动弹网格特征点，并基于这些特征点来定义计算区域。对一组给定的网格特征点，总是存在着一组与之唯一对应的Delaunay图三角形，边界会在映射回来时自动保持，原始网格的拓扑结构会在移动和映射过程中保持不变。

步骤2，步骤1中最终生成的Delaunay图三角形集合与气弹计算的计算域能够完全重合，因此三角形区域能够包含网格中的任意一个节点。对于一个给定三角形单元内的网格节点，只有当它所有的面积比系数大于等于零时，该节点才包含于这个三角形单元。根据1.1节中定义的面积比系数，可建立Delaunay图三角形与任意网格节点的映射关系。此时，Delaunay图中的特征点能够真实准确地描述模型物面的边界运动。

步骤3，根据气动弹性计算中气动模型的形变规律，首先对与网格边界节点一一对应的Delaunay图三角形进行变形，然后再根据步骤2中建立的Delaunay图三角形单元顶点坐标与边界网格节点坐标之间的对应关系，更新变形后网格节点的坐标。注意：在此过程中，动弹网格的面积比系数不变。

如图2～图4所示，利用Delaunay图映射方法对大位移的混合网格进行变形(注：Delaunay图映射网格变形方法是非数值迭代方法)，变形后的效果较好，即在网格质量得到了较好的保证的同时，保证了计算效率和计算的稳定性。

2 非结构运动网格在静气弹求解中的应用

2.1基于非结构网格的流场求解

本文CFD解算控制方程是基于ALE描述的守恒型三维可压缩欧拉方程，其积分形式如下:

式中：Ω为控制体，∂Ω表示控制体单元的边界；dV是体积微元；n是控制体边界外法向单位向量；dS是面积微元；守恒变量项Q为

其中：ρ表示流体的密度;u,v,w分别为x,y和z轴方向的流体速度;e为单位体积流体的总内能。

对流通量矢量项F(Q)在x,y和z轴的分量为：

逆变速度矢量Ψ=Ui+Vj+Wk,沿x,y和z轴的分量U,V和W定义为：

式中：xt,yt和zt分别为网格沿x,y和z轴方向的网格速度。

此外对于理想气体压强p为：

式中：γ为理想气体比热比，这里取γ=1.4。

本文采用格心格式的有限体积法对CFD解算控制方程进行空间上的离散，采用经典四步龙格库塔方法进行虚拟时间上的离散，使用壁面函数对无粘求解进行了粘性修正，并且使用焓阻尼方法、当地时间步长、隐式残值光顺等加速收敛技术[5]，实现对三维气动弹性问题的流场求解。

2.2静气动弹性求解流程

本文使用CFD控制方程和结构静平衡方程耦合求解的方法对静气动弹性问题进行求解。利用气动力方程和结构静平衡方程耦合迭代求解，得到收敛的结果。

结构影响系数法是普遍适用于静气动弹性问题求解的一种方法[6]。结构影响系数法的主要原理描述如下。

结构网格节点(x,y)处的变形为

式中：W(x,y)为节点(x,y)处的法向变形；Czz(x,y,ξ,η)为机翼(ξ,η)处作用单位的力在节点(x,y)处发生的位移；F(ξ,η)为作用在(ξ,η)处的载荷，分为气动力和惯性力两部分，即

F(ξ,η)=L(ξ,η)-Nm(ξ,η)g

其中:L(ξ,η)为作用在(ξ,η)处的气动力，由CFD控制方程求解出的气动网格流场值进行三维面插值得到；N为过载系数；m(ξ,η)为(ξ,η)的质量；g为重力加速度。

静气动弹性计算流程如图5所示。

本文所采用的气弹计算算例，是跨音速航空标模M6机翼的静气动弹性变形计算。计算条件：马赫数Ma=0.84，初始攻角α0=3.06°，飞行高度20 000m。

静气动弹性的计算收敛结果如图6～7所示。从图6中可以看出，M6机翼刚轴沿展向均有由于静气动弹性所引起的垂向形变，这主要是由于机翼的升、阻力联合作用使得M6机翼产生了向上弯曲(根部固连于对称面)；从图7可以看出，由于气动弹性的作用，机翼在展向各个位置均产生了逐渐增大的负扭转角，这一负扭转角使展向各位置剖面的当地迎角减小，如图8和图9所示，机翼的升力分布在机翼产生静气动弹性变形后变化很大，同时机翼表面特别是机翼前、后缘的压力分布也发生了较大变化，弹性机翼相对于刚性机翼来说，展向各剖面的升力系数均有显著下降(主要由于机翼在气动力作用下发生弯曲和扭转变形，从而产生负的附加攻角)。静气动弹性变形给M6弹性机翼的跨音速升力特性带来了负面的影响。

图中，X/C为归一化当地翼型弦长，Y/b为归一化机翼展长。

3 并行化计算

通常航空气动弹性问题的计算，效率都比较低，需要耗费大量的计算时间和占用大量CPU[7]。本文使用MPI并行编程环境对气动弹性解算程序进行并行化程序改造，并使用4台单核8线程IBM工作站和8台单核8线程IBM工作站进行集群化并行计算，与8核CPU的IBM刀片服务器的并行计算精度以及并行效率对比，对比结果见表1。

其中， 定义网格并行加速比如下：假设tseq是使用一个求解单元计算所需的时间，tp是用p个求解单元计算所需时间；并行加速比定义为Sp=tseq/tp，并行效率定义为Ep=Sp/p。表1中的时间为计算程序进行一步时间迭代所需要的时间，单位为s。

由表1中可以看出，想要使用并行编程来获得较高的并行加速比和并行效率，需要从以下几个方面考虑：

a.采用良好的并行算法，在编程中最大限度地增加程序并行度，比如说采用显式方法来替代隐式方法。

b.减少通信规模，提高通信带宽，根据表1中的数据可以看出，集群计算的网络通讯速率要小于多核服务器内部系统总线的传输速率，从而使得网络通信成为集群并行计算实现高效计算的瓶颈。

4 结束语

本文所采用的Delaunay图映射的动弹性网格方法，能够大幅度提高动弹网格的变形能力，并在变形过程中较好地保证了网格质量和计算效率。本文通过对航空标模M6机翼的静气动弹性计算可以看出，该方法能够较好解决航空气动弹性计算中关键的网格变形问题，同时耦合气动力方程和静平衡方程的流固耦合结算方法也具有较高的结算精度和较好的适应性。本文所开发的并行化静气弹求解程序，能够较好地解决航空复杂构型的静气动弹性问题，在静气动弹性耦合数值模拟和气弹剪裁等方面，也拥有较为广阔的应用前景。

[1] Blom F J. Considerations on the spring analogy[J]．Journal of Aircraft, 2000,32: 647-668.

[2] Liu F, Cai J, Zhu Y．Calculation of wing flutter by a coupled CFD-CSD method [R]．AIAA-2000-0907,2000.

[3] 范锐军，冯朝辉．大展弦比无人机的静气弹问题计算及分析[J]．力学季刊， 2009(4)：182-188．

[4] 沈克扬，张锡华．弹性机翼的跨音速压强分布计算[J]．空气动力学学报，1984，1(3)：221～230．

[5] 王军利．基于非结构动网格的非定常气动力计算[D]．西安：西北工业大学，2004．

[6] 史爱明．非结构动网格下非定常气动力计算和嗡鸣研究[D]．西安：西北工业大学，2003．

[7] 范锐军，白俊强．气动弹性及并行化求解研究[J]．弹箭与制导学报，2006，6(1)：95-102．

Application of improved dynamic unstructured grids in aeroelastic model

HU Fan1, FAN Ruijun2

(1.Shaanxi Science and Technology Resource Center, Shaanxi Xi'an, 710077, China)

(2.Xi'an Modern Control Technology Research Institute, Shaanxi Xi'an, 710065, China)

Aeroelastic model is significant for large amount of airplanes in modern aerodynamics computing. This paper presents a strategy for generating 3D unstructured grids and the dynamic grids based on Delaunay graph mapping method. On the base of above, it performs a set of flow field solver based on Euler equations into establishmeng of the static aeroelastic cases of the M6 standard model coupled with structure dynamic equation, shows the design of the parallel computing method and numerical example.

dynamics grids deformation technique based on Delaunay graph mapping dynamics unstructured grids; Euler equations; CFD/CSD coupled technology; static aeroelastics

10.3969/j.issn.2095-509X.2015.03.013

2015-02-10

国家自然科学基金/青年基金项目资助(11302175)

胡凡(1975—)，女，陕西咸阳人，陕西省科技资源统筹中心工程师，主要研究方向为工业设计。

V211.3

A

2095-509X(2015)03-0060-05