王昭力，曾涛，周志宏,*，陈宇

1. 四川大学 建筑与环境学院，成都 610065

2. 宜宾四川大学产业技术研究院，宜宾 644000

飞机穿过含有过冷水滴的云层时，水滴撞击到飞机表面可能凝结成冰，结冰会改变飞机升力部件的外形并影响其气动性能，导致飞机的操纵性和稳定性恶化；
其他部位结冰也可能导致传感器失灵、操纵面卡死、发动机失火等不良影响，严重时直接导致空难事故[1-2]。飞机结冰研究主要有飞行试验、风洞试验和数值模拟3种研究手段，其中数值模拟具有成本低、效率高等优点，在飞行器防除冰系统设计、适航取证、空难事故分析等阶段被广泛应用。飞机结冰是水滴在空气绕流作用下，撞击到飞机表面发生凝结的现象，这一物理过程的数值模拟主要包含水滴轨迹及撞击到物体表面的收集率、传热相变及水膜溢流、冰形增长3个主要的计算模块，其中，水滴收集率计算是整个结冰计算的基础，是其他模块的输入依据。同时，水滴收集率也是防冰系统设计的重要输入依据。

水滴收集率计算主要有基于水滴粒子的拉格朗日方法和基于两相流观点的欧拉方法。最先发展起来的水滴撞击特性的计算方法是拉格朗日轨迹法，在远场释放一定数量水滴粒子，通过求解水滴受力方程，获得这些水滴的运动轨迹，通过统计撞在物面上的水滴数量或撞击面积比等方式，获得飞机表面的水滴收集率，国内杨倩[3]、易贤[4]等在这些方面做了许多工作。

水滴轨迹计算的欧拉两相流法是将空气和水滴视为两相流，水滴相分布在空气相中，通过在全流域内求解水滴相的控制方程得到水含量的分布以及水滴的速度，再依据物面处的水含量及法向速度计算水滴收集率，相比拉格朗日法，欧拉法在处理复杂外形水滴撞击特性计算时更有优势。欧拉两相流计算方法的问题在于，控制方程的高度非线性容易带来数值振荡。Bourgault等[5]建立了两相流动的控制方程并使用有限单元法进行求解，后又在此基础上添加了间断捕捉项来解决欧拉法数值不稳定的问题，并将此算法植入到FENSAP中[6]。Iuliano等[7]对多段翼的水滴流场进行了计算，并对比了欧拉法与拉格朗日法之间的区别以及空气黏性对水滴撞击特性的影响。Tong等[8]分析了欧拉法数值不稳定的原因，在原始的水滴控制方程基础上添加了扩散项用来消除迭代过程中的奇异性，该扩散项能根据网格粗细和水含量梯度自动调整。国内张大林等[9]在欧拉法中对流项使用MUSCL格式计算了二维NACA0012翼型的水滴收集率，林贵平等[10]在水滴的质量守恒方程中添加自适应的扩散项来提高使用对流项高阶离散格式时的稳定性，易贤等[11]使用欧拉法发展了适合于三维、复杂构型的水滴收集率计算程序，其对流项的离散是采用迎风插值和线性插值相结合的方法。招启军等[12]为解决尾流等区域的密度脉冲现象引起的稳定性和收敛性问题，提出并建立了遮蔽区扩散模型。陈海昕等基于欧拉法，发展了一套结冰模拟软件[13]。

由于一阶格式的精度不够，而单纯的高阶格式又容易引起数值振荡和密度脉冲，目前针对空气流场对流离散格式的研究，已经有较多的研究[14-17]，但绝大部分未应用于水滴流场的数值求解中。水滴流场和空气流场的求解有所不同，水滴流场的控制方程相比常规的空气流场没有黏性项，因此其对数值振荡或者奇异值更加敏感；
此外，针对空气流场的对流项离散格式的研究大多是集中于超高速工况，由于超高速会使空气流场产生激波，但激波前后流动参数形成的间断性质与水滴流场中遮蔽区的界面间断性质不同。因此，研究水滴流场的对流项离散格式是有意义的。

本文在不添加扩散项的前提下，通过在二阶精度格式中引入系列通量限制器来解决收敛性的问题，以发展飞机结冰计算中水滴流场的高阶精度计算格式。

1.1 水滴流场控制方程

针对飞机结冰过程中水滴-空气两相流欧拉法的建立，由于云雾中水滴含量很小，通常有如下假设：

1) 将水滴视为是分布在空气中的连续相，采用液态水含量(Liquid Water Content, LWC)表示单位体积中的水质量，单位为g/m3。引入水滴体积分数α，定义为控制单元内水滴体积占总体积的比例，自由来流的α可采用如下公式计算：

(1)

式中：ρ为水的密度，运动过程中水滴密度不发生改变。

2) 常规水滴条件(水滴平均粒径<40 μm)下，水滴为球形，单个水滴的尺寸用水滴当量直径deq描述。

3) 水滴在运动过程中不会发生变形、融合和破碎，且撞击物面后不发生飞溅。

4) 水滴的动量方程中不考虑黏性和压力项。

5) 由于水滴体积占比很小，往往在10-6级别，因此，可近似假设水滴和空气是单向耦合的，即只考虑水滴受到的空气阻力，水滴对空气流场的影响可以忽略，且水滴与空气不发生任何能量或质量的交换。

6) 作用在水滴上的力只有阻力和重力，且阻力是稳态的。

基于上述假设，可以获得如下微分形式的水滴流场的控制方程：

(2)

(3)

使用中心有限体积方法，对积分形式的水滴控制方程在控制体内进行积分，并使用高斯定理，写成向量的形式，可得到：

(4)

式中：

(5)

式中：ua和u分别为空气和水滴的速度向量；
F为水滴受到的体积力(只有重力)；
u、v、w和ua、va、wa分别为水滴速度向量u和空气速度向量ua的分量；
Fx、Fy、Fz为体积力F的分量；
τm为动量反应时间：

(6)

f为阻力函数：

(7)

由于常见的飞机飞行时的雷诺数超出了斯托克斯公式计算水滴黏性阻力的理论允许范围，因此采用修正的阻力函数[11]：

(8)

Rer为相对雷诺数：

(9)

式中：μa为空气的动力黏性系数。

1.2 空气流场计算

欧拉法的求解是将空气和水滴视为两相流，且只考虑空气相对水滴相的单向耦合作用。因此在计算水滴流场前需先计算空气流场。

本文的空气流场采用中国空气动力研究与发展中心计算空气动力研究所开发的NNW-PHengLEI软件进行计算。计算采用SSTk-ω模型，通量分裂采用Roe格式，对流项通量限制器使用的是Minmod限制器。

1.3 水滴流场控制方程的TVD格式

水滴流场的控制方程是一个带有源项的双曲型非线性对流方程，因为不具有压力梯度和扩散项来及时抹平水含量的高梯度，导致了数值不稳定性，因此其求解的关键在对流项的离散，需要同时考虑数值解的准确性和稳定性。

将水滴流场控制方程的积分形式半离散化：

(10)

(FC·S)i,j+1/2,k-(FC·S)i,j-1/2,k+

(FC·S)i,j,k+1/2-(FC·S)i,j,k-1/2

(11)

式中：S为控制体表面的面积矢量；
FC为控制体表面的对流矢量。

FC先使用通量矢量分裂(Flux Vector Splitting, FVS)方法分裂成非负和非正两组：

FC(Q)=FC+(QL)+FC-(QR)

(12)

式中：QL和QR分别为控制体表面左右两侧的守恒变量[ααuαvαw]T。

由于水滴流场控制方程的Jacobian系数矩阵有4重特征值Vc，Vc为逆变速度，是流场速度与网格线切向矢量的点积，ξ、η、ζ分别为曲线坐标系下的3个方向：

Vc=ukx+vky+wkzk=ξ,η,ζ

(13)

(14)

控制体表面的值是根据控制体上的均值插值得到的，不同的插值方式决定了对流项的离散格式的空间精度和收敛性。

1.3.1 控制体表面的对流通量插值

对于任意的物理量q，一阶迎风格式(First-Order Upwinding, FOU)可以写为

qi+1/2=qi

(15)

普通的低阶格式，包括一阶迎风格式、混合格式、指数格式等拥有良好的稳定性和有界性，但精度较低，往往导致过度的数值扩散。而单纯的二阶格式会使流场的间断附近产生非物理振荡，有可能导致迭代的发散。为了提高精度，研究者们提出了一系列线性对流格式，这些线性对流格式可以用van Leer提出的“k-格式”来概括：

qi+1/2=qi+

(16)

当k取不同的值时能够获得不同的高阶对流格式，这些格式能在一个维度上使用5点模块来计算，且在计算线性对流方程时都是线性的。但是“k-格式”的所有成员都不具备有界性，Godunov定理阐明了任何截断误差高于一阶的线性对流格式都不会是单调的[18]，因此在计算有间断或梯度很大的情况时会出现严重的非物理数值振荡，从而导致计算的发散。

为了解决低阶格式的低精度和高阶格式的不稳定的问题，需要构造非线性的对流格式。最常用和有效的一种构造非线性对流格式的方法就是“通量限制器”的使用。

重构过程的插值可以写成如下“标准半离散通量限制形式”，在均匀网格下：

(17)

式中：ri+1/2为梯度率，表示中心梯度和迎风梯度的比值，本质上是针对变量平滑度的局部测量：

(18)

以ξ方向为例，在非均匀结构网格中水滴流场控制体表面左右两侧的守恒变量可以写为

(19)

(20)

(21)

式中：ψ(r)即为通量限制器函数，不同的通量限制器函数决定了对流项离散格式的精度和收敛性。

1.3.2 水滴流场发散原因

对于不使用通量限制器的高阶格式计算水滴流场时，导致发散的原因一般有两个：

1) 水滴流场中，虽然不存在空气流场中的激波现象，但是由于物面遮挡效应所导致的无水区边缘的α的梯度很大，使用高阶精度格式计算时也会出现数值振荡。

2) 由于欧拉法中水滴速度的单值性引发水滴流动交汇处出现奇异值,如图1所示。在物体后部一般会有水滴交汇的现象出现，对于拉格朗日法，每个水滴是独立运动的，其轨迹可以发生交叉和重叠，此处α的值为两股交汇水流的α之和；
但对于欧拉法，由于水滴被视为是连续相，流线不能出现交叉或重叠，两股水滴流动分别从上至下和从下至上汇聚成一股向同一方向流动，导致了此处的α异常增大，若离散格式缺乏耗散项及时将此处的异常流动抹平，则极其容易导致计算的发散[8]。

图1 水流交汇示意图

因此需要引入限制准则，并在对流项的离散格式中添加通量限制器，以此来抑制水滴流场间断处的波动和抹平欧拉法中速度单值性带来的奇异值。

1.3.3 TVD限制准则

总变差减小(Total-Variation Diminishing，TVD)的概念最早是由Harten提出的[19]。总变差是度量函数波动状况的一种指标，对于离散的数值解来说，第n个时间步的总变差的定义为

(22)

而TVD的特征为

TV(qn+1)≤TV(qn)

(23)

总变差能够衡量数值解的振荡情况，满足TVD的格式能够保证数值解的单调和有界。Sweby[20]在显式离散的背景下推导了用通量限制器函数表达的Sweby’s TVD限制准则：

(24)

Sweby’s TVD限制准则不依赖于CFL数，对于求解半离散对流方程的恒定状态解来说，不依赖CFL数的限制准则更加适用。而在空气流场恒定的情况下，水滴流场就是定常的，其求解过程是使用“伪时间步”推进的。“k-格式”系列所包括的区域与Sweby’s TVD限制区域重合的部分，既可以保证总残差减小又能保证空间二阶精度。

国内外有许多的完全或部分满足上述限制准则的格式被相继开发出来，可以统称为TVD格式，这些格式大致可以分为两类：逐段线性(Piecewise-Linear)格式和平滑(Gradually-switching Smooth)格式，本文将5种经典的通量限制器应用于水滴流场控制方程的对流项离散，这些通量限制器的ψ-r图像如图2所示。

图2 5种通量限制器的ψ-r图像

1) Minmod

ψ(r)=max(0,min(r,1))

(25)

Minmod是一种最简单的限制器，使用时，实际上是只用迎风梯度或中心梯度来修正一阶迎风量。当在极值处时r≤0，离散格式直接降为一阶迎风格式，而远离极值处时，离散格式会在中心差分和二阶迎风格式中切换。张涵信[16]提出的NND格式就是使用了此限制器，并阐述了中心差分与二阶迎风格式的组合使用何以能限制间断(梯度陡峭)处上下游的数值振荡。

2) MUSCL

(26)

MUSCL[20-21]限制器是在Sweby’s TVD限制区域内使用Fromm格式，因此在光滑区域具有Fromm格式的精度和收敛性。

3) Superbee

ψ(r)=max(0,min(2r,1),min(r,2))

(27)

Superbee[22]限制器在ψ-r图中是二阶Sweby’s TVD限制区域的上界，此格式是被设计用来详细刻画流场间断的，但是由于在高曲率区域大量使用一阶顺风格式导致了光滑区域的过度压缩[20]。和Minmod限制器一样，Superbee限制器在光滑区域也是使用中心差分与二阶迎风格式的组合，不过判断条件与其相反。

4) Harmonic

(28)

Harmonic[23-24]格式实质上是用中心梯度和迎风梯度的调和平均值来修正一阶迎风项，其也是基于Fromm格式的斜率限制器，但只在r=1处具有和Fromm格式相同的斜率。

5) Van Albada

(29)

Van Albada限制器在整个定义域上都是光滑、连续且可微的，这为离散格式带来了良好的收敛性。与其他限制器不同的是，此限制器在r≤0时并不为0，即完全不使用一阶迎风格式。

将上述通量限制器代入到式(19)、式(20)中即可得到对应的对流项离散格式。

1.4 时间推进格式

对于水滴流场的半离散式(10)，令

(30)

略去下标i,j,k，则式(10)可以写为

(31)

式中：R称为方程的残差。

采用显示的四步Runge-Kutta法来进行时间推进，具体格式如下：

(32)

式中：α1=1/4;α2=1/3;α3=1/2;α4=1;Δt为时间步长。

1.5 边界条件和初始条件

水滴流场的求解采用有限体积法，流场的边界条件通过设置虚网格格心值的方式给出，边界面上的通量基于虚网格格心值插值得到，水滴流场的边界条件分为自由来流边界和物面边界。

自由来流边界上α如1.1节中假设1)所述，速度向量根据攻角设置；
物面边界的设置需要判断撞击区和无水区，u为水滴速度向量，n为物面的外法向量。当u·n≤0时，即为撞击区，此时水滴撞上物面，并被全部捕捉后结冰，计算时将水滴视为穿过此边界，该区域虚网格的α和u均与贴近物面的第1层网格的值保持相同；
当u·n>0时,即为无水区，此时相当于水滴从物面上流出，实际上此区域是被物体遮挡，没有水滴的流入，理论上此边界的α应设置为0，为了避免计算过程中出现数值错误，该区域虚网格的α取物面上第1层网格的α的10-9倍，u与其保持一致。

针对物体之后在迭代过程中会存在奇异值的问题，将初始流场的α设置成一个接近于0的小值可以有效的解决此问题，相当于水滴从自由来流边界渐渐传输到整个计算域，物体之后的无水区一直是保持没有水滴的状态，从根本上避免了水滴流动交汇产生奇异值的问题。

但上述初始条件的计算较慢，本文的初始条件设置如下：水滴流场全场的初始α用工况设置中的LWC计算，初始速度沿攻角方向，速度大小由来流马赫数给出。TVD格式在这种初始条件下物体之后虽然存在密度脉冲，但不会引起计算的发散。

1.6 水滴收集率计算

水滴流场的求解主要是为了获得物面水滴收集率，为结冰或防冰计算提供输入，水滴收集率的定义是物面表面微元的实际水收集量与最大可能的水收集量间的比值，用于衡量物面收集水滴能力的大小，水滴收集率的计算采用如下公式：

(33)

式中：αn为物面附近的水滴体积分数；
α∞为无穷远处的水滴体积分数；
u∞为无穷远处水滴的速度矢量。

2.1 NACA0012翼型

对于NACA0012翼型，计算了如表1所示的4种工况和如表2所示的3种粗细的网格。

图3为Case1工况下一阶迎风格式和5种带不同限制器的TVD格式计算得到的水滴收集率β与实验数据的对比。图中：s为到驻点的距离；
c为机翼的弦长。图4为Case1工况下计算过程中α的收敛曲线，可以看出，相比于一阶迎风格式，具有限制器的TVD格式均能更好的计算水滴收集率。一阶格式的误差主要出现在驻点附近，驻点附近会出现计算的水滴收集率偏大的现象，结合图5中水滴体积分数的分布，可看出，一阶格式的计算结果会在驻点附近产生水含量的集中。

表1 NACA0012翼型计算工况

表2 NACA0012翼型计算所用网格参数

TVD格式的计算结果与实验数据[25]相比，在驻点附近偏小，撞击极限偏小，且有轻微的不对称现象。收集率和撞击极限偏小的原因可能在于实验过程无法保证所有水滴的直径均相同，而是在一个区间内分布；
而不对称的来源可能是网格的不对称和重力的作用效应，但小水滴的重力和空气对水滴的作用力相比是一个小量，小水滴的重力几乎可以忽略不计，因为有限体积方法对流项的离散需要用到每个单元网格的面法向量，所以网格的不对称可能会造成机翼上下区域水含量

图3 Case1的水滴收集率

图4 Case1中α的迭代收敛曲线

图5 Case1的水滴体积分数云图

轻微出入，因此导致水滴收集率的不对称。

图6为Case2工况下网格适应性分析的结果。对于文中所使用的5种TVD格式，中网格和细网格的计算结果基本一致，只在驻点处轻微有一点出入。实际上粗网格的计算结果趋势也与中网格和细网格的一致，只是因为驻点处的点不够密，无法精细的描述出物面驻点处的水滴收集率分布。总体来说TVD格式网格适应性较好，在普通密度的网格中就能计算出良好的水滴收集率分布结果。

图7、图8分别为Case3、Case4下，使用粗细不同的两套网格计算得到的水滴收集率，网格的参数如表2所示，计算结果与Lewice和DROP3D软件计算的[26]进行对比。不论是粗网格还是细网格，一阶迎风格式均在驻点处有水含量过度集中导致的水滴收集率偏大。当水滴直径较小时(20 μm)，Minmod限制器也出现了驻点处水滴收集率偏大的情况。对比粗网格的计算结果，使用细网格时，一阶格式与TVD格式的计算结果更加接近，几乎只在驻点附近不同。说明当网格比较稀疏时，TVD格式能够较好的提升计算水滴收集系数的精度。使用细网格计算时，Case4状态下，Superbee格式计算水滴流场的迭代过程发生了发散，主要原因在是翼型遮挡部位的无水区边缘产生了无法抹平的数值波动。

图6 Case2中用TVD格式在不同尺寸网格下计算的水滴收集系数

图7 Case3的水滴收集系数

图8 Case4的水滴收集系数

2.2 NACA23012翼型

NACA23012翼型的计算工况如表3所示。

表3 NACA23012翼型计算工况

从图9水滴收集率的计算结果来看，一阶格式的计算结果在驻点处还是偏大，所采用的5种TVD格式的计算结果相差无几。TVD格式间的差别主要体现在对间断的分辨率的差别上，由于水滴收集率只涉及物面前缘的流场，这一部分并没有流动的间断，因此所采用的5种TVD格式的计算结果非常接近。

对于Case5(图9和图10)，除了撞击极限外，一阶格式的计算结果偏大，而TVD格式的计算结果与实验数据基本一致。撞击极限偏小的原因可能是由于本文程序是按照单一的水滴尺寸计算的，而实验和Lewice中的水滴并不是单一尺寸，而是呈现一定的分布的，其中包含有一定量大尺度粒径的水滴，这部分水滴撞击会导致撞击极限更大；
对于Case8(图11和图12)，一阶格式在驻点处收集率偏大，其余部分与TVD格式基本一致。而与实验结果相比，不论是一阶格式还是TVD格式，均在驻点处偏大，且撞击极限偏小，但总体来说，本文程序的计算结果较Lewice[26]更接近实验数据[27]。

图9 Case5的水滴收集系数

图10 Case5中α的迭代收敛曲线

图11 Case8的水滴收集系数

图12 Case8中α的迭代收敛曲线

从α的收敛曲线上看，残差变化平稳后，不同格式的残差从大到小的顺序为Superbee、MUSCL、Van Albada、Minmod、Harmonic、一阶迎风格式。对于Case8，使用Superbee格式计算时发生了发散，发散的原因是在无水区边缘产生了数值波动且无法及时抹平。

图13为使用Minmod限制器时在Case7下计算得到水滴体积分数云图，可以看到尾部有明显的高密度脉冲。图14是在其中取y/c=0.036截面得到的α随横坐标x/c的变化趋势，可以看出，单值性导致的物体后方的高密度区域会随着deq的增加而愈发趋于明显和陡峭。

图15是使用不同的对流项离散格式后取截面后得到的α随横坐标x/c的变化，使用TVD格

图13 Case7下使用Minmod格式的水滴体积分数云图

图14 Minmod格式下不同粒径尾部大密度区域对比

图15 Case7下不同限制器尾部大密度区域对比

式时在物体后方产生的高密度区域要比使用一阶格式时更加的明显，陡峭程度从大至小依次为：Superbee、MUSCL、Van Albada、Minmod、Harmonic、一阶迎风格式，这与残差从大至小的顺序一致。由此可知，在残差变化平稳后，影响水滴流场迭代求解的因素就是物体之后的密度脉冲。

图16为使用另一种初始条件，即将初始α设置成一个接近于0的小值时，在Case7下使用Minmod限制器计算得到的水滴体积分数云图。可以看到此种初始条件的设置可以有效地解决物体之后密度脉冲的问题。

2.3 30P30N三段翼

30P30N三段翼的计算工况如表4所示。图17为30P30N三段翼的外形示意图，机翼由前缘缝翼、翼身和后缘襟翼3部分组成。机翼弦长c的定义为当前缘缝翼和后缘襟翼均未偏转时整体的长度，在此算例中，前缘缝翼和后缘襟翼的偏转角度均为30°。前缘缝翼的上下两端容易发生边界层分离和转捩。

图16 在初始水滴体积分数接近于0时Case7下使用Minmod格式的水滴体积分数云图

Table 4 Calculation conditions for 30P30N three-element airfoil

图17 三段翼30P30N的外形示意图

图18～图20分别为前缘缝翼、翼身和后缘襟翼上的水滴收集系数分布曲线。使用Van Albada格式计算时，在前缘缝翼的上端出现了不可抹平的奇异值。对于其他几种计算格式，在前缘缝翼上，TVD格式在驻点附近和上端撞击极限处跟一阶迎风格式有所出入；
在翼身上，5种格式的计算结果基本一致；
在后缘襟翼上，4种TVD格式的计算结果几乎没有出入，而一阶迎风格式相比下水滴收集率的分布整体向左侧偏移。

图21为Case9中的收敛曲线，平均残差从大到小依次为Superbee、MUSCL、Minmod、Harmonic、一阶迎风格式，顺序与之前计算的几种工况一致。TVD格式中的Van Albada格式无法处理此算例中边界层分离的情况，而其他的格式对复杂构型的适应性较好。

图18 Case9中前缘缝翼的水滴收集系数

图19 Case9中翼身的水滴收集系数

图20 Case9中后缘襟翼的水滴收集系数

图21 Case9中α的迭代收敛曲线

将TVD格式应用到水滴流场求解的对流项离散中，通过计算分析，有如下结论：

1) 基于水滴收集率计算和实验数据的对比，TVD格式较一阶迎风格式的精度有较好的提升，特别是在驻点附近，但所采用的5种TVD格式计算的精度无明显差别。

2) TVD格式有较好的网格适应性，在较稀网格条件下，两种格式的计算精度差距更大。

3) 在欧拉法水滴流场迭代计算，残差的主要来源之一是物体之后的密度脉冲，随着水滴粒径的增大，高密度区域范围和密度都随之增加。将初始条件设为接近于0的小值时可以避免密度脉冲。若使用本文中的的初始条件设置方式，不同的TVD格式对此区域的抹平能力不同，从强至弱依次为：Harmonic、Minmod、Van Albada、MUSCL、Superbee。

4) Superbee格式的稳定性较差，不仅体现在物体尾部高密度区域抹平能力较差，而且在部分工况下，无法处理水滴流场间断导致的数值振荡。

5) 针对复杂构型，Van Albada格式在计算过程中发生了发散，而其他格式都其表现出良好的适应性。

6) 在不修改水滴流场控制方程的前提下，探索用通量限制器构造的TVD格式求解水滴流场的可能性，其计算结果也证明此类格式是有效的。

猜你喜欢 驻点对流水滴 齐口裂腹鱼集群行为对流态的响应农业工程学报(2022年11期)2022-08-22利用水滴来发电红领巾·探索(2020年6期)2020-10-26水滴轮的日常拆解与保养办法垂钓(2019年2期)2019-09-10完全催化壁驻点高超声速流动加热地面模拟方法研究空气动力学学报(2018年2期)2018-04-25超临界压力RP-3壁面结焦对流阻的影响北京航空航天大学学报(2016年6期)2016-11-16基于ANSYS的自然对流换热系数计算方法研究现代计算机(2016年11期)2016-02-28利用远教站点，落实驻点干部带学党的生活·党员电教与远程教育(2016年3期)2016-02-262300名干部进村“串户”办实事源流(2015年8期)2015-09-16“一线为民”的庐阳探索决策(2015年2期)2015-09-10航天器相对运动水滴型悬停轨道研究航天器工程(2014年4期)2014-03-11