在气动弹性和流固耦合这类问题的分析中,由于结构外形在外部气动力的作用下发生了变形,这就需要采用某种方法根据结构表面的变形使气动网格中的空间节点做出相应调整,即所谓的动网格方法[1]。如何提高动网格方法的计算效率和通用性正是目前研究的热点。
在动网格领域,常用的方法有无限插值法(TFI)和径向基函数法(RBF)。传统的TFI方法一般仅能用于结构化网格,但计算效率高,能够支持中等程度的变形。虽然国内学者对TFI方法的缺点进行了诸多的改进[1, 2],该方法也在很多问题中得到了应用[3, 4],但在三维算例中,要么弹身为方形截面[3, 4],要么根部无变形[1, 2]。即这些改进依然没有解决一般截面弹身组合体仅舵面变形的问题。最近出现的径向基函数法(RBF)在用于网格变形时,有着通用性好,容易编程实现和在变形连续时对原始网格的正交性保持性好的优点,从而在数值计算中得到了应用[5, 6, 7]。但是,文献[7]中仅考虑了舵面的刚性运动,尚未考虑舵面同时存在刚性运动和弹性运动时的动网格策略。
综上所述,目前常用的动网格技术都不能解决舵身组合体只有舵面运动,且同时发生刚性和弹性运动时的网格变形问题。其根本原因在于传统的动网格方法无法保持运动后的组合体几何外形:即使对于单个舵面这样的简单外形,由于实际计算过程中往往通过设置对称面边界进行半模计算以减少计算量,而原始的RBF和TFI方法无法保持对称面位置。如不引入新的方法,单纯的RBF方法、TFI方法或是RBF和TFI结合的方法必然会至少破坏一个部件(翼面／舵面／弹身／机身)的几何外形,甚至会将结合部位的各部件几何外形全部破坏掉。仅仅在极少数的情况下(弹身为方形截面或运动部件的根部不变形),才能保持组合体的几何外形。
为了解决该问题,本文采用有限元插值投影来保持运动后的组合体外形,并结合RBF和TFI动网格方法,形成保形动网格策略(SPGDM),然后将其应用于CFD／CSD耦合气动弹性分析系统中,计算了舵身组合体(弹身为圆形截面)中舵面在同时考虑刚性运动与弹性运动时的动响应。
1 保形动网格策略较新的RBF动网格方法[5, 6, 7]与传统的TFI动网格方法[1, 2, 3, 4, 8]已经在许多气动网格发生刚性／弹性变形的问题中得到了应用。但是这2种方法要进一步应用于复杂外形的气动网格变形则还存在各自的问题。
而本文提出的保形动网格策略通过将目前已有的成熟技术(基于有限元插值的投影技术,RBF插值技术和TFI动网格技术)有机地结合起来,以可以接受的计算效率,来解决组合体外形非定常CFD网格变形中的几何外形保持问题。针对原始物面变形不连续这一本质问题,采用合理的投影策略来获得连续的物面网格变形。其中，RBF插值的实现见文献[9]，径向基函数选择Wendland’s C2函数[9],采用贪心算法进行选点；TFI动网格的实现见文献[2]。
由于本文采用的基于有限元形函数的曲面插值有其特殊的处理步骤,故将其与本文提出的保形动网格策略分别进行介绍。
1.1 基于有限元形函数的曲面插值本文所采用的有限元插值,借鉴了文献[10]和文献[11]中等参单元逆变换的思想,采用数值迭代算法获得当地坐标,然后进行物理量的正向插值。
对一般的三维等参单元来说,坐标变换就是将各单元从整体坐标系x-y-z下变换到当地坐标系ξ-η-ζ下。其显式正变换为:
式中,xi,yi,zi为单元中第i个节点的位置坐标,n为单元的节点个数,Ni为第i个节点对应的形函数,x,y,z为插值点处的位置坐标,ξ,η,ζ为插值点处的当地坐标。一般物理量u拟采用与坐标变换相同的函数来进行插值:
式中,ui为单元中第i个节点处的物理量。
虽然(7)式～(9)式本身并不复杂,但是由于Ni一般并非ξ,η,ζ的线性函数,其逆变换一般并不存在解析形式。因此需要采用数值迭代算法进行逆变换。令b=(x,y,z)T,α=(ξ,η,ζ)T,则(7)式～(9)式可以表示为:b=f(α)。这是一个3阶非线性方程组,可以采用非线性方程组的迭代方法求解。本文选择牛顿迭代法来求解,具体求解步骤如下:
1) 给定迭代初值α0;
2) 计算第k步的误差向量:ek=f(αk)-b;
3) 如果‖ek‖足够小,则转步骤7);
4) 计算方程组雅克比矩阵:Jk= 
;
5) 解线性方程组:JkΔα=ek;
6) 更新迭代值:αk+1=αk+Δα,k=k+1,转步骤2);
7) 结束迭代:α=αk。
求出α(即:ξ,η,ζ)之后,将其代入(12)式就能获得u在另外一套网格节点上的值。
由于本文将上述方法用于三维曲面插值,此时α为二维,而b为三维,因此(5)式中需采用最小二乘法,即在等式两边都左乘JKT。此外为了提高曲面插值的精度和光滑性,宜采用2阶单元。
上面给出的逆变换求解方法对于插值点在网格单元内部或者距离并不远的情况下,收敛非常迅速,对本文使用的算例,一般10步之内就能收敛。但是,如果插值点与网格单元相距较远时就不容易收敛。所以,为了保证逆变换求解的正确性和效率,应先确定插值点所在网格单元,具体方法见文献[12]。
1.2 保形动网格策略的实现步骤在实现了上文的RBF技术、基于有限元形函数的插值和TFI动网格方法之后,本文提出的保形动网格策略具体步骤和示意图如下(实线表示CFD网格的物面边界,虚线表示投影基面网格,为保证相贯,投影基面应有适当延伸):
 |
| 图 1 原始CFD物面 |
 |
| 图 2 原始投影基面 |
1) 采用RBF方法使投影基面网格发生弹性变形,然后使其发生刚性变形,如图 3所示;
 |
| 图 3 投影基面变形 |
2) 采用基于有限元插值的投影方法将上一步的物面网格的边界向各投影基面网格反复投影,直到相贯线位置不再变化,如图 4所示;
 |
| 图 4 计算相贯 |
3) 采用TFI技术将物面网格的边界变形均匀地插到整个物面上,如图 5所示;
 |
| 图 5 物面网格插值 |
4) 采用基于有限元插值的投影方法将物面的内点向各投影基面网格投影获得当前步的物面网格,如图 6所示。
 |
| 图 6 物面网格投影 |
5) 采用TFI技术将物面网格的变形均匀地插往全场网格,获得当前步的全场网格。
本文实现的保形动网格策略的实际应用效果见3.2.2节。
2 CFD／CSD耦合方法 2.1 气动力计算三维可压缩非定常守恒积分形式的欧拉方程[8]为:
式中,t为时间,Q为守恒变量矢量,Fc为对流通量项,Ω和∂Ω是控制体和其边界。对(1)式采用有限性体积法进行求解;对流通量项Fc采用Vanleer格式计算,该格式在超声速流域的稳定性较好;限制器选用鲁棒性较好的minmod型限制器;时间推进格式采用双时间推进,通过引入子迭代将LUSGS格式的时间离散精度扩展到2阶;边界条件在远场处根据速度方向判断是流入还是流出,流入时赋远场值,流出时由内场做一次外插;在物面处对速度采用无穿透条件,令温度和压强在物面法向梯度为零。 2.2 结构动力学求解
气动弹性耦合系统中的结构动力学部分采用模态叠加法求解。模态叠加法本质是将各阶结构模态解耦后的降阶方法。利用模态解耦后的结构动力学方程为:
式中:
当存在刚性扰动时,将刚性运动产生的惯性力作为外力加到右端项P中,然后对(20)式采用4阶龙格库塔方法进行时间推进求解。
2.3 耦合技术
气动弹性耦合计算采用松耦合策略,如图 7所示,主要步骤有:
 |
| 图 7 气动弹性耦合示意图 |
1) 气动力计算;
2) 将气动力(压强)插往结构网格;
3) 进行结构动力学计算;
4) 将结构变形插往气动网格,并进行网格变形。
其中,第2步和第4步的插值是在流固耦合边界上进行的(三维问题就是面插值)。在气动系统方面,耦合边界上的速度在边界法向上为零;在结构动力学方面,耦合边界上采用给定压强的边界条件。气动力的插值采用RBF技术,气动网格变形采用本文提出的保形动网格策略予以实现。
3 数值计算为了考察保形动网格策略的保形效果及在CFD／CSD耦合计算中的适用性,以圆截面弹身的舵身组合体为对象进行了16°舵偏和扰动下的动响应计算。其中,气动舵为全动舵,通过舵轴与弹身相连。在对舵身组合体进行CFD／CSD耦合计算时,将弹身考虑为刚性,仅考虑气动舵本身的弹性效应。
3.1 模型描述 3.1.1 几何外形及气动网格气动舵的翼型为双弧形,相对厚度为0.05;平面形状为梯形,根弦弦长0.8 m,梢弦弦长0.36 m,前缘后掠角34.02°。舵身组合体为半模,图 8中给出了组合体的的网格拓扑和物面网格。全场网格共21个块,单元总数为373 480。
 |
| 图 8 舵身组合体气动网格 |
气动舵的结构模型如图 9所示。该气动舵由舵面和舵轴组成,舵面和舵轴的材料参数见表 1,舵面采用实体单元进行建模,舵轴采用半径为10 mm的圆截面梁进行建模。舵轴位于舵面根部1／2弦线处。在舵轴根部加载固支约束,在舵轴与舵面相连的节点通过RBE2与周围8个点相连。
 |
| 图 9 气动舵结构网格 |
为了在舵身组合体的耦合计算中保证弹身和舵面的几何外形都不被破坏,本文采用了保形动网格策略,其中所需的投影基面如图 10所示。
 |
| 图 10 投影基面网格 |
采用Nastran计算获得气动舵的前4阶结构模态,其中:1阶模态为舵轴弯曲模态,2阶模态为舵轴扭转模态,3阶模态为舵面弯曲模态,4阶模态为舵面扭转模态。从图 11中可以看出,舵轴模态的频率要比舵面模态的频率低得多:舵轴模态频率不超过60 Hz,舵面模态的频率不低于299 Hz。这说明舵轴刚度要远小于舵面刚度。
 |
| 图 11 舵面结构模态 |
对舵面施加16°的舵偏,采用常规TFI方法,保持舵面形状不变,获得变形后的网格如图 12所示。
 |
| 图 12 TFI方法所获得的组合体变形后网格 |
 |
| 图 13 保形动网格策略获得的组合体变形后网格 |
从图 12中可以看到,由于舵面变形的拖拽效应,导致弹身发生了巨大的变形,而这种变形在实际导弹舵偏时完全不存在。由表 2知变形后的弹身网格仅有7.06％的网格点位于真实的弹身上。采用图 12中的气动网格进行定常CFD计算(Ma=2,H=3 km)所获得的压力系数如图 14所示,从图 14可以看出,由于弹身被拉坏,产生了非物理的压缩区和膨胀区。可想而知,在组合体的CFD／CSD耦合分析中采用常规TFI方法必然会造成巨大的误差,结果的可信度将非常有限。
 |
| 图 14 TFI变形后Cp云图 |
而采用保形动网格策略,在舵偏16°时,可同时保持舵面和弹身形状不变,变形后的网格如图 13所示。在这一过程中所使用的4个投影基面网格如图 10所示,为了保证投影的精度,这些投影基面都采用了8节点四边形单元。为了考察保形动网格策略是否有效,表 2给出了弹身节点的保形效果,从表中可以看出气动网格弹身节点在变形之后100％保持在了弹身上。从图 13中可以清晰地看出,采用了保形动网格策略之后,在常规TFI方法中巨大的弹身变形得以完全避免。从图 15中可以看出由于变形之后的弹身网格依然维持了原始外形,因此没有产生非物理的压缩区和膨胀区,从而说明保形动网格策略能够保持组合体的CFD／CSD耦合分析的可信度。
 |
| 图 15 SPGDM变形后Cp云图 |
采用本文的CFD／CSD耦合方法,在Ma=2,H=3 km的来流状态下,对舵身组合体进行扰动下的动响应分析。CSD部分采用3.2.1节中计算出的前4阶结构模态。耦合时间步长取为0.0001 s,内迭代步数为50步。气动舵上加载的舵偏扰动为:δ(t)=Δδasin(2πωt)。其中,Δδa=2°,ω=40 Hz。
表 3中给出了在1个非定常耦合时间步内,不同动网格方法计算时间(计算机CPU主频3.4GHz,内存2G)的对比。
从表 3中可以看出,虽然保形动网格策略本身的计算时间(0.73 s)相对TFI方法(0.05 s)大大增加,但在单个CFD／CSD耦合时间步内仅使计算时间(32.92 s)增加了0.68 s,不到总时间的2.1％。考虑到保形策略能保证耦合分析的可信度,仅增加2.1％的计算时间完全可以接受。
为了考察该种扰动下的气动弹性效应,特别设置2个算例:case 1中来流状态为Ma=2,H=3 km,case 2则忽略外界气流的影响,仅考虑结构本身在扰动激励下的响应。两算例各阶模态广义位移的对比如图 16所示。
 |
| 图 16 广义位移响应对比 |
从该图中可以看出,case 1和case 2的响应主要位于1阶模态和2阶模态上;但由于气动和结构的耦合效应,case 1的广义位移振荡幅值远超case 2,并且幅值随时间衰减,在达到动稳定之后,1阶模态和2阶模态的振荡幅值分别是case 2的3.2倍和1.75倍。
为了考察保形动网格策略在组合变形下的效果,特别选择case 1中t=0.105 s(N=1 050,刚性舵偏角为1.8°,各阶广义位移为0.238,0.073,-0.591×10-3,0.174×10-4)和case 2中t=0.356 s(N=3 560,刚性舵偏角为2°,各阶广义位移分别为-0.022,0.043,-0.102×10-6,-0.155×10-5)这2个典型时刻,其物面网格变形如图 11所示。从图 17中可以看出,在气动弹性耦合效应下,case 1的变形远大于case 2;即使在case 1如此巨大的变形之下,保形动网格策略依然保证了弹身和舵面的几何形状不被破坏,同时还保持了原始网格上较为合理的网格分布。
 |
| 图 17 典型时刻物面网格变形示意图 |
本文结合RBF技术、基于有限元形函数的插值和TFI动网格方法,形成了能够同时考虑刚性变形与弹性变形的保形动网格策略并将其应用于舵身组合体CFD／CSD耦合计算中。通过计算分析,获得结论如下:
1) 通过使舵身组合体刚性舵偏16°,验证了本文提出的保形动网格策略能够避免常规TFI方法导致的巨大弹身变形,实现了100％的保形率,从而没有产生非物理的压缩区和膨胀区,说明保形动网格策略能够保证组合体的CFD／CSD耦合分析的可信度。
2) 在舵身组合体CFD／CSD耦合计算中发现,采用保形动网格策略仅使单个非定常步内总时间(32.92 s)增加了0.68 s,不到总时间的2.1％,其计算效率可以接受。
3) 在Ma=2,H=3 km的来流和振幅为2°频率为40 Hz舵偏扰动的共同作用下,舵身组合体气动弹性效应明显,其1阶模态和2阶模态的广义位移幅值为无来流时的3.2倍和1.75倍。
(4) 在CFD／CSD耦合分析中产生的组合变形下,保形动网格策略依然保证了弹身和舵面的几何形状不被破坏,同时还保持了原始网格上较为合理的网格分布。
| [1] | 黄礼铿, 高正红, 左英桃. 一种快速稳健的并行多块结构动网格方法[J]. 计算力学学报, 2012,29(3):363-367 Huang L K, Gao Z H, Zuo Y T. A Fast and Robust Paralelizable Moving Mesh Algorithm for Multi-Block Structured Grids[J]. Chinese Journal of Computational Mechanics, 2012,29(3):363-367 (in Chinese) |
| Cited By in Cnki (5) | Click to display the text | |
| [2] | 张兵, 韩景龙. 带旋转修正的弹簧-TFI混合动网格方法[J]. 航空学报, 2011, 32(10):1815-1823 Zhang Bing, Han Jinglong. Spring-TFI Hybrid Dynamic Mesh Method with Rotation Correction[J]. Acta Aeronoutica et Astronautica Sinica, 2011, 32(10):1815-1823 (in Chinese) |
| Cited By in Cnki (4) | Click to display the text | |
| [3] | 陈坚强, 陈琦, 袁先旭,等. 舵面操纵动态响应的数值模拟研究[J]. 力学学报, 2013(2):302-306 Chen Jianqiang, Chen Qi, Yuan Xianxu, et al. Numerical Simulation Study on Dynamic Response under Rudder Control[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013(2):302-306 (in Chinese) |
| Cited By in Cnki (2) | Click to display the text | |
| [4] | 董琳琳. 超限插值法在舵面振荡偏转数值模拟中的应用[J]. 弹箭与制导学报, 2011(4):207-209 Dong Linlin. The Application of Transfinite Interpolation to Numerical Simulation of Control Surface Oscillating Deflection[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2011(4):207-209 (in Chinese) |
| Cited By in Cnki (2) | Click to display the text | |
| [5] | 谢亮, 徐敏, 张斌,等. 基于径向基函数的高效网格变形算法研究[J]. 振动与冲击, 2013, 32(10):141-145 Xie Liang, Xu Min, Zhang Bin, et al. Space Points Reduction in Grid Deforming Method Based on Radial Basis Functions[J]. Journal of Vibration and Shock, 2013, 32(10):141-145 (in Chinese) |
| Cited By in Cnki (1) | Click to display the text | |
| [6] | 谢亮, 徐敏, 安效民,等. 基于径向基函数的网格变形及非线性气动弹性时域仿真研究[J]. 航空学报, 2013, 34(7):1501-1511 Xie L, Xu M, An X M, et al. Research of Mesh Deforming Method Based on Radial Basis Functions and Nonlinear Aeroelastic Simulation[J]. Acta Aeronoutica et Astronautica Sinica, 2013, 34(7):1501-1511 (in Chinese) |
| Cited By in Cnki (5) | Click to display the text | |
| [7] | Andreas K M. Aircraft Control Surface Deflection Using RBF-Based Mesh Deformation[J]. International Journal for Numerical Methods in Engineering, 2011, 88:986-1007 |
| Click to display the text | |
| [8] | 谢亮, 徐敏, 李杰,等. 基于CFD/CSD耦合的颤振与动载荷分析方法[J]. 振动与冲击, 2012, 31(3):106-110 Xie L, Xu M, Li J, et al. Flutter and Dynamic Analysis Based on CFD/CSD Coupling Method[J]. Journal of Vibration and Shock, 2012, 31(3):106-110 (in Chinese) |
| Cited By in Cnki (3) | Click to display the text | |
| [9] | 王刚,雷博琪,叶正寅. 一种基于径向基函数的非结构混合网格变形技术[J]. 西北工业大学学报, 2011, 29(5):783-788 Wang G, Lei B Q, Ye Z Y. An Efficient Deformation Technique for Hybrid Unstructured Grid Using Radial Basis Functions[J]. Journal of Northwestern Polytechnical University, 2011, 29(5):783-788 (in Chinese) |
| Cited By in Cnki (7) | Click to display the text | |
| [10] | 包劲青, 杨强, 官福海. 线性等参单元逆变换的几何二分法及其修正[J]. 计算力学学报, 2010, 27(5):770-774 Bao J Q, Yang Q, Guan F H. Geometric Dichotomy Method for Inverse Isoparametric Mapping in Linear Finite Element and Its Correction Measure[J]. Chinese Journal of Computational Mechanics, 2010, 27(5):770-774 (in Chinese) |
| Cited By in Cnki (1) | Click to display the text | |
| [11] | 周艳国, 陈胜宏, 张雄等. 改进的等参逆变换算法在耦合场分析中的应用[J]. 岩土力学, 2008, 29(11):3170-3173 Zhou Y G, Chen S H, Zhang X, et al. Improved Inverse Isoparametric Mapping Method and Its Application to Coupled Analysis[J]. Rock and Soil Mechanics, 2008, 29(11):3170-3173 (in Chinese) |
| Cited By in Cnki (3) | Click to display the text | |
| [12] | Wang J Z, Parthasarathy V. A Fully Automated Chimera Methodology for Multiple Moving Body Problems[J]. International Journal for Numerial Methods in Fluids, 2000, 33:919-938 |
| Click to display the text |