随着计算机运算速度的大幅提高和计算流体力学(Computational Fluid Dynamics，CFD) 方法的日趋完善[1]，快速求解Navier-Stokes方程已成为可能但快速的工程方法在飞行器设计中仍有很大的应用前景。无黏Euler方程CFD方法与壁面黏性力快速工程计算方法相结合能够快速预测流场特性，并保证良好的计算精度[2]。因此目前工程应用中，特别是针对复杂飞行器外形绕流情况，常常采用该方法快速获得满足工程设计精度要求的气动性能[3-4]。
壁面黏性力工程计算方法通常基于边界层理论，以自由来流参数作为边界层外缘参数进行计算，并需要针对飞行器几何形状进行简化和近似。比如已经获得广泛应用的等效平板层流Blassius近似方法和湍流Van Driest Ⅱ方法[5]，工程应用于一定马赫数和雷诺数范围内均具有较好的预测精度[6]。而对于高超声速高温气体作用下的壁面黏性力预测，参考温度法则提供了一种基本满足工程要求的预测能力[7-8]。然而，对于复杂的飞行器外形，特别是在高超声速情况下，自由来流受到飞行器各部件的干扰非常严重，使得壁面边界层外缘的流场参数与自由来流参数相差很大，因而上述壁面黏性力工程计算方法从理论上就会产生较大的偏差。而事实上，基于无黏Euler方程的CFD方法在进行无黏流场计算时，就已经比较准确的获得了边界层外缘的流场参数。鉴于此，本文从无黏Euler方程CFD方法获得的壁面流场参数出发，建立了基于当地流线长度和流场参数的参考温度法，应用于高超声速飞行器壁面黏性力的预测。将本文计算方法获取的壁面黏性力计算结果与完全Navier-Stokes方程CFD方法计算结果进行了比较，分析表明该方法相比于基于自由来流参数的经典参考温度计算方法具有更高的精度，能够更好的满足工程应用需求。
1 基于无黏流场的壁面流线计算方法无黏Euler方程计算中，壁面采用无穿透边界条件，即速度方向平行于壁面。通过壁面的三个速度分量对时间的积分得到壁面流线上的坐标值。假设P=(x, y, z)T为流线上的点，该点上的速度为：
|
(1) |
根据流线[9]的定义：
|
(2) |
初值条件为：x(t0)=x0, y(t0)=y0, z(t0)=z0
采用4阶显示Runge-Kutta方法求解该初值问题，具体公式如下：
|
(3) |
其中
|
(4) |
h为步长，通常可以取一个较小的常数，反映时间推进的间隔。
由于壁面流线的特殊性，流线计算推进时需要解决一些特殊问题[10]，主要包括：时间步长的合理选取；流线推进点在壁面的投影；壁面速度场的插值。下面分别进行讨论。
1.1 变时间步长处理方法时间步长h的选取决定了流线计算的精度和效率。如果时间步长取值过大，必然造成流线计算的偏差很大；如果时间步长取值过小，流线计算精度虽然很高，但会极大的增加计算时间。通常，Runge-Kutta方法中时间步长取某一合理的常数，但其合理性也只是从整个流场大概来说的。另一方面，无黏流场计算获得的壁面流动参数与来流状态和气动外形密切相关，因此合理的时间步长还需要根据具体的流动特性来确定。鉴于此，本文提出了流场自适应的变时间步长处理方法。
流场自适应的变时间步长处理方法基本思路是：从流线点Pn推进到下一个点Pn+1时，4阶Runge-Kutta方法中的时间步长h取点Pn所在当地单元流动特征时间Tc。该特征时间定义为单元四条边中以各自平均速度通过的最小时间(见图 1)，即:
|
(5) |
其中
|
eij为点i到点j的单位矢量。
|
| 图 1 当地单元特征时间定义示意图 Fig. 1 Schematic diagram of specific time at the local cell |
该时间基本反映了流动穿越该当地单元的特征时间，因此能够保证每一个壁面单元都至少存在一个流线点，从而保证了流线精度；同时也保证同一个单元内不会存在过多的流线点使得流线计算时间太长。流场自适应的变时间步长处理方法有效解决了流线计算精度和计算效率的矛盾。
1.2 流线推进点在壁面的投影壁面流线与空间流线的很大差别就是其二维特性，即流线点总是在壁面上，实质是近似流线，而并非完全精确的流线。尽管可以保证壁面上各点的速度方向平行于壁面切向，但壁面的不光滑必然使得壁面所有点的速度方向不会在同一个平面上，这就导致流线沿着速度方向推进时必然要向空间发展，而不会再位于壁面上：如图 2所示，流线点Pn沿其速度方向V(平行于该点所在壁面单元) 推进到点Pn+1*，该点已经不在任何壁面单元内，而是在空间某个位置。因此要得到壁面流线，必须将推进到空间的流线点Pn+1*投影到壁面上以获得Pn+1。
|
| 图 2 壁面流线推进示意图 Fig. 2 Schematic diagram of streamline tracing on the wall surface |
然而空间点投影到壁面也存在很多细节上需要处理的问题。对于壁面四边形结构网格，其单元由四个节点构成(图 3)，而此四个节点并非一定共面，因此空间点往四边形单元的投影存在不确定性。而不共线的三点能够构成一个平面，因此可以将四边形单元拆解为两个三角形平面单元(①和②)；于是空间点P*往壁面的投影可以简化为空间点往壁面所有三角形平面单元的投影，取其中投影距离最小的点为真实投影点。如果空间点Pn+1*在所有壁面单元内都无法找到投影点，则认为流线在点Pn终止。
|
| 图 3 空间点投影到壁面单元的示意图 Fig. 3 Schematic diagram of spatial point projecting on the surface cell |
1.3 壁面流场参数的插值
流线推进获得壁面投影点的坐标后，其该点的流场参数是未知的。而已知的仅仅是单元节点上的值。因此需要采用插值方法获得单元内投影点的流场参数。本文采用了反距离权重的插值方法。如图 4所示，三角形单元由节点1、2、3构成，其上的流场参数V1、V2、V3为已知量。Pn为单元内的投影点，d1、d2、d3分别为Pn到三个节点的距离。Pn点流场参数Vn采用反距离权重的插值方法表示为：
|
(6) |
其中权重因子表示为：
|
(7) |
V为任意流场参数，包括速度、压力和密度等，均采用该方法进行插值。
|
| 图 4 反距离权重插值示意图 Fig. 4 Schematic diagram of inverse distance weighted interpolation |
1.4 流线计算的验证算例
采用流场后处理商业软件Tecplot对本文的流线计算方法进行了考核验证。图 5给出了典型升力体外形壁面流场所确定的流线，图中紫色圆点是本文方法计算得到的壁面流线，蓝线为Tecplot软件绘出的壁面流线。图 6为典型高升阻比外形在高超声速流动情况下壁面的流线，同样图中用紫色圆点代表本文方法计算得到的壁面流线，蓝线代表Tecplot软件绘出的壁面流线。通过比较发现：本文计算方法获得的流线与Tecplot软件绘制结果完全一致。
|
| 图 5 典型升力体外形壁面流线图(Ma=8, α=10°) Fig. 5 Surface streamlines on a typical lifting body configuration (Ma=8, α=10°) |
|
| 图 6 高升阻比外形高超声速流动状态下壁面流线图(Ma=20, α=10°, β=5°) Fig. 6 Surface streamlines on a high lift-to-drag ratio configuration at hypersonic flow condition (Ma=20, α=10°, β=5°) |
2 基于当地流场参数的参考温度法 2.1 参考温度法
参考温度法被认为是高超声速高温气体流动下黏性气动力预测的一种较成熟的近似工程算法[7]。在高超声速飞行器选型设计中，该方法得到了广泛应用[11-13]。该方法基于不可压平板边界层流动理论，考虑了可压缩修正和温度影响。
首先讨论层流黏性摩擦力系数的参考温度计算方法。对于不可压平板层流，Blasius解得到当地摩擦力系数为：
|
(8) |
式中，雷诺数
参考温度法引入可压缩修正，可以将当地摩擦力系数写成：
|
(9) |
式中雷诺数
ρ*和μ*为边界层内参考温度T*下的密度和黏性系数。参考温度为边界层外缘马赫数Me和壁面温度与边界层外缘温度比值
|
(10) |
假设边界层内气体参数仍然满足完全气体状态方程，参考温度T*下的密度ρ*可表示为：
|
(11) |
根据边界层内压力的分布特点，通常取压力p*≈pe，于是
|
(12) |
黏性系数μ*由温度T*唯一确定，采用粘温指数函数关系式进行计算[11]，即：
|
(13) |
通常取ω=0.75。
对于湍流情况，也采用与层流参考温度法相同的处理。根据平板不可压湍流常用的摩擦力系数关系式：
|
(14) |
考虑可压缩修正，可以得到：
|
(15) |
式中雷诺数Rex*的计算与层流情况完全相同。
对于层、湍流同时存在的情况，引入转捩雷诺数来判定流动是层流还是湍流，并分段计算黏性力。转捩雷诺数可用如下的近似公式计算[11]，
|
(16) |
其中Rext为转捩雷诺数。当边界层外缘雷诺数Reex>Rext时认为流动为湍流，摩擦力系数采用公式(15) 计算；而当Reex≤Rext时，认为流动为层流，摩擦力系数采用公式(9) 进行计算。
2.2 黏性力计算通常平板流动边界层外缘的物理量近似认为等于自由来流的物理量。但是对于复杂外形的高超声速流动，由于可压缩气流的复杂作用(激波压缩和气流膨胀等)，边界层外缘流场物理量与自由来流的值差别很大，因此如果用自由来流的参数作为边界层外缘的流动参数将使当地摩擦力系数的计算结果产生很大的偏差。采用无黏Euler方程CFD方法可以快速得到无黏流场的特性，将其壁面流动参数作为对应当地边界层外缘的流动参数则能够比较准确的反映流动真实特性。
全体表面的黏性作用力需要通过对所有壁面单元的摩擦力进行求和获得，即:
|
(17) |
其中：Fv为总的黏性作用力；
fi为壁面单元i的黏性摩擦力；
cfi为壁面单元i的黏性摩擦力系数；
QAi为壁面单元i的当地动压；
Ai为壁面单元i的面积。
单元i的黏性摩擦力系数cfi采用基于当地流场参数的参考温度法进行计算，见公式(9) 和(15)，即:
|
(18) |
采用当地位置流线长度计算雷诺数，即:
|
(19) |
式中s为从头部附近发出的到达单元i的流线总长度(图 7)，由第1节中的流线计算方法获得。
|
| 图 7 飞行器表面总的黏性作用力计算示意图 Fig. 7 Schematic diagram of total skin viscous forces computing on the vehicle surface |
单元i的当地动压QAi为参考温度T*下的动压：
|
(20) |
于是单元i的黏性摩擦力可以表示为：
|
(21) |
该摩擦力fi的方向是流线的切向，也就是当地壁面速度方向。令切向单位矢量为τi=(τx, τy, τz)，见图 7。于是单元i的摩擦力fi在直角坐标系(x, y, z) 中的三个分量为：
|
(22) |
于是飞行器表面总的黏性作用力在直角坐标系的三个分量可以表示为：
|
(23) |
无量纲化的黏性作用力系数表示为：
|
(24) |
其中ρ∞、u∞分别为自由来流的密度和速度，AR为参考面积。
3 典型算例验证为考察上述黏性力快速计算方法的准确性，针对高升阻比外形(图 6) 的高空高超声速流动情况进行了研究。其中以Navier-Stokes方程CFD数值计算方法[14-15]获得的完全气体模型下的黏性力作为准确值，并采用基于自由来流条件的经典平板边界层参考温度法和本文快速计算方法分别进行了黏性力的计算。本文的研究选用了马赫数10和20两个速度情况，高度也选取了50km和60km两组情况，迎角最大值达到20°。通过对Navier-Stokes方程CFD数值计算结果的分析表明，各计算状态下飞行器表面未发生明显流动分离现象。
图 8给出了三种方法所获得的黏性轴向力分量结果的比较。可以看到，在马赫数10情况下，经典参考温度法和本文计算方法所获得的结果相差不大；马赫数20情况下，经典参考温度法计算结果与准确值存在较大偏差，而本文计算方法获得的结果与准确值更接近；在50km和60km两种不同高度情况下，上述规律完全一致。
|
| 图 8 高升阻比外形轴向力黏性分量计算结果比较 Fig. 8 Computation results comparison about skin viscous parts of axial forces on a high lift-to-drag ratio configuration |
同时我们发现在所研究的马赫数和迎角范围内，马赫数越高，迎角越大，经典参考温度法计算结果偏差越大，而本文计算方法与准确值更接近。这是因为马赫数越高，迎角越大，自由来流受到飞行器影响就越大，边界层外缘实际流动参数与自由来流之间的差异则更大；采用本文方法以Euler方程CFD方法获得的准确流动参数作为边界层外缘参数，相比于经典参考温度法以自由来流作为边界层外缘参数，能够更精确的反映控制黏性底层流动的外围流动特性，因此本文方法计算结果相比于经典参考温度法会更加准确。
总之，相比于经典参考温度法，本文计算方法的结果与Navier-Stokes方程CFD数值计算结果的一致性更好，存在约10~20%的偏差，能够满足工程初步设计的要求。
另一方面，从计算效率来看，对于单个计算状态，采用Navier-Stokes方程CFD方法所需要的计算时间高出Euler方程CFD方法计算时间一个量级，而本文基于Euler方程CFD方法计算获得的无黏流场所建立的黏性力快速计算方法在个人台式计算机单个CPU上运行所需的计算时间通常只有几秒或几十秒，相比于CFD方法的计算时间可以忽略不计。因此本文计算方法与Navier-Stokes方程CFD方法的效率之比，相当于Euler方程CFD方法与Navier-Stokes方程CFD方法的效率之比，即达到10倍的量级。可见本文方法在计算效率方面具有明显的优势。
4 结论本文从工程应用的角度出发，基于无黏Euler方程CFD方法获得的壁面流场参数，建立了一种高超声速飞行器壁面黏性力的快速工程计算方法。所提出的流场自适应变时间步长流线计算方法大大提高了流线计算效率；建立的基于当地流场参数的参考温度法，对壁面黏性力的计算相比于经典参考温度法具有更高的精确性，并得到了Navier-Stokes方程CFD方法的验证，能够满足工程设计的要求。
由于壁面流线建立的前提是流动在壁面没有发生分离，因此本文的计算方法在理论上仅仅适用于没有发生流动分离的情况。但飞行器高超声速飞行环境下，在一定迎角范围内通常不会发生比较明显的分离流动；而且在工程设计中也需要尽量避免大迎角高超声速飞行情况。因而，对于高超声速飞行器工程设计来讲，本文建立的高超声速壁面黏性力快速计算方法仍然具有良好的准确性和实用性。
| [1] |
Yan C, Yu J, et al. On the achievements and prospects for the methods of computational fluid dynamics[J].
Advances In Mechanics, 2011, 41(5):562–589.
(in Chinese) 阎超, 于剑, 等. CFD模拟方法的发展成就与展望[J]. 力学进展, 2011, 41(5) : 562–589. |
| [2] | Zhu Z Q, et al. Application computational fluid dynamics[M]. Beijing: Press of Beijing University of Aeronautics and Astronautics, 1998. |
| [3] |
Wang R H, Yin G L, et al. Fast CFD tools for civil aircraft conceptual design and optimization use[J].
Aircraft Design, 2012, 32(5):31–35.
(in Chinese) 王如华, 尹贵鲁, 等. 快速CFD计算工具在民机概念优化设计中的应用[J]. 飞机设计, 2012, 32(5) : 31–35. |
| [4] |
Xu R F, Deng Y J, et al. Aerodynamic optimization design and its requirement to CFD[J].
Aeronautical Science & Technology, 2011(2):50–52.
(in Chinese) 许瑞飞, 邓一菊, 等. 气动优化设计及其对CFD的需求[J]. 航空科学技术, 2011(2) : 50–52. |
| [5] | Van Driest E R. On turbulent flow near a wall[J]. Journal of the Aeronautical Sciences, 1956, 23(11):1007–1011. DOI:10.2514/8.3713 |
| [6] | Fleeman E L. Tactical missile design, second edition[M]. Virginia: AIAA Inc., 2006. |
| [7] | Anderson J D. Hypersonic and high temperature gas dynamics, second edition[M]. Virginia: AIAA Inc., 2006. |
| [8] | Gnoffo P A. Computational fluid dynamics technology for hypersonic applications[R]. NASA 20040013407, 2004. |
| [9] | Legendre R. Streamline in a steady flow[R]. NASA ADA395523, 1966. |
| [10] | Diachin D P, Herzog J. Analytic streamline calculations on linear tetrahedral[R]. AIAA-97-1975:733-742, 1997. |
| [11] | Corda C, Anderson J D. Viscous optimized hypersonic waveriders designed from axisymmetric flow fields[R]. AIAA-88-0369, 1988. |
| [12] | Frederick F, Terry C, et al. The development of Waveriders from axisymmetric flowfield[R]. AIAA 2007-847. |
| [13] |
Zhang D, Tang S, et al. Engineering calculation method of viscous force for air-breathing hypersonic vehicle[J].
Journal of Solid Rocket Technology, 2013, 36(3):291–295.
(in Chinese) 张栋, 唐硕, 等. 吸气式高超声速飞行器黏性力工程计算方法[J]. 固体火箭技术, 2013, 36(3) : 291–295. |
| [14] |
Gong A L, Liu Z, et al. Numerical study on hypersonic double-cone separated flow[J].
Acta Aerodynamica Sinica, 2014, 32(6):767–771.
(in Chinese) 龚安龙, 刘周, 等. 高超声速激波/边界层干扰流动数值模拟研究[J]. 空气动力学学报, 2014, 32(6) : 767–771. |
| [15] | Liu Z. Delayed detached eddy simulation for static and forced oscillating airfoil at high angle of attack[R]. IAC 2013, 2013. |