基于网格速度法的非定常流场模拟和动导数计算

基于网格速度法的非定常流场模拟和动导数计算
郭东;徐敏;陈士橹
【摘要】Sections 1 through 3 of the full paper explain and evaluate the computation method mentioned in the title, which we believe is effective. Their core consists of; the field velocity or grid velocity approach provides a u-nique methodology for directly calculating aerodynamic responses to step change in flow conditions; the grid time metrics include the velocity caused by the impulsive change in angle of attack but the mesh is not moved accordingly; this approach avoids numerical instabilities and decouples the step change in the angle of attack from a pitch rate; based on this approach, a technique is presented to model longitudinal unsteady flow phenomenon by superposing the step change in the angle of attack upon the impulsive change in pitch rate. In Figs. 4 and 5, numerical results are validated by comparison with experimental results for NACA 0012 airfoil under forced oscillations. To validate further the applicability for the present method, pitch damping derivatives, calculated from the load history of the unsteady flow around a standard research configuration, known as the Basic Finner Missile, are presented in Figs. 10 and 11. Predicted results show indeed good agreement with available wind tunnel data.%基于网格速度法的思想,将迎角的突变和俯仰角速率的突变叠加起来,发展了一套在固定网格上模拟飞行器俯仰振荡非定常流场的方法,该方法不需要实时更新网格,减少了计算时间和所需内存,避免了负体积的出现.首先计算了NACA0006的阵风响应和NACA0012翼型的俯仰振荡,所得结果与实验值和动网格方法符合较好,
这表明该方法能够准确模拟此类非定常问题；最后将该方法应用于国外动导数计算标模Basic Finner Missile(BFM)俯仰振荡运动的数值模拟,并计算了其在马赫数
1.58 ～
2.5的静、动稳定性导数,计算结果与风洞实验值基本吻合,体现了该方法的
正确性.
【期刊名称】《西北工业大学学报》
【年(卷),期】2012(030)005
【总页数】5页(P784-788)
【关键词】计算流体力学;非定常流场;网格速度法;动导数
【作者】郭东;徐敏;陈士橹
【作者单位】西北工业大学航天学院,陕西西安710072;西北工业大学航天学院,陕
西西安710072;西北工业大学航天学院,陕西西安710072
【正文语种】中文
【中图分类】V211.3
动导数是飞行器动态气动特性设计中的关键参数。

精确导数计算工具是飞行器设计部门所必需的。

随着计算机技术和计算流体力学(CFD)的快速发展，CFD技术已经成为求解飞行器的静、动导数的一种重要手段［1-3］。

非定常流场的准确求解是气动迟滞效应时域仿真的关键，直接影响到静、动导数的计算精度。

目前，利用CFD求解静、动导数时，非定常流场都是采用刚性或者柔性动网格方法来实现的，然而动网格方法需要在每个时间步更新网格，这不仅增加了计算时间和所需内存，还可能出现负体积。

为了求解飞行器迎角或者俯仰角速率的阶跃响应，Baeder等人提出了“网格速度法”，将迎角的阶跃变化，采用对整个流场区域施加一个垂直方向的速度来进行模拟。

在理论上，这只改变了迎角，而对俯仰角速率没有任何影响。

对于角速率的阶跃响应，可以通过对网格点施加一个与旋转中心距离呈线性关系的网格速度来描述，即v(x，y)=˙α·r(x，y)，其中v(x，y)是网格点对应的速度，r(x，y)是网格点相对
于旋转中心的距离，从而实现了迎角和俯仰角速率时间历程的解耦，由此得到的结果与精确值吻合较好，而且还不会产生数值振荡［4，5］。

该方法被广泛应用于
分析飞行器的阵风响应特性［6-8］。

本文在“网格速度法”的基础上，将迎角的突变和俯仰角速率的突变叠加起来，发展了一套在固定网格上模拟飞行器俯仰振荡非定常流场的方法。

相对“动网格方法”而言，该方法具有以下特点:1)不用实时更新网格，从而减少了计算时间和所需内存;2)不需要考虑网格的密度、类型，避免了负体积的出现;3)可以研究单纯的迎角
或者俯仰角速率阶跃响应和阵风响应。

首先通过研究NACA 0006的阵风响应来
验证了本文所使用程序的可靠性，然后对NACA 0012翼型俯仰振荡进行了数值
研究，通过CT5算例与实验值和动网格方法所得结果进行比较，结果吻合较好，
从而表明了该方法和使用的CFD程序能够比较准确模拟此类非定常问题;最后利用该方法数值模拟了国外动导数计算标模Basic Finner Missile(BFM)的非定常俯仰
振荡，并计算了其在马赫数1.58˜2.5的静、动导数，计算结果和风洞实验结果吻
合较好，体现了方法的正确性。

1 理论基础及数值方法
1.1 流场控制方程
采用有限体积法离散Euler方程
式中，Q=(ρ，ρu，ρv，ρw，ρe)T 为解向量;∂Ω 为流体域Ω的边界;n为边界的
外法向矢量，F为对流矢通量。

无粘通量利用Roe通量差分分裂格式离散，采用双时间推进，引入子迭代，利用二阶LUSGS隐式格式推进。

在计算过程中，严格地满足几何守恒定律［9］。

1.2 网格速度法
在计算区域内的流场速度可以表示为
式中，u、v、w分别是来流速度沿坐标轴的分量，xτ、yτ、zτ是相应的网格速度分量。

当流场中还有外来的速度场(u'，v'，w')时，流场速度变为
网格速度方法通过改变网格的速度来模拟这些速度场。

修正后的网格速度可以表示为
从而，整个流场的速度为
1.3 俯仰振荡模拟方法
对于俯仰振荡运动，本文将其视为迎角的突变和角速率的突变两个部分的叠加。

迎角的突变可以用流场在z方向的速度突增来模拟。

此时相当于在流场中注入了一个外部速度场，其3个速度分量分别为
式中，α0是初始迎角，Δα是迎角的增量。

对于角速率的变化，可以通过对网格点施加一个与旋转中心距离呈线性关系的网格速度来模拟
综上，结合方程(5)和(6)即可得到俯仰振荡运动中网格速度为
式中，(xi，yi，zi)是网格点的坐标，(xc，yc，zc)是转轴中心。

基于此，就可以在不使用动网格技术的前提下，模拟飞行器的俯仰振荡运动。

1.4 静、动态稳定性参数辨识
对单自由度强迫俯仰振荡，给定简谐振动(无量纲化)形式如下
式中，α0为初始迎角，Δα 为振幅，k= ωLref/(2V∞)为减缩频率，t*是无量纲时间，对动态俯仰力矩系数在起始攻角处展开(保留至一阶导数项)如下
式中，Cm0为起始攻角处的静态俯仰力矩系数为俯仰静导数为俯仰动导数。

参考文献［3］，利用积分法有
2 算例验证
为了验证本文的算法和所使用CFD程序的可靠性，首先计算了有实验结果的非定常气动力标准算例。

2.1 NACA 0006翼型迎角的阶跃响应
对于迎角的阶跃变化，采用对整个流场区域施加一个垂直方向的速度来进行模拟。

假设翼型初始状态分别为 Ma=0.3，0.5，0.65，0.8，初始迎角为0°，受到垂直阵风的作用，取wg=0.08V∞，导致翼型迎角突然增加约4.6°。

取真实时间步长dt=1.0×10-4s，定义无量纲时间，其中c是参考长度。

图1给出了不同马赫数在此阵风作用下翼型的升力系数随无量纲时间的变化历程。

Lomax［10］推导了线性可压缩流中平板翼型迎角突增α时，阵风响应初始阶段升力系数的理论解
图1 迎角阶跃响应下的升力响应
图2给出了该理论解与本文计算结果之间的比较，从图中可以看出，在Ma=0.3
时，计算结果与理论解吻合得较好，随着马赫数增加，计算结果与理论解的偏差增大，其原因是由于随着马赫数增加，流场的非线性特性增加造成的。

2.2 NACA0012强迫俯仰振荡
AGARD CT5是常用的非定常算例，其运动规律为α = α0+ Δαsin(2kt*)，
α0=0.016°，Δα =2.51°，k=0.0814翼型绕1/4弦长处做俯仰运动，图3给出了AGARD CT5算例［11］验证所采用的C型网格，网格维数为199×40。

图4给
出了升力和力矩系数随时间的变化曲线，图5为不同迎角时翼型上下表面的压强
分布，从图中可以发现，计算结果与实验值［12］和动网格方法结果吻合较好，
动网格采用刚性旋转法，这表明本文提出的方法和采用的仿真程序在计算此类非定常问题时具有较高的可靠性。

图2 升力响应计算结果与理论解比较
图3 CT5算例网格示意图
图4 AGARD CT5算例计算结果与实验数据和动网格方法结果比较
图5 不同迎角状态物面CP曲线
3 BFM静、动导数计算研究
Basic Finner Missile(BFM)是国际上用于动导数计算程序验证标准模型［2］，其外形尺寸和计算网格如图6和图7所示。

重心位置为5.0d，Ma=1.58，1.75，1.89，2.10，2.50，减缩频率为 k=0.003 16，初始迎角α0=1.5°，振幅Δα
=1.5°。

图8和图9分别给出了Ma=1.58时俯仰力矩系数的响应曲线和迟滞曲线，所得结果与动网格方法吻合较好。

图10和图11分别给出了不同马赫数下静导数和动导数与实验数据的比较。

从中
可以看到，计算结果和实验数据相吻合，计算结果略小于实验值，静导数最大误差10.089%，动导数最大误差4.335%，这说明了本文方法的可信度。

图6 BFM外形尺寸
图7 BFM计算网格示意图
图8 Ma=1.58时俯仰力矩系数的时间响应
图9 Ma=1.58时俯仰力矩迟滞曲线
图10 不同马赫数下静导数比较
图11 不同马赫数下动导数比较
4 结论
本文基于“网格速度法”的思想，将迎角的突变和俯仰角速率的突变叠加起来，发展了一套在固定网格上模拟飞行器俯仰振荡非定常流场的方法，研究后得到以下结论:
1)通过研究NACA 0006的阵风响应和NACA 0012翼型俯仰振荡运动，该方法所得结果与实验值和“动网格方法”所得结果基本一致，这表明该方法可以较为准确地模拟此类非定常流场。

2)通过数值模拟国外动导数计算标模BFM的跨音速非定常纵向俯仰振荡，将计算得到的静、动导数与风洞实验结果进行比较，所得结果与风洞实验值吻合较好，体现了该方法的正确性。

3)相对“动网格方法”而言，该方法具有以下特点:①不用实时更新网格，从而减少了计算时间和所需内存;②不需要考虑网格的密度、类型，避免了负体积的出现;③可以研究单纯的迎角或者俯仰角速率阶跃响应和阵风响应。

参考文献:
［1］牟斌，刘伟，霍章华.球椎体高超声速扰流的俯仰阻尼导数的数值计算.国防科技大学学报，2000，22(4):5-8 Mou Bin，Liu Wei，Huo
Zhanghua.Numerical Calculation of Damping-in-Pitch Derivatives for Hypersonic Flow over Sphere-Core.Journal of National University of Defense Technology，2000，22(4):5-8(in Chinese)
［2］Oktay Erdal，Akay Hasan U.CFD Predictions of Dynamic Derivatives for Missiles.AIAA-2002-0276
［3］袁先旭，张涵信，谢昱飞.基于CFD方法的俯仰静、动导数数值计算.空气动力学学报，2005，23(4):458-463 Yuan Xianxu，Zhang Hanxin，Xie Yufei.The Pitching Static/Dynamic Derivatives Computation Based on CFD Methods.Acta Aerodynamic Sinica，2005，23(4):458-463(in Chinese) ［4］Parameswaran Vasudev，Baeder James D.Indicial Aerodynamics in Compressible Flow-Direct Computational Fluid Dynamic Calculations.Journal of Aircraft，1997，34(1):131-133
［5］Rajneesh Singh，James D Baeder.Direct Calculation of Three-Dimensional Indicial Lift Response Using Computational Fluid Dynamics.Journal of Aircraft，1997，34(4):465-471
［6］詹浩，钱炜祺.薄翼型阵风响应的数值模拟.航空学报，2007，28(3):527-530 Zhan Hao，Qian Weiqi.Numerical Simulation of Gust Response for Thin Airfoil.Acta Aeronautica et Astronautica Sinica.2007，28(3):527-530(in Chinese)
［7］Yang Guowei，Obayashi Shigeru.Numerical Analyses of Discrete Gust Response for an Aircraft.Journal of Aircraft，2004，41(6):1353-1359 ［8］詹浩，钱炜祺.弹性机翼阵风响应数值计算方法.计算力学学报，2009，
26(2):270-275 Zhan Hao，Qian Weiqi.Numerical Simulation on Gust Response of Elastic Wing.Chinese Journal of Computational Mechanics，2009，26(2):270-275(in Chinese)
［9］Sitaraman Jayanarayanan，Baeder James D.Field Velocity Approach and Geometric Conservation Law for Unsteady Flow Simulations.AIAA
Journal，2006，44(9):2084-2094
［10］Lomax H T.Indicial Aerodynamics.AGARD Manual of Aeroelasticity:Part 2，NATO.Advisory Group for Aerospace Research and Development，1968，Chap.6
［11］Li Jie，Huang Shouzhi，Jiang Shengju，Li Fengwei.Unsteady Viscous Flow Simulations by a Fully Implicit Method with Deforming Grid.AIAA-2005-1221
［12］Landon R H.NACA0012 Oscillatory and Transient Pitching.AGARD Report 702，Compendium of Unsteady Aerodynamic Measurements，1982。